Apply diversion to both Qup and Quc

CMakeLists.txt

@@ -205,13 +205,13 @@ if (CMAKE_Fortran_COMPILER_ID MATCHES "GNU.*")
     if (CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 9)
       set (CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fallow-argument-mismatch")
     endif()
-    set (CMAKE_Fortran_FLAGS_RELEASE "-O2")
+    set (CMAKE_Fortran_FLAGS_RELEASE "-O2 -march=native")
     set (CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g -fbacktrace")
 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/OrchestratorLayer/config.f90

@@ -170,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

src/Routing/Diversions/module_diversions.F

@@ -106,10 +106,10 @@ subroutine init_diversions(diversions_file, timeslice_path)
 
     end subroutine
 
-    subroutine calculate_diversion(src_link, dst_link, qlink_src, qlink_dst)
+    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
+        real, intent(inout) :: qlink_src_previous, qlink_src_current
         real, intent(inout) :: qlink_dst
 
         integer :: i
@@ -126,17 +126,17 @@ subroutine calculate_diversion(src_link, dst_link, qlink_src, qlink_dst)
                 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)
+                    call gage_assisted_diversion(src_link, diversions(i), qlink_src_previous, qlink_src_current)
                 end if
             end if
         end do
 
     end subroutine
 
-    subroutine gage_assisted_diversion(src_link, diversion, qlink_src)
+    subroutine gage_assisted_diversion(src_link, diversion, qlink_src_prev, qlink_src_curr)
         integer(kind=int64), intent(in) :: src_link
         type(diversion_t), intent(in) :: diversion
-        real, intent(inout) :: qlink_src
+        real, intent(inout) :: qlink_src_prev, qlink_src_curr
 
         real :: div_gage_flow, fraction
 
@@ -155,15 +155,32 @@ subroutine gage_assisted_diversion(src_link, diversion, qlink_src)
             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
+            div_gage_flow = qlink_src_curr * 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, "from link id =", src_link
+            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
-        qlink_src = qlink_src - div_gage_flow
+        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
     end subroutine
 
 end module

src/Routing/Diversions/module_diversions_timeslice.F

@@ -9,6 +9,8 @@ module module_diversions_timeslice
     integer, parameter :: PATH_MAX = 4096
     type(glob_t) :: timeslice_files
 
+    character(*), parameter :: free = '(*(g0,1x))'
+
     contains
 
     integer function init_timeslices(timeslice_path) result(ierr)
@@ -51,6 +53,9 @@ real function get_flow_for_gage(gage) result(flow)
 
             found = findloc(gage_ids, gage)
             if (found(1) /= 0) then
+#ifdef HYDRO_D
+                print free, "DEBUG: Reading diversion discharge for gage " // trim(adjustl(gage)) // " from " // trim(timeslice_files%filenames(i))
+#endif
                 ierr = ierr + nf90_inq_varid(ncid, 'discharge', varid)
                 ierr = ierr + nf90_get_var(ncid, varid, discharge, start=found, count=(/1/))
                 if (ierr /= 0) call hydro_stop("Error occurred reading gage data from " // trim(timeslice_files%filenames(i)))
@@ -64,7 +69,7 @@ real function get_flow_for_gage(gage) result(flow)
             ierr = nf90_close(ncid)
         end do
 
-        print *, "WARNING: Gage discharge not found in any timeslice file, falling back to fractional diversion"
+        print free, "WARNING: Gage discharge not found in any timeslice file, falling back to fractional diversion"
     end function
 
 end module

src/Routing/module_channel_routing.F

@@ -1961,11 +1961,12 @@ 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))
+
 #ifdef WRF_HYDRO_NUDGING
          call nudge_apply_upstream_muskingumCunge( Qup,  Quc,  nudge(k),  k )
 #endif
 
-         call calculate_diversion(LINKID(k), -1_int64, Quc, tmpQLINK(k,2))
          call SUBMUSKINGCUNGE(&
               tmpQLINK(k,2), velocity(k), qloss(k), LINKID(k),     Qup,        Quc, QLINK(k,1), &
               QLateral(k),   DTRT_CH,     So(k), CHANLEN(k),                  &