From f8c3a37920782564a9dbd467f00a984a14343593 Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D." Date: Thu, 23 Jan 2025 14:36:05 -0600 Subject: [PATCH] apply precip to maximum top width of the channel ensure cross sections are defined by at least 2 points ensure precip is greater than or equal to zero --- doc/mf6io/mf6ivar/dfn/chf-cxs.dfn | 2 +- src/Model/SurfaceWaterFlow/swf-cxs.f90 | 93 +++++++++++++++++++++----- src/Model/SurfaceWaterFlow/swf-pcp.f90 | 46 ++++++++++--- 3 files changed, 114 insertions(+), 27 deletions(-) diff --git a/doc/mf6io/mf6ivar/dfn/chf-cxs.dfn b/doc/mf6io/mf6ivar/dfn/chf-cxs.dfn index 259cf04920c..0fe273aa84c 100644 --- a/doc/mf6io/mf6ivar/dfn/chf-cxs.dfn +++ b/doc/mf6io/mf6ivar/dfn/chf-cxs.dfn @@ -56,7 +56,7 @@ tagged false in_record true reader urword longname number of points used to define cross section -description integer value that defines the number of points used to define the define the shape of a section. NXSPOINTS must be greater than zero or the program will terminate with an error. NXSPOINTS defines the number of points that must be entered for the reach in the CROSSSECTIONDATA block. The sum of NXSPOINTS for all sections must equal the NPOINTS dimension. +description integer value that defines the number of points used to define the define the shape of a section. NXSPOINTS must be greater than 1 or the program will terminate with an error. NXSPOINTS defines the number of points that must be entered for the reach in the CROSSSECTIONDATA block. The sum of NXSPOINTS for all sections must equal the NPOINTS dimension. # --------------------- chf cxs crosssectiondata --------------------- diff --git a/src/Model/SurfaceWaterFlow/swf-cxs.f90 b/src/Model/SurfaceWaterFlow/swf-cxs.f90 index 55981f8e9cc..3e3a6cd1b10 100644 --- a/src/Model/SurfaceWaterFlow/swf-cxs.f90 +++ b/src/Model/SurfaceWaterFlow/swf-cxs.f90 @@ -9,7 +9,7 @@ module SwfCxsModule use MemoryHelperModule, only: create_mem_path use MemoryManagerModule, only: mem_allocate use SimVariablesModule, only: errmsg - use SimModule, only: store_error + use SimModule, only: store_error, store_error_filename, count_errors use NumericalPackageModule, only: NumericalPackageType use BaseDisModule, only: DisBaseType @@ -41,6 +41,7 @@ module SwfCxsModule procedure :: log_dimensions procedure :: source_packagedata procedure :: log_packagedata + procedure :: check_packagedata procedure :: source_crosssectiondata procedure :: log_crosssectiondata procedure :: cxs_da @@ -51,6 +52,7 @@ module SwfCxsModule procedure :: get_conveyance => cxs_conveyance procedure :: get_hydraulic_radius procedure :: get_wetted_top_width + procedure :: get_maximum_top_width procedure :: write_cxs_table end type SwfCxsType @@ -288,44 +290,47 @@ end subroutine allocate_arrays !> @brief Copy options from IDM into package !< subroutine source_packagedata(this) - ! -- modules + ! modules use KindModule, only: LGP use MemoryManagerExtModule, only: mem_set_value use SimVariablesModule, only: idm_context use SwfCxsInputModule, only: SwfCxsParamFoundType - ! -- dummy + ! dummy class(SwfCxsType) :: this - ! -- locals + ! locals character(len=LENMEMPATH) :: idmMemoryPath type(SwfCxsParamFoundType) :: found - ! - ! -- set memory path + + ! set memory path idmMemoryPath = create_mem_path(this%name_model, 'CXS', idm_context) - ! - ! -- update defaults with idm sourced values + + ! update defaults with idm sourced values call mem_set_value(this%idcxs, 'IDCXS', idmMemoryPath, & found%idcxs) call mem_set_value(this%nxspoints, 'NXSPOINTS', idmMemoryPath, & found%nxspoints) - ! - ! -- ensure idcxs was found + + ! ensure idcxs was found if (.not. found%idcxs) then write (errmsg, '(a)') 'Error in PACKAGEDATA block: IDCXS not found.' call store_error(errmsg) end if - ! - ! -- ensure nxspoints was found + + ! ensure nxspoints was found if (.not. found%nxspoints) then write (errmsg, '(a)') 'Error in PACKAGEDATA block: NXSPOINTS not found.' call store_error(errmsg) end if - ! - ! -- log values to list file + + ! log values to list file if (this%iout > 0) then call this%log_packagedata(found) end if - ! - ! -- Calculate the iacross index array using nxspoints + + ! Check to make sure package data is valid + call this%check_packagedata() + + ! Calculate the iacross index array using nxspoints call calc_iacross(this%nxspoints, this%iacross) end subroutine source_packagedata @@ -341,6 +346,41 @@ subroutine calc_iacross(nxspoints, iacross) end do end subroutine calc_iacross + !> @brief Check packagedata + !< + subroutine check_packagedata(this) + ! dummy arguments + class(SwfCxsType) :: this !< this instance + ! local variables + integer(I4B) :: i + + ! Check that all cross section IDs are in range + do i = 1, size(this%idcxs) + if (this%idcxs(i) <= 0 .or. this%idcxs(i) > this%nsections) then + write (errmsg, '(a, i0, a)') & + 'IDCXS values must be greater than 0 and less than NSECTIONS. & + &Found ', this%idcxs(i), '.' + call store_error(errmsg) + end if + end do + + ! Check that nxspoints are greater than one + do i = 1, size(this%nxspoints) + if (this%nxspoints(i) <= 1) then + write (errmsg, '(a, i0, a, i0, a)') & + 'NXSPOINTS values must be greater than 1 for each cross section. & + &Found ', this%nxspoints(i), ' for cross section ', this%idcxs(i), '.' + call store_error(errmsg) + end if + end do + + ! write summary of package error messages + if (count_errors() > 0) then + call store_error_filename(this%input_fname) + end if + + end subroutine check_packagedata + !> @brief Write user packagedata to list file !< subroutine log_packagedata(this, found) @@ -746,4 +786,25 @@ function get_wetted_top_width(this, idcxs, width, depth) result(r) end if end function get_wetted_top_width + function get_maximum_top_width(this, idcxs, width) result(r) + ! modules + use SwfCxsUtilsModule, only: get_saturated_topwidth + ! dummy + class(SwfCxsType) :: this + integer(I4B), intent(in) :: idcxs !< cross section id + real(DP), intent(in) :: width !< width in reach + ! local + real(DP) :: r !< calculated hydraulic radius + integer(I4B) :: i0 + integer(I4B) :: i1 + integer(I4B) :: npts + integer(I4B) :: icalcmeth + call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth) + if (npts == 0) then + r = width + else + r = get_saturated_topwidth(npts, this%xfraction(i0:i1)) + end if + end function get_maximum_top_width + end module SwfCxsModule diff --git a/src/Model/SurfaceWaterFlow/swf-pcp.f90 b/src/Model/SurfaceWaterFlow/swf-pcp.f90 index 80dc8402e85..28c85938b96 100644 --- a/src/Model/SurfaceWaterFlow/swf-pcp.f90 +++ b/src/Model/SurfaceWaterFlow/swf-pcp.f90 @@ -11,7 +11,7 @@ module SwfPcpModule use MemoryHelperModule, only: create_mem_path use BndModule, only: BndType use BndExtModule, only: BndExtType - use SimModule, only: store_error, store_error_filename + use SimModule, only: store_error, store_error_filename, count_errors use SimVariablesModule, only: errmsg use ObsModule, only: DefaultObsIdProcessor use TimeArraySeriesLinkModule, only: TimeArraySeriesLinkType @@ -50,6 +50,7 @@ module SwfPcpModule procedure :: log_pcp_options procedure :: read_initial_attr => pcp_read_initial_attr procedure :: bnd_rp => pcp_rp + procedure :: bnd_ck => pcp_ck procedure :: bnd_cf => pcp_cf procedure :: bnd_fc => pcp_fc procedure :: bnd_da => pcp_da @@ -266,6 +267,35 @@ subroutine pcp_rp(this) end if end subroutine pcp_rp + !> @brief Ensure precipitation is positive + !< + subroutine pcp_ck(this) + ! dummy + class(SwfPcpType), intent(inout) :: this + ! local + character(len=30) :: nodestr + integer(I4B) :: i, nr + character(len=*), parameter :: fmterr = & + &"('Specified stress ',i0, & + &' precipitation (',g0,') is less than zero for cell', a)" + + ! Ensure precipitation rates are positive + do i = 1, this%nbound + nr = this%nodelist(i) + if (nr <= 0) cycle + if (this%precipitation(i) < DZERO) then + call this%dis%noder_to_string(nr, nodestr) + write (errmsg, fmt=fmterr) i, this%precipitation(i), trim(nodestr) + call store_error(errmsg) + end if + end do + + ! write summary of package error messages + if (count_errors() > 0) then + call store_error_filename(this%input_fname) + end if + end subroutine pcp_ck + !> @brief Formulate the HCOF and RHS terms !! !! Skip if no precipitation. Otherwise, calculate hcof and rhs @@ -280,9 +310,7 @@ subroutine pcp_cf(this) real(DP) :: qpcp real(DP) :: area real(DP) :: width_channel - real(DP) :: wetted_top_width - real(DP) :: depth - real(DP) :: dummy + real(DP) :: top_width real(DP), dimension(:), pointer :: reach_length ! Return if no precipitation @@ -309,15 +337,13 @@ subroutine pcp_cf(this) ! Determine the water surface area if (this%dis%is_2d()) then + ! this is for overland flow case area = this%dis%get_area(node) else if (this%dis%is_1d()) then - ! this is for channel case; could be improved with newton formulation + ! this is for channel case idcxs = this%dfw%idcxs(node) - call this%dis%get_flow_width(node, node, 0, width_channel, dummy) - depth = this%xnew(node) - this%dis%bot(node) - wetted_top_width = this%cxs%get_wetted_top_width(idcxs, width_channel, & - depth) - area = reach_length(node) * wetted_top_width + top_width = this%cxs%get_maximum_top_width(idcxs, width_channel) + area = reach_length(node) * top_width end if ! calculate volumetric precipitation flow in L^3/T