diff --git a/CMakeLists.txt b/CMakeLists.txt index 7e353b010..03698f70c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,13 +6,13 @@ endif() # set project name and version numbers -project(WRF_Hydro LANGUAGES Fortran) -set(WRF_Hydro_VERSION_MAJOR 5) -set(WRF_Hydro_VERSION_MINOR 3) -set(WRF_Hydro_VERSION_PATCH 0) -set(National_Water_Model_VERSION_MAJOR 3) -set(National_Water_Model_VERSION_MINOR 0) -set(National_Water_Model_VERSION_PATCH beta) +project (WRF_Hydro LANGUAGES Fortran C) +set (WRF_Hydro_VERSION_MAJOR 5) +set (WRF_Hydro_VERSION_MINOR 4) +set (WRF_Hydro_VERSION_PATCH 0) +set (National_Water_Model_VERSION_MAJOR 3) +set (National_Water_Model_VERSION_MINOR 1) +set (National_Water_Model_VERSION_PATCH beta) # set cmake to work with MPI Fortran find_package(MPI REQUIRED) @@ -211,7 +211,7 @@ 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/CMakeLists.txt b/src/CMakeLists.txt index ed12b2dd2..753db9580 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,7 @@ add_subdirectory("Debug_Utilities") add_subdirectory("Routing/Overland") add_subdirectory("Routing/Subsurface") add_subdirectory("Routing/Reservoirs") +add_subdirectory("Routing/Diversions") add_subdirectory("Data_Rec") add_subdirectory("Routing") add_subdirectory("HYDRO_drv") @@ -112,13 +113,15 @@ if (HYDRO_LSM MATCHES "NoahMP") if (WRF_HYDRO_NUDGING STREQUAL "1") target_link_libraries(wrfhydro.exe hydro_nudging) + target_link_libraries(wrfhydro.exe hydro_routing_diversions) add_dependencies(wrfhydro.exe hydro_nudging) + add_dependencies(wrfhydro.exe hydro_routing_diversions) endif() # bash commands to copy namelists to the Run directory set(BASH_CP_HRLDAS_NML "if [[ ! -f ${CMAKE_BINARY_DIR}/Run/namelist.hrldas ]]\; then cp ${PROJECT_SOURCE_DIR}/src/template/NoahMP/namelist.hrldas ${CMAKE_BINARY_DIR}/Run \; fi\;") set(BASH_CP_HYDRO_NML "if [[ ! -f ${CMAKE_BINARY_DIR}/Run/hydro.namelist ]]\; then cp ${PROJECT_SOURCE_DIR}/src/template/HYDRO/hydro.namelist ${CMAKE_BINARY_DIR}/Run \; fi\;") - + add_custom_command(TARGET wrfhydro.exe POST_BUILD COMMAND mkdir -p ${CMAKE_BINARY_DIR}/Run COMMAND cp ${PROJECT_SOURCE_DIR}/tests/ctests/run_dir_makefile.mk ${CMAKE_BINARY_DIR}/Run/Makefile diff --git a/src/Data_Rec/module_namelist.F90 b/src/Data_Rec/module_namelist.F90 index c72a4275d..083aee4ec 100644 --- a/src/Data_Rec/module_namelist.F90 +++ b/src/Data_Rec/module_namelist.F90 @@ -33,6 +33,7 @@ subroutine read_rt_nlst(nlst) character(len=256) :: route_lake_f="" character(len=256) :: route_direction_f="" character(len=256) :: route_order_f="" + character(len=256) :: diversions_file="" logical :: reservoir_persistence_usgs logical :: reservoir_persistence_usace character(len=256) :: reservoir_parameter_file="" @@ -105,8 +106,8 @@ subroutine read_rt_nlst(nlst) RT_OPTION, CHANRTSWCRT, channel_option, & SUBRTSWCRT,OVRTSWCRT,AGGFACTRT, dtrt_ter,dtrt_ch,dxrt,& GwSpinCycles, GwPreCycles, GwSpinUp, GwPreDiag, GwPreDiagInterval, gwIhShift, & - GWBASESWCRT, gwChanCondSw, gwChanCondConstIn, gwChanCondConstOut , & - route_topo_f,route_chan_f,route_link_f, compound_channel, route_lake_f, & + GWBASESWCRT, gwChanCondSw, gwChanCondConstIn, gwChanCondConstOut, & + route_topo_f,route_chan_f,route_link_f, compound_channel, route_lake_f, diversions_file, & reservoir_persistence_usgs, reservoir_persistence_usace, reservoir_parameter_file, reservoir_usgs_timeslice_path, & reservoir_usace_timeslice_path, reservoir_observation_lookback_hours, reservoir_observation_update_time_interval_seconds, & reservoir_rfc_forecasts, reservoir_rfc_forecasts_time_series_path, reservoir_rfc_forecasts_lookback_hours, & @@ -248,6 +249,7 @@ subroutine read_rt_nlst(nlst) nlst%DEEPGWSPIN = DEEPGWSPIN nlst%SOLVEG_INITSWC = SOLVEG_INITSWC nlst%reservoir_obs_dir = "testDirectory" + nlst%diversions_file = diversions_file nlst%reservoir_persistence_usgs = reservoir_persistence_usgs nlst%reservoir_persistence_usace = reservoir_persistence_usace nlst%reservoir_parameter_file = reservoir_parameter_file diff --git a/src/Data_Rec/module_namelist_inc.F90 b/src/Data_Rec/module_namelist_inc.F90 index c34b9d768..52922903c 100644 --- a/src/Data_Rec/module_namelist_inc.F90 +++ b/src/Data_Rec/module_namelist_inc.F90 @@ -42,6 +42,7 @@ module module_namelist_inc character(len=256) :: route_chan_f="" character(len=256) :: route_link_f="" character(len=256) :: route_lake_f="" + character(len=256) :: diversions_file="" logical :: reservoir_persistence_usgs logical :: reservoir_persistence_usace character(len=256) :: reservoir_parameter_file="" diff --git a/src/HYDRO_drv/module_HYDRO_drv.F90 b/src/HYDRO_drv/module_HYDRO_drv.F90 index a5a6a37df..0efbc46ef 100644 --- a/src/HYDRO_drv/module_HYDRO_drv.F90 +++ b/src/HYDRO_drv/module_HYDRO_drv.F90 @@ -29,6 +29,7 @@ module module_HYDRO_drv #endif use module_hydro_stop, only: HYDRO_stop use module_UDMAP, only: get_basn_area_nhd + use module_channel_diversions, only: init_diversions use netcdf implicit none @@ -1583,6 +1584,11 @@ subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp) if(nlst(did)%CHANRTSWCRT .ne. 0) call init_stream_nudging #endif +!#ifdef WRF_HYDRO_DIVERSIONS +! TODO: should this check to make sure we have nudging on too? [RC] + call init_diversions(nlst(did)%diversions_file, nlst(did)%timeSlicePath) +!#endif + ! if (trim(nlst_rt(did)%restart_file) == "") then ! output at the initial time diff --git a/src/OrchestratorLayer/config.F90 b/src/OrchestratorLayer/config.F90 index e600e1d28..20bc49bcc 100644 --- a/src/OrchestratorLayer/config.F90 +++ b/src/OrchestratorLayer/config.F90 @@ -124,6 +124,7 @@ module config_base character(len=256) :: route_link_f="" character(len=256) :: route_lake_f="" integer :: lake_option + character(len=256) :: diversions_file="" logical :: reservoir_persistence_usgs logical :: reservoir_persistence_usace character(len=256) :: reservoir_parameter_file="" @@ -170,7 +171,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 @@ -554,6 +555,7 @@ subroutine init_namelist_rt_field(did) integer :: channel_loss_option = 0 character(len=256) :: route_lake_f="" integer :: lake_option !0: lakes off 1: level pool 2: passthrough, 3: reservoir da + character(len=256) :: diversions_file="" logical :: reservoir_persistence_usgs logical :: reservoir_persistence_usace character(len=256) :: reservoir_parameter_file="" @@ -628,7 +630,7 @@ subroutine init_namelist_rt_field(did) SUBRTSWCRT,OVRTSWCRT,AGGFACTRT, dtrt_ter,dtrt_ch,dxrt,& GwSpinCycles, GwPreCycles, GwSpinUp, GwPreDiag, GwPreDiagInterval, gwIhShift, & GWBASESWCRT, gwChanCondSw, gwChanCondConstIn, gwChanCondConstOut , & - route_topo_f,route_chan_f,route_link_f, compound_channel, channel_loss_option, lake_option, route_lake_f, & + route_topo_f,route_chan_f,route_link_f, compound_channel, channel_loss_option, lake_option, route_lake_f, diversions_file, & route_direction_f,route_order_f,gwbasmskfil, & geo_finegrid_flnm, gwstrmfil,GW_RESTART,RSTRT_SWC,TERADJ_SOLAR, sys_cpl, & order_to_write , rst_typ, rst_bi_in, rst_bi_out, gwsoilcpl, & @@ -876,6 +878,8 @@ subroutine init_namelist_rt_field(did) nlst(did)%route_link_f = route_link_f nlst(did)%route_lake_f = route_lake_f + nlst(did)%diversions_file = diversions_file + nlst(did)%reservoir_persistence_usgs = reservoir_persistence_usgs nlst(did)%reservoir_persistence_usace = reservoir_persistence_usace nlst(did)%reservoir_parameter_file = reservoir_parameter_file diff --git a/src/Routing/CMakeLists.txt b/src/Routing/CMakeLists.txt index 0fbf3d557..ae80c4dad 100644 --- a/src/Routing/CMakeLists.txt +++ b/src/Routing/CMakeLists.txt @@ -29,4 +29,5 @@ target_link_libraries(hydro_routing hydro_routing_reservoirs_hybrid hydro_data_rec hydro_routing_reservoirs_rfc + hydro_routing_diversions ) diff --git a/src/Routing/Diversions/CMakeLists.txt b/src/Routing/Diversions/CMakeLists.txt new file mode 100644 index 000000000..c2b688564 --- /dev/null +++ b/src/Routing/Diversions/CMakeLists.txt @@ -0,0 +1,10 @@ +add_library(hydro_routing_diversions STATIC + module_diversions.F90 + module_diversions_timeslice.F90 +) + +add_dependencies(hydro_routing_diversions hydro_orchestrator) +add_dependencies(hydro_routing_diversions fortglob) + +target_link_libraries(hydro_routing_diversions PUBLIC hydro_orchestrator) +target_link_libraries(hydro_routing_diversions PUBLIC fortglob) diff --git a/src/Routing/Diversions/module_diversions.F90 b/src/Routing/Diversions/module_diversions.F90 new file mode 100644 index 000000000..11e2bad9c --- /dev/null +++ b/src/Routing/Diversions/module_diversions.F90 @@ -0,0 +1,177 @@ +module module_channel_diversions + use netcdf + use iso_fortran_env, only: int8, int16, int64 + use ieee_arithmetic, only: ieee_is_nan + + use module_diversions_timeslice, only: get_flow_for_gage, init_timeslices + use module_hydro_stop, only: hydro_stop + + implicit none + + type diversion_t + character(len=128) :: name + + character(len=16) :: da_src, da_dest + integer(kind=int8) :: type_div, type_src, type_dest + integer(kind=int64) :: id_src, id_dest, src_index, dest_index + real :: capacity, fraction + integer(kind=int16) :: lookback + + real :: persisted_flow_src, persisted_flow_dest + end type + + logical :: diversions_active = .false. + integer :: ndivs = 0 + + type(diversion_t), allocatable :: diversions(:) + + character(*), parameter :: free = '(*(g0,1x))' + +contains + subroutine init_diversions(diversions_file, timeslice_path) + character(*), intent(in) :: diversions_file + character(*), intent(in) :: timeslice_path + + integer :: g, i, ierr = 0 + character(len=20) :: istr + character(len=256) :: char_tmp + + integer :: ncid, dimid + integer :: name_vid, type_div_vid, type_src_vid, type_dest_vid, da_src_vid, da_dest_vid + integer :: id_src_vid, id_dest_vid, capacity_vid, fraction_vid, lookback_vid + + if (len_trim(diversions_file) > 0) then + print *, "Loading diversions data from " // trim(diversions_file) + ierr = nf90_open(trim(diversions_file), NF90_NOWRITE, ncid) + if (ierr /= 0) call hydro_stop("Could not open diversions file: " // trim(diversions_file)) + ierr = nf90_inq_dimid(ncid, "diversion", dimid) + if (ierr /= 0) call hydro_stop("Error reading diversions file: " // trim(diversions_file)) + ierr = nf90_inquire_dimension(ncid, dimid, len=ndivs) + if (ierr /= 0) call hydro_stop("Error reading diversions file: " // trim(diversions_file)) + + write (istr, *) ndivs + print *, "Diversions file has " // trim(adjustl(istr)) // " diversions" + + ! get fields + ierr = 0 + ierr = ierr + nf90_inq_varid(ncid, "Diversion_Name", name_vid) + ierr = ierr + nf90_inq_varid(ncid, "DivType", type_div_vid) + ierr = ierr + nf90_inq_varid(ncid, "FromType", type_src_vid) + ierr = ierr + nf90_inq_varid(ncid, "ToType", type_dest_vid) + ierr = ierr + nf90_inq_varid(ncid, "DA_Src", da_src_vid) + ierr = ierr + nf90_inq_varid(ncid, "DA_Dest", da_dest_vid) + ierr = ierr + nf90_inq_varid(ncid, "DivFrom", id_src_vid) + ierr = ierr + nf90_inq_varid(ncid, "DivTo", id_dest_vid) + ierr = ierr + nf90_inq_varid(ncid, "DivCap", capacity_vid) + ierr = ierr + nf90_inq_varid(ncid, "DivFrac", fraction_vid) + ierr = ierr + nf90_inq_varid(ncid, "Lookback", lookback_vid) + + if (ierr /= 0) call hydro_stop("Error occurred accessing diversion file variables") + + if (ndivs > 0) then + diversions_active = .true. + + ! Read the timeslice data + ierr = init_timeslices(timeslice_path) + if (ierr /= 0) call hydro_stop("No timeslice files available when initializing diversions") + + allocate(diversions(ndivs)) + do i = 1, ndivs + associate (div => diversions(i)) + div = diversion_t('', '', '', -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1) + + ierr = 0 + !ierr = ierr + nf_get_var1(ncid, name_vid, i, div%name) !! can't read string with Fortran :-( + + ierr = ierr + nf90_get_var(ncid, da_src_vid, div%da_src, start=(/i/), count=(/15/)) + div%persisted_flow_src = get_flow_for_gage(div%da_src) + ierr = ierr + nf90_get_var(ncid, da_dest_vid, div%da_dest, start=(/i/), count=(/15/)) + div%persisted_flow_dest = get_flow_for_gage(div%da_dest) + + ierr = ierr + nf90_get_var(ncid, type_div_vid, div%type_div, start=(/i/)) + ierr = ierr + nf90_get_var(ncid, type_src_vid, div%type_src, start=(/i/)) + ierr = ierr + nf90_get_var(ncid, type_dest_vid, div%type_dest, start=(/i/)) + ierr = ierr + nf90_get_var(ncid, id_src_vid, div%id_src, start=(/i/)) + ierr = ierr + nf90_get_var(ncid, id_dest_vid, div%id_dest, start=(/i/)) + ierr = ierr + nf90_get_var(ncid, capacity_vid, div%capacity, start=(/i/)) + ierr = ierr + nf90_get_var(ncid, fraction_vid, div%fraction, start=(/i/)) + ierr = ierr + nf90_get_var(ncid, lookback_vid, div%lookback, start=(/i/)) + + if (ierr /= 0) call hydro_stop("Error occurred reading diversion variables from diversion file") + end associate + end do + + end if + end if + + end subroutine + + 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_out + real, intent(in) :: qlink_src_in + real, intent(out) :: diversion_quantity_out, diversion_quantity_in + + integer :: i + + diversion_quantity_out = 0 + diversion_quantity_in = 0 + + ! bail if we're inactive + if (.not. diversions_active) return + + ! link to gage + ! look to see what type of diversion it is + ! call into sub-procedure to handle type=1, type=2, type=3, etc + + do i = 1, ndivs + 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_in, diversions(i), qlink_src_in, diversion_quantity_out) + ! 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) + integer(kind=int64), intent(in) :: src_link + type(diversion_t), intent(in) :: diversion + real, intent(in) :: qlink_src + real, intent(out) :: div_gage_flow + + 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). + ! + ! If it's not a valid flow, we try to use the Fraction property of the diversion, + ! and if -that's- not available, we just leave the flow untouched. + + div_gage_flow = diversion%persisted_flow_dest + if (ieee_is_nan(div_gage_flow)) then + fraction = diversion%fraction + if (fraction == -1) then + print free, "WARNING: No fractional diversion value specified for diversion at gage '" // trim(adjustl(diversion%da_dest)) // "', skipping" + fraction = 0 + 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 + end if + end subroutine + +end module diff --git a/src/Routing/Diversions/module_diversions_timeslice.F90 b/src/Routing/Diversions/module_diversions_timeslice.F90 new file mode 100644 index 000000000..908d8919f --- /dev/null +++ b/src/Routing/Diversions/module_diversions_timeslice.F90 @@ -0,0 +1,75 @@ +module module_diversions_timeslice + use fortglob, only: glob_t, globfiles + use module_hydro_stop, only: hydro_stop + + use netcdf + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan + implicit none + + integer, parameter :: PATH_MAX = 4096 + type(glob_t) :: timeslice_files + + character(*), parameter :: free = '(*(g0,1x))' + + contains + + integer function init_timeslices(timeslice_path) result(ierr) + character(*), intent(in) :: timeslice_path + character(len=PATH_MAX) :: tslice_glob + + ierr = 0 + tslice_glob = trim(adjustl(timeslice_path)) // '/' // "*.15min.usgsTimeSlice.ncdf" + timeslice_files = globfiles(tslice_glob) + + if (timeslice_files%nfiles == 0) ierr = 1 + end function + + real function get_flow_for_gage(gage) result(flow) + character(len=15), intent(in) :: gage + + integer :: i, ierr=0, ncid, dimid, varid, num_stns, found(1) + real :: discharge(1) + character(len=15), allocatable :: gage_ids(:) + + flow = ieee_value(flow, ieee_quiet_nan) + + if (gage(1:4) == 'None') then + return + end if + + ! start looking at files, going backward from most recent + do i = timeslice_files%nfiles, 1, -1 + ierr = nf90_open(trim(timeslice_files%filenames(i)), NF90_NOWRITE, ncid) + + ! look for gage + ierr = ierr + nf90_inq_dimid(ncid, 'stationIdInd', dimid) + ierr = ierr + nf90_inquire_dimension(ncid, dimid, len=num_stns) + + allocate(gage_ids(num_stns)) + ierr = ierr + nf90_inq_varid(ncid, 'stationId', varid) + ierr = ierr + nf90_get_var(ncid, varid, gage_ids) + + if (ierr /= 0) call hydro_stop("Error occurred reading gage data from " // trim(timeslice_files%filenames(i))) + + 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))) + + flow = discharge(1) + deallocate(gage_ids) + return + end if + + deallocate(gage_ids) + ierr = nf90_close(ncid) + end do + + 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.F90 b/src/Routing/module_channel_routing.F90 index c9d64962a..4be8c7d01 100644 --- a/src/Routing/module_channel_routing.F90 +++ b/src/Routing/module_channel_routing.F90 @@ -509,10 +509,10 @@ subroutine SUBMUSKINGCUNGE( & !comment out to prevent excessive file sizes when running model !print*, "qloss,dx,WP,WPk,depth,ChannK,qdc,ql,dt,D", qloss,dx,WP,WPk,depth,ChannK,qdc,ql,dt,D if((qloss*dt)/D > ((ql*dt)/D - C4)) then - qloss = ql - C4*(D/dt) - if (qloss < 0) then - print*, 'WARNING CHANNEL LOSS IS NEGATIVE',qloss - endif + qloss = ql - C4*(D/dt) + if ((qloss < 0) .and. (ChannK /= 0)) then + print*, 'WARNING CHANNEL LOSS IS NEGATIVE',qloss + endif endif ! ---------------------------------------------------------------- @@ -1645,7 +1645,7 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & nudgeWAdvance, & nudge_apply_upstream_muskingumCunge #endif - + use module_channel_diversions, only: calculate_diversion implicit none @@ -1756,6 +1756,10 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & real, allocatable,dimension(:) :: tmpAssimilatedValue character(len=256), allocatable,dimension(:) :: tmpAssimilatedSourceFile + ! diversions + real :: div_src, div_dst + character(*), parameter :: free = '(*(g0,1x))' + #ifdef MPP_LAND if(my_id .eq. io_id) then #endif @@ -2005,6 +2009,32 @@ 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) ) + ! HANDLE DIVERSIONS + + 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_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 + + if (div_dst /= 0) then + ! apply observed value to downstream +#ifdef HYDRO_D + print free, "DEBUG: diverting", div_dst, "to link id =", LINKID(k) !, "on processor", my_id +#endif + Qup = div_dst + Quc = div_dst + tmpQLINK(k,1) = div_dst + tmpQLINK(k,2) = div_dst + end if + #ifdef WRF_HYDRO_NUDGING call nudge_apply_upstream_muskingumCunge( Qup, Quc, nudge(k), k ) #endif diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index 923017f10..5d838d291 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -18,3 +18,5 @@ add_library(hydro_utils STATIC module_version.F90 module_hydro_stop.F90 ) + +add_subdirectory(fortglob) \ No newline at end of file diff --git a/src/utils/fortglob/CMakeLists.txt b/src/utils/fortglob/CMakeLists.txt new file mode 100644 index 000000000..bf41788c3 --- /dev/null +++ b/src/utils/fortglob/CMakeLists.txt @@ -0,0 +1,4 @@ +add_library(fortglob STATIC + libfortglob.c + fortglob.f90 +) \ No newline at end of file diff --git a/src/utils/fortglob/fortglob.f90 b/src/utils/fortglob/fortglob.f90 new file mode 100644 index 000000000..658c537d2 --- /dev/null +++ b/src/utils/fortglob/fortglob.f90 @@ -0,0 +1,76 @@ +module fortglob + use iso_c_binding + implicit none + + private + public :: glob_t, globfiles + + integer, parameter :: PATH_MAX = 4096 + + type :: glob_t + integer :: nfiles + character(len=PATH_MAX), allocatable :: filenames(:) + end type + + type, bind(c) :: globrec_c + integer (c_size_t) :: nfiles + type(c_ptr) :: filenames + type(c_ptr) :: globptr + end type globrec_c + + interface + function globfiles_c(pattern) bind(c) + use iso_c_binding + character(kind=c_char), dimension(*), intent(in) :: pattern + type(c_ptr) :: globfiles_c + end function globfiles_c + + subroutine freeglobrec(recptr) bind (c) + use iso_c_binding + type(c_ptr), intent(in) :: recptr + end subroutine + + function strlen(s) bind(c) + use iso_c_binding + type(c_ptr), intent(in), value :: s + integer(c_size_t) :: strlen + end function strlen + end interface + + contains + + function globfiles(pattern) + use iso_c_binding + + type(glob_t) :: globfiles + character(len=*), intent(in) :: pattern + + integer :: i, j, slen + character(len=PATH_MAX) :: pattern_c + type(c_ptr) :: glob_c_ptr + type(globrec_c), pointer :: glob_f_ptr + type(c_ptr), pointer :: c_strs(:) + character(kind=c_char), pointer, dimension(:) :: temp_fstrp + + pattern_c = pattern + slen = len_trim(pattern_c) + pattern_c(slen+1:slen+1) = c_null_char + glob_c_ptr = globfiles_c(pattern_c) + call c_f_pointer(glob_c_ptr, glob_f_ptr) + globfiles%nfiles = glob_f_ptr%nfiles + + allocate(globfiles%filenames(globfiles%nfiles)) + + call c_f_pointer(glob_f_ptr%filenames, c_strs, [glob_f_ptr%nfiles]) + do i = 1, globfiles%nfiles + slen = strlen(c_strs(i)) + call c_f_pointer(c_strs(i), temp_fstrp, [slen]) + globfiles%filenames(i) = '' + do j = 1, slen + globfiles%filenames(i)(j:j) = temp_fstrp(j) + end do + end do + + call freeglobrec(glob_f_ptr%globptr) + end function +end module diff --git a/src/utils/fortglob/libfortglob.c b/src/utils/fortglob/libfortglob.c new file mode 100644 index 000000000..5a188746b --- /dev/null +++ b/src/utils/fortglob/libfortglob.c @@ -0,0 +1,19 @@ +#include "libfortglob.h" + +globrec * globfiles_c(const char * pattern) { + globrec *record = malloc(sizeof(globrec)); + glob_t *globbuf = malloc(sizeof(glob_t)); + + glob(pattern, 0, NULL, globbuf); + record->nfiles = globbuf->gl_pathc; + record->filenames = globbuf->gl_pathv; + record->__globptr = globbuf; + + return record; +} + +void freeglobrec(glob_t ** rec) { + glob_t *glob = *rec; + globfree(glob); + free(glob); +} diff --git a/src/utils/fortglob/libfortglob.h b/src/utils/fortglob/libfortglob.h new file mode 100644 index 000000000..aec96808f --- /dev/null +++ b/src/utils/fortglob/libfortglob.h @@ -0,0 +1,13 @@ +#include +#include +#include +#include + +typedef struct globrec { + size_t nfiles; + char ** filenames; + glob_t * __globptr; +} globrec; + +globrec * globfiles_c(const char * pattern); +void freeglobrec(glob_t **rec); diff --git a/tests/utils/xrcmp.py b/tests/utils/xrcmp.py index 6cc6d781b..4c22bb1e0 100644 --- a/tests/utils/xrcmp.py +++ b/tests/utils/xrcmp.py @@ -11,6 +11,7 @@ import math from multiprocessing import Pool +import os import pathlib import sys # import time @@ -139,6 +140,7 @@ def xrcmp( n_cores: int = 1, chunks={}, exclude_vars: list = [], + interactive: bool = False ) -> int: if exclude_vars is None: @@ -147,7 +149,7 @@ def xrcmp( # Delete log file first # Should write a log file that says nothing yet determined? log_file = pathlib.Path(log_file) - if log_file.exists(): + if log_file.exists() and not interactive: log_file.unlink() # Dont chunk, this is just a meta-data read. @@ -183,8 +185,8 @@ def xrcmp( diff_var_names = sorted(all_stats.keys()) if not diff_var_names: - with open(log_file, 'w') as opened_file: - opened_file.write("Files are identical\n") + with open(log_file, 'a') as opened_file: + opened_file.write(f"Files are identical: {os.path.basename(ref_file)}\n") return 0 # Formatting: @@ -259,6 +261,11 @@ def xrcmp( header_dict = {name: name for name in stat_names} the_header = header_string.format(**header_dict) + if len(all_stats) > 0: + print(the_header, end='') + for key in all_stats.keys(): + print(var_string.format(**all_stats[key]), end='') + with open(log_file, 'w') as opened_file: opened_file.write(the_header) for key in all_stats.keys(): @@ -294,29 +301,35 @@ def parse_arguments(): default=1, help="Chunks as integer." ) + parser.add_argument( + "-i", "--interactive", action='store_true', default=False, + help="Run in interactive mode for local analysis." + ) + args = parser.parse_args() can_file = args.candidate ref_file = args.reference log_file = args.log_file chunks = args.chunks n_cores = args.n_cores + interactive = args.interactive if chunks == 1: chunks = {} # No chunking elif chunks == 0: chunks = None # This will use the conus_chunks_dict - return can_file, ref_file, log_file, chunks, n_cores - + return can_file, ref_file, log_file, chunks, n_cores, interactive if __name__ == "__main__": - can_file, ref_file, log_file, chunks, n_cores = parse_arguments() + can_file, ref_file, log_file, chunks, n_cores, inter = parse_arguments() ret = xrcmp( can_file=can_file, ref_file=ref_file, log_file=log_file, n_cores=n_cores, - chunks=chunks + chunks=chunks, + interactive=inter ) sys.exit(ret) diff --git a/trunk/NDHMS b/trunk/NDHMS deleted file mode 120000 index e057607ed..000000000 --- a/trunk/NDHMS +++ /dev/null @@ -1 +0,0 @@ -../src/ \ No newline at end of file