From e1848e20f9195ab25e56bd4dc9cb72f18a7f8cdb Mon Sep 17 00:00:00 2001
From: Ryan Cabell <rcabell@ucar.edu>
Date: Tue, 12 Nov 2024 16:46:00 -0700
Subject: [PATCH] Bypass nudging and set diversion destination to observed
 value

---
 src/Routing/Diversions/module_diversions.F | 35 +++++++++-----------
 src/Routing/module_channel_routing.F       | 37 ++++++++--------------
 2 files changed, 29 insertions(+), 43 deletions(-)

diff --git a/src/Routing/Diversions/module_diversions.F b/src/Routing/Diversions/module_diversions.F
index 35c2f7998..11e2bad9c 100644
--- a/src/Routing/Diversions/module_diversions.F
+++ b/src/Routing/Diversions/module_diversions.F
@@ -106,16 +106,16 @@ subroutine init_diversions(diversions_file, timeslice_path)
 
     end subroutine
 
-    subroutine calculate_diversion(src_link_in, dst_link_out, qlink_src_in, diversion_quantity_out)
+    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_link_out
+        ! integer(kind=int64), intent(out) :: dst_out
         real, intent(in) :: qlink_src_in
-        real, intent(out) :: diversion_quantity_out
+        real, intent(out) :: diversion_quantity_out, diversion_quantity_in
 
         integer :: i
 
-        dst_link_out = -99
         diversion_quantity_out = 0
+        diversion_quantity_in = 0
 
         ! bail if we're inactive
         if (.not. diversions_active) return
@@ -130,11 +130,21 @@ subroutine calculate_diversion(src_link_in, dst_link_out, qlink_src_in, diversio
                     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
+                    ! 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)
@@ -162,21 +172,6 @@ subroutine gage_assisted_diversion(src_link, diversion, qlink_src, div_gage_flow
             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
diff --git a/src/Routing/module_channel_routing.F b/src/Routing/module_channel_routing.F
index e30b522de..ce8b594bf 100644
--- a/src/Routing/module_channel_routing.F
+++ b/src/Routing/module_channel_routing.F
@@ -1713,8 +1713,7 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT,  &
        character(len=256), allocatable,dimension(:) :: tmpAssimilatedSourceFile
 
        ! diversions
-       integer(kind=int64) :: div_dest
-       real :: div_qty
+       real :: div_src, div_dst
        character(*), parameter :: free = '(*(g0,1x))'
 
 #ifdef MPP_LAND
@@ -1968,35 +1967,27 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT,  &
 
          ! HANDLE DIVERSIONS
 
-         call calculate_diversion(LINKID(k), div_dest, Quc, div_qty)
-         if (div_qty /= 0) then
+         call calculate_diversion(LINKID(k), Quc, div_src, div_dst)
+
+         if (div_src /= 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."
+            print free, "DEBUG: diverting", div_src, "of", Quc, "from link id =", LINKID(k), "on processor", my_id
+            if (div_src > Quc) &
+               print free, "DEBUG WARNING: diverted flow (", div_src, ") exceeds total flow, zeroing."
 #endif
+            Quc = max(0.0, Quc - div_src)
+            Qup = max(0.0, Qup - div_src)
+         end if
 
-            ! 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
+         if (div_dst /= 0) then
+            ! apply observed value to downstream
 #ifdef HYDRO_D
-                     print free, "Found diversion destination ID", LINKID(kk),"in local processor at index", kk, ", replacing", tmpQLINK(kk,2), "with", div_qty
+            print free, "DEBUG: diverting", div_dst, "to link id =", LINKID(k), "on processor", my_id
 #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
+            tmpQLINK(k,2) = div_dst
          end if
 
-
 #ifdef WRF_HYDRO_NUDGING
          call nudge_apply_upstream_muskingumCunge( Qup,  Quc,  nudge(k),  k )
 #endif