Skip to content

Commit

Permalink
apply precip to maximum top width of the channel
Browse files Browse the repository at this point in the history
ensure cross sections are defined by at least 2 points
ensure precip is greater than or equal to zero
  • Loading branch information
langevin-usgs committed Jan 23, 2025
1 parent 9dd17c2 commit f8c3a37
Showing 3 changed files with 114 additions and 27 deletions.
2 changes: 1 addition & 1 deletion doc/mf6io/mf6ivar/dfn/chf-cxs.dfn
Original file line number Diff line number Diff line change
@@ -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 ---------------------

93 changes: 77 additions & 16 deletions src/Model/SurfaceWaterFlow/swf-cxs.f90
Original file line number Diff line number Diff line change
@@ -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
46 changes: 36 additions & 10 deletions src/Model/SurfaceWaterFlow/swf-pcp.f90
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit f8c3a37

Please sign in to comment.