diff --git a/CMakeLists.txt b/CMakeLists.txt index e46ebabe2..769737d49 100644 --- a/CMakeLists.txt +++ b/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") diff --git a/src/OrchestratorLayer/config.f90 b/src/OrchestratorLayer/config.f90 index bb0f27f64..bbf39c5cf 100644 --- a/src/OrchestratorLayer/config.f90 +++ b/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 diff --git a/src/Routing/Diversions/module_diversions.F b/src/Routing/Diversions/module_diversions.F index 9e118de75..3fec5dfcf 100644 --- a/src/Routing/Diversions/module_diversions.F +++ b/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 diff --git a/src/Routing/Diversions/module_diversions_timeslice.F b/src/Routing/Diversions/module_diversions_timeslice.F index d5e96e20b..908d8919f 100644 --- a/src/Routing/Diversions/module_diversions_timeslice.F +++ b/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 \ No newline at end of file diff --git a/src/Routing/module_channel_routing.F b/src/Routing/module_channel_routing.F index 8f8d676f4..9e0d4d773 100644 --- a/src/Routing/module_channel_routing.F +++ b/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), &