sfrmodule Module Reference

This module contains the SFR package methods. More...

Data Types
type	sfrtype

Functions/Subroutines
subroutine, public	sfr_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname)
	@ brief Create a new package object More...
subroutine	sfr_allocate_scalars (this)
	@ brief Allocate scalars More...
subroutine	sfr_allocate_arrays (this)
	@ brief Allocate arrays More...
subroutine	sfr_read_dimensions (this)
	@ brief Read dimensions for package More...
subroutine	sfr_options (this, option, found)
	@ brief Read additional options for package More...
subroutine	sfr_ar (this)
	@ brief Allocate and read method for package More...
subroutine	sfr_read_packagedata (this)
	@ brief Read packagedata for the package More...
subroutine	sfr_read_crossection (this)
	@ brief Read crosssection block for the package More...
subroutine	sfr_read_connectiondata (this)
	@ brief Read connectiondata for the package More...
subroutine	sfr_read_diversions (this)
	@ brief Read diversions for the package More...
subroutine	sfr_rp (this)
	@ brief Read and prepare period data for package More...
subroutine	sfr_ad (this)
	@ brief Advance the package More...
subroutine	sfr_cf (this)
	@ brief Formulate the package hcof and rhs terms. More...
subroutine	sfr_fc (this, rhs, ia, idxglo, matrix_sln)
	@ brief Copy hcof and rhs terms into solution. More...
subroutine	sfr_fn (this, rhs, ia, idxglo, matrix_sln)
	@ brief Add Newton-Raphson terms for package into solution. More...
subroutine	sfr_cc (this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
	@ brief Convergence check for package. More...
subroutine	sfr_cq (this, x, flowja, iadv)
	@ brief Calculate package flows. More...
subroutine	sfr_ot_package_flows (this, icbcfl, ibudfl)
	@ brief Output package flow terms. More...
subroutine	sfr_ot_dv (this, idvsave, idvprint)
	@ brief Output package dependent-variable terms. More...
subroutine	sfr_ot_bdsummary (this, kstp, kper, iout, ibudfl)
	@ brief Output advanced package budget summary. More...
subroutine	sfr_da (this)
	@ brief Deallocate package memory More...
subroutine	define_listlabel (this)
	@ brief Define the list label for the package More...
logical function	sfr_obs_supported (this)
	Determine if observations are supported. More...
subroutine	sfr_df_obs (this)
	Define the observation types available in the package. More...
subroutine	sfr_bd_obs (this)
	Save observations for the package. More...
subroutine	sfr_rp_obs (this)
	Read and prepare observations for a package. More...
subroutine	sfr_process_obsid (obsrv, dis, inunitobs, iout)
	Process observation IDs for a package. More...
subroutine	sfr_set_stressperiod (this, n, ichkustrm, crossfile)
	Set period data. More...
subroutine	sfr_solve (this, n, h, hcof, rhs, update)
	Solve reach continuity equation. More...
subroutine	sfr_update_flows (this, n, qd, qgwf)
	Update flow terms. More...
subroutine	sfr_adjust_ro_ev (this, qc, qu, qi, qr, qro, qe, qfrommvr)
	Adjust runoff and evaporation. More...
subroutine	sfr_calc_qd (this, n, depth, hgwf, qgwf, qd)
	Calculate downstream flow term. More...
subroutine	sfr_calc_qsource (this, n, depth, qsrc)
	Calculate sum of sources. More...
subroutine	sfr_calc_qman (this, n, depth, qman)
	Calculate streamflow. More...
subroutine	sfr_calc_qgwf (this, n, depth, hgwf, qgwf, gwfhcof, gwfrhs)
	Calculate reach-aquifer exchange. More...
integer(i4b) function	sfr_gwf_conn (this, n)
	Determine if a reach is connected to a gwf cell. More...
subroutine	sfr_calc_cond (this, n, depth, cond, hsfr, h_temp)
	Calculate reach-aquifer conductance. More...
subroutine	sfr_calc_constant (this, n, d1, hgwf, qgwf, qd)
subroutine	sfr_calc_steady (this, n, d1, hgwf, qu, qi, qfrommvr, qr, qe, qro, qgwf, qd)
subroutine	sfr_calc_div (this, n, i, qd, qdiv)
	Calculate diversion flow. More...
subroutine	sfr_calc_reach_depth (this, n, q1, d1)
	Calculate the depth at the midpoint. More...
subroutine	sfr_calc_xs_depth (this, n, qrch, d)
	Calculate the depth at the midpoint of a irregular cross-section. More...
subroutine	sfr_check_conversion (this)
	Check unit conversion data. More...
subroutine	sfr_check_reaches (this)
	Check reach data. More...
subroutine	sfr_check_connections (this)
	Check connection data. More...
subroutine	sfr_check_diversions (this)
	Check diversions data. More...
subroutine	sfr_check_ustrf (this)
	Check upstream fraction data. More...
subroutine	sfr_setup_budobj (this)
	Setup budget object for package. More...
subroutine	sfr_fill_budobj (this)
	Copy flow terms into budget object for package. More...
subroutine	sfr_setup_tableobj (this)
	Setup stage table object for package. More...
real(dp) function	calc_area_wet (this, n, depth)
	Calculate wetted area. More...
real(dp) function	calc_perimeter_wet (this, n, depth)
	Calculate wetted perimeter. More...
real(dp) function	calc_surface_area (this, n)
	Calculate maximum surface area. More...
real(dp) function	calc_surface_area_wet (this, n, depth)
	Calculate wetted surface area. More...
real(dp) function	calc_top_width_wet (this, n, depth)
	Calculate wetted top width. More...
subroutine	sfr_activate_density (this)
	Activate density terms. More...
subroutine	sfr_activate_viscosity (this)
	Activate viscosity terms. More...
subroutine	sfr_calculate_density_exchange (this, n, stage, head, cond, bots, flow, gwfhcof, gwfrhs)
	Calculate density terms. More...

Variables
character(len=lenftype)	ftype = 'SFR'
	package ftype string More...
character(len=lenpackagename)	text = ' SFR'
	package budget string More...

Detailed Description

This module contains the overridden methods for the streamflow routing (SFR) package. Most of the methods in the base Boundary Package are overridden.

Function/Subroutine Documentation

calc_area_wet()

real(dp) function sfrmodule::calc_area_wet ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth  )

private

Function to calculate the wetted area for a SFR package reach.

Returns

wetted area

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth

Definition at line 5747 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_area_wet !< wetted area
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local variables
    integer(I4B) :: npts
    integer(I4B) :: i0
    integer(I4B) :: i1
    !
    ! -- Calculate wetted area
    npts = this%ncrosspts(n)
    i0 = this%iacross(n)
    i1 = this%iacross(n + 1) - 1
    if (npts > 1) then
      calc_area_wet = get_cross_section_area(npts, this%station(i0:i1), &
                                             this%xsheight(i0:i1), depth)
    else
      calc_area_wet = this%station(i0) * depth
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

calc_perimeter_wet()

real(dp) function sfrmodule::calc_perimeter_wet ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth  )

private

Function to calculate the wetted perimeter for a SFR package reach.

Returns

wetted perimeter

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth

Definition at line 5779 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_perimeter_wet !< wetted perimeter
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local variables
    integer(I4B) :: npts
    integer(I4B) :: i0
    integer(I4B) :: i1
    !
    ! -- Calculate wetted perimeter
    npts = this%ncrosspts(n)
    i0 = this%iacross(n)
    i1 = this%iacross(n + 1) - 1
    if (npts > 1) then
      calc_perimeter_wet = get_wetted_perimeter(npts, this%station(i0:i1), &
                                                this%xsheight(i0:i1), depth)
    else
      calc_perimeter_wet = this%station(i0) ! no depth dependence in original implementation
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

calc_surface_area()

real(dp) function sfrmodule::calc_surface_area ( class(sfrtype)  this, integer(i4b), intent(in)  n  )

private

Function to calculate the maximum surface area for a SFR package reach.

Returns

surface area

Parameters

	this	SfrType object
[in]	n	reach number

Definition at line 5811 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_surface_area !< surface area
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    ! -- local variables
    integer(I4B) :: npts
    integer(I4B) :: i0
    integer(I4B) :: i1
    real(DP) :: top_width
    !
    ! -- Calculate surface area
    npts = this%ncrosspts(n)
    i0 = this%iacross(n)
    i1 = this%iacross(n + 1) - 1
    if (npts > 1) then
      top_width = get_saturated_topwidth(npts, this%station(i0:i1))
    else
      top_width = this%station(i0)
    end if
    calc_surface_area = top_width * this%length(n)
    !
    ! -- return
    return

Here is the call graph for this function:

calc_surface_area_wet()

real(dp) function sfrmodule::calc_surface_area_wet ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth  )

private

Function to calculate the wetted surface area for a SFR package reach.

Returns

wetted surface area

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth

Definition at line 5843 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_surface_area_wet !< wetted surface area
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local variables
    real(DP) :: top_width
    !
    ! -- Calculate wetted surface area
    top_width = this%calc_top_width_wet(n, depth)
    calc_surface_area_wet = top_width * this%length(n)
    !
    ! -- return
    return

calc_top_width_wet()

real(dp) function sfrmodule::calc_top_width_wet ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth  )

private

Function to calculate the wetted top width for a SFR package reach.

Returns

wetted top width

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth

Definition at line 5866 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_top_width_wet !< wetted top width
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local variables
    integer(I4B) :: npts
    integer(I4B) :: i0
    integer(I4B) :: i1
    real(DP) :: sat
    !
    ! -- Calculate wetted top width
    npts = this%ncrosspts(n)
    i0 = this%iacross(n)
    i1 = this%iacross(n + 1) - 1
    sat = scubicsaturation(dem5, dzero, depth, dem5)
    if (npts > 1) then
      calc_top_width_wet = sat * get_wetted_topwidth(npts, &
                                                     this%station(i0:i1), &
                                                     this%xsheight(i0:i1), &
                                                     depth)
    else
      calc_top_width_wet = sat * this%station(i0)
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

define_listlabel()

subroutine sfrmodule::define_listlabel

(

class(sfrtype), intent(inout)

this

)

Method defined the list label for the SFR package. The list label is the heading that is written to iout when PRINT_INPUT option is used.

Parameters

[in,out]

this

SfrType object

Definition at line 2774 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    !
    ! -- create the header list label
    this%listlabel = trim(this%filtyp)//' NO.'
    if (this%dis%ndim == 3) then
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
    elseif (this%dis%ndim == 2) then
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
    else
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
    end if
    write (this%listlabel, '(a, a16)') trim(this%listlabel), 'STRESS RATE'
    if (this%inamedbound == 1) then
      write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
    end if
    !
    ! -- return
    return

sfr_activate_density()

subroutine sfrmodule::sfr_activate_density ( class(sfrtype), intent(inout)  this )

private

Method to activate addition of density terms for a SFR package reach.

Parameters

[in,out]

this

SfrType object

Definition at line 5902 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    !
    ! -- Set idense and reallocate denseterms to be of size MAXBOUND
    this%idense = 1
    call mem_reallocate(this%denseterms, 3, this%MAXBOUND, 'DENSETERMS', &
                        this%memoryPath)
    do i = 1, this%maxbound
      do j = 1, 3
        this%denseterms(j, i) = dzero
      end do
    end do
    write (this%iout, '(/1x,a)') 'DENSITY TERMS HAVE BEEN ACTIVATED FOR SFR &
      &PACKAGE: '//trim(adjustl(this%packName))
    !
    ! -- return
    return

sfr_activate_viscosity()

subroutine sfrmodule::sfr_activate_viscosity

(

class(sfrtype), intent(inout)

this

)

Method to activate addition of viscosity terms for exchange with groundwater along a SFR package reach.

Parameters

[in,out]

this

SfrType object

Definition at line 5933 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    !
    ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND
    this%ivsc = 1
    call mem_reallocate(this%viscratios, 2, this%MAXBOUND, 'VISCRATIOS', &
                        this%memoryPath)
    do i = 1, this%maxbound
      do j = 1, 2
        this%viscratios(j, i) = done
      end do
    end do
    write (this%iout, '(/1x,a)') 'VISCOSITY HAS BEEN ACTIVATED FOR SFR &
      &PACKAGE: '//trim(adjustl(this%packName))
    !
    ! -- return
    return

sfr_ad()

subroutine sfrmodule::sfr_ad

(

class(sfrtype)

this

)

Advance data in the SFR package. The method sets advances time series, time array series, and observation data.

Parameters

this	SfrType object

Definition at line 1909 of file gwf-sfr.f90.

    ! -- modules
    use timeseriesmanagermodule, only: var_timeseries
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: iaux
    !
    ! -- Most advanced package AD routines have to restore state if
    !    the solution failed and the time step is being retried with a smaller
    !    step size.  This is not needed here because there is no old stage
    !    or storage effects in the stream.
    !
    ! -- Advance the time series manager
    call this%TsManager%ad()
    !
    ! -- check upstream fractions if time series are being used to
    !    define this variable
    if (var_timeseries(this%tsManager, this%packName, 'USTRF')) then
      call this%sfr_check_ustrf()
    end if
    !
    ! -- update auxiliary variables by copying from the derived-type time
    !    series variable into the bndpackage auxvar variable so that this
    !    information is properly written to the GWF budget file
    if (this%naux > 0) then
      do n = 1, this%maxbound
        do iaux = 1, this%naux
          if (this%noupdateauxvar(iaux) /= 0) cycle
          this%auxvar(iaux, n) = this%rauxvar(iaux, n)
        end do
      end do
    end if
    !
    ! -- reset upstream flow to zero and set specified stage
    do n = 1, this%maxbound
      this%usflow(n) = dzero
      if (this%iboundpak(n) < 0) then
        this%stage(n) = this%sstage(n)
      end if
    end do
    !
    ! -- pakmvrobj ad
    if (this%imover == 1) then
      call this%pakmvrobj%ad()
    end if
    !
    ! -- For each observation, push simulated value and corresponding
    !    simulation time from "current" to "preceding" and reset
    !    "current" value.
    call this%obs%obs_ad()
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_adjust_ro_ev()

subroutine sfrmodule::sfr_adjust_ro_ev ( class(sfrtype)  this, real(dp), intent(inout)  qc, real(dp), intent(in)  qu, real(dp), intent(in)  qi, real(dp), intent(in)  qr, real(dp), intent(inout)  qro, real(dp), intent(inout)  qe, real(dp), intent(in)  qfrommvr  )

private

Method to adjust runoff and evaporation for a SFR package reach based on the total reach flow.

Parameters

	this	SfrType object
[in,out]	qc	total reach volumetric flow
[in]	qu	upstream reach volumetric flow
[in]	qi	reach volumetric inflow
[in]	qr	reach volumetric rainfall
[in,out]	qro	reach volumetric runoff
[in,out]	qe	reach volumetric evaporation
[in]	qfrommvr	reach volumetric flow from mover

Definition at line 3625 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    real(DP), intent(inout) :: qc !< total reach volumetric flow
    real(DP), intent(in) :: qu !< upstream reach volumetric flow
    real(DP), intent(in) :: qi !< reach volumetric inflow
    real(DP), intent(in) :: qr !< reach volumetric rainfall
    real(DP), intent(inout) :: qro !< reach volumetric runoff
    real(DP), intent(inout) :: qe !< reach volumetric evaporation
    real(DP), intent(in) :: qfrommvr !< reach volumetric flow from mover
    ! -- local variables
    real(DP) :: qt
    !
    ! -- adjust runoff or evaporation if sum of sources is negative
    if (qc < dzero) then
      !
      ! -- calculate sources without evaporation
      qt = qu + qi + qr + qro + qfrommvr
      !
      ! -- runoff exceeds sources of water for reach
      if (qt < dzero) then
        if (qro < dzero) then
          qro = -(qu + qi + qr + qfrommvr)
          qe = dzero
        end if
        !
        ! -- evaporation exceeds sources of water for reach
      else
        if (qe > dzero) then
          qe = qu + qi + qr + qro + qfrommvr
        end if
      end if
      qc = qu + qi + qr - qe + qro + qfrommvr
    end if

sfr_allocate_arrays()

subroutine sfrmodule::sfr_allocate_arrays

(

class(sfrtype), intent(inout)

this

)

Allocate and initialize array for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 351 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    !
    ! -- allocate character array for budget text
    allocate (this%csfrbudget(this%bditems))
    call mem_allocate(this%sfrname, lenboundname, this%maxbound, &
                      'SFRNAME', this%memoryPath)
    !
    ! -- variables originally in SfrDataType
    call mem_allocate(this%iboundpak, this%maxbound, 'IBOUNDPAK', &
                      this%memoryPath)
    call mem_allocate(this%igwfnode, this%maxbound, 'IGWFNODE', this%memoryPath)
    call mem_allocate(this%igwftopnode, this%maxbound, 'IGWFTOPNODE', &
                      this%memoryPath)
    call mem_allocate(this%length, this%maxbound, 'LENGTH', this%memoryPath)
    call mem_allocate(this%width, this%maxbound, 'WIDTH', this%memoryPath)
    call mem_allocate(this%strtop, this%maxbound, 'STRTOP', this%memoryPath)
    call mem_allocate(this%bthick, this%maxbound, 'BTHICK', this%memoryPath)
    call mem_allocate(this%hk, this%maxbound, 'HK', this%memoryPath)
    call mem_allocate(this%slope, this%maxbound, 'SLOPE', this%memoryPath)
    call mem_allocate(this%nconnreach, this%maxbound, 'NCONNREACH', &
                      this%memoryPath)
    call mem_allocate(this%ustrf, this%maxbound, 'USTRF', this%memoryPath)
    call mem_allocate(this%ftotnd, this%maxbound, 'FTOTND', this%memoryPath)
    call mem_allocate(this%ndiv, this%maxbound, 'NDIV', this%memoryPath)
    call mem_allocate(this%usflow, this%maxbound, 'USFLOW', this%memoryPath)
    call mem_allocate(this%dsflow, this%maxbound, 'DSFLOW', this%memoryPath)
    call mem_allocate(this%depth, this%maxbound, 'DEPTH', this%memoryPath)
    call mem_allocate(this%stage, this%maxbound, 'STAGE', this%memoryPath)
    call mem_allocate(this%gwflow, this%maxbound, 'GWFLOW', this%memoryPath)
    call mem_allocate(this%simevap, this%maxbound, 'SIMEVAP', this%memoryPath)
    call mem_allocate(this%simrunoff, this%maxbound, 'SIMRUNOFF', &
                      this%memoryPath)
    call mem_allocate(this%stage0, this%maxbound, 'STAGE0', this%memoryPath)
    call mem_allocate(this%usflow0, this%maxbound, 'USFLOW0', this%memoryPath)
    !
    ! -- reach order and connection data
    call mem_allocate(this%isfrorder, this%maxbound, 'ISFRORDER', &
                      this%memoryPath)
    call mem_allocate(this%ia, this%maxbound + 1, 'IA', this%memoryPath)
    call mem_allocate(this%ja, 0, 'JA', this%memoryPath)
    call mem_allocate(this%idir, 0, 'IDIR', this%memoryPath)
    call mem_allocate(this%idiv, 0, 'IDIV', this%memoryPath)
    call mem_allocate(this%qconn, 0, 'QCONN', this%memoryPath)
    !
    ! -- boundary data
    call mem_allocate(this%rough, this%maxbound, 'ROUGH', this%memoryPath)
    call mem_allocate(this%rain, this%maxbound, 'RAIN', this%memoryPath)
    call mem_allocate(this%evap, this%maxbound, 'EVAP', this%memoryPath)
    call mem_allocate(this%inflow, this%maxbound, 'INFLOW', this%memoryPath)
    call mem_allocate(this%runoff, this%maxbound, 'RUNOFF', this%memoryPath)
    call mem_allocate(this%sstage, this%maxbound, 'SSTAGE', this%memoryPath)
    !
    ! -- aux variables
    call mem_allocate(this%rauxvar, this%naux, this%maxbound, &
                      'RAUXVAR', this%memoryPath)
    !
    ! -- diversion variables
    call mem_allocate(this%iadiv, this%maxbound + 1, 'IADIV', this%memoryPath)
    call mem_allocate(this%divreach, 0, 'DIVREACH', this%memoryPath)
    call mem_allocate(this%divflow, 0, 'DIVFLOW', this%memoryPath)
    call mem_allocate(this%divq, 0, 'DIVQ', this%memoryPath)
    !
    ! -- cross-section data
    call mem_allocate(this%ncrosspts, this%maxbound, 'NCROSSPTS', &
                      this%memoryPath)
    call mem_allocate(this%iacross, this%maxbound + 1, 'IACROSS', &
                      this%memoryPath)
    call mem_allocate(this%station, this%ncrossptstot, 'STATION', &
                      this%memoryPath)
    call mem_allocate(this%xsheight, this%ncrossptstot, 'XSHEIGHT', &
                      this%memoryPath)
    call mem_allocate(this%xsrough, this%ncrossptstot, 'XSROUGH', &
                      this%memoryPath)
    !
    ! -- initialize variables
    this%iacross(1) = 0
    do i = 1, this%maxbound
      this%iboundpak(i) = 1
      this%igwfnode(i) = 0
      this%igwftopnode(i) = 0
      this%length(i) = dzero
      this%width(i) = dzero
      this%strtop(i) = dzero
      this%bthick(i) = dzero
      this%hk(i) = dzero
      this%slope(i) = dzero
      this%nconnreach(i) = 0
      this%ustrf(i) = dzero
      this%ftotnd(i) = dzero
      this%ndiv(i) = 0
      this%usflow(i) = dzero
      this%dsflow(i) = dzero
      this%depth(i) = dzero
      this%stage(i) = dzero
      this%gwflow(i) = dzero
      this%simevap(i) = dzero
      this%simrunoff(i) = dzero
      this%stage0(i) = dzero
      this%usflow0(i) = dzero
      !
      ! -- boundary data
      this%rough(i) = dzero
      this%rain(i) = dzero
      this%evap(i) = dzero
      this%inflow(i) = dzero
      this%runoff(i) = dzero
      this%sstage(i) = dzero
      !
      ! -- aux variables
      do j = 1, this%naux
        this%rauxvar(j, i) = dzero
      end do
      !
      ! -- cross-section data
      this%ncrosspts(i) = 0
      this%iacross(i + 1) = 0
    end do
    !
    ! -- initialize additional cross-section data
    do i = 1, this%ncrossptstot
      this%station(i) = dzero
      this%xsheight(i) = dzero
      this%xsrough(i) = dzero
    end do
    !
    !-- fill csfrbudget
    this%csfrbudget(1) = '        RAINFALL'
    this%csfrbudget(2) = '     EVAPORATION'
    this%csfrbudget(3) = '          RUNOFF'
    this%csfrbudget(4) = '      EXT-INFLOW'
    this%csfrbudget(5) = '             GWF'
    this%csfrbudget(6) = '     EXT-OUTFLOW'
    this%csfrbudget(7) = '        FROM-MVR'
    this%csfrbudget(8) = '          TO-MVR'
    !
    ! -- allocate and initialize budget output data
    call mem_allocate(this%qoutflow, this%maxbound, 'QOUTFLOW', this%memoryPath)
    call mem_allocate(this%qextoutflow, this%maxbound, 'QEXTOUTFLOW', &
                      this%memoryPath)
    do i = 1, this%maxbound
      this%qoutflow(i) = dzero
      this%qextoutflow(i) = dzero
    end do
    !
    ! -- allocate and initialize dbuff
    if (this%istageout > 0) then
      call mem_allocate(this%dbuff, this%maxbound, 'DBUFF', this%memoryPath)
      do i = 1, this%maxbound
        this%dbuff(i) = dzero
      end do
    else
      call mem_allocate(this%dbuff, 0, 'DBUFF', this%memoryPath)
    end if
    !
    ! -- allocate character array for budget text
    allocate (this%cauxcbc(this%cbcauxitems))
    !
    ! -- allocate and initialize qauxcbc
    call mem_allocate(this%qauxcbc, this%cbcauxitems, 'QAUXCBC', &
                      this%memoryPath)
    do i = 1, this%cbcauxitems
      this%qauxcbc(i) = dzero
    end do
    !
    ! -- fill cauxcbc
    this%cauxcbc(1) = 'FLOW-AREA       '
    !
    ! -- allocate denseterms to size 0
    call mem_allocate(this%denseterms, 3, 0, 'DENSETERMS', this%memoryPath)
    !
    ! -- allocate viscratios to size 0
    call mem_allocate(this%viscratios, 2, 0, 'VISCRATIOS', this%memoryPath)
    !
    ! -- return
    return

sfr_allocate_scalars()

subroutine sfrmodule::sfr_allocate_scalars

(

class(sfrtype), intent(inout)

this

)

Allocate and initialize scalars for the SFR package. The base model allocate scalars method is also called.

Parameters

[in,out]

this

SfrType object

Definition at line 282 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    use memoryhelpermodule, only: create_mem_path
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    !
    ! -- call standard BndType allocate scalars
    call this%BndType%allocate_scalars()
    !
    ! -- allocate the object and assign values to object variables
    call mem_allocate(this%iprhed, 'IPRHED', this%memoryPath)
    call mem_allocate(this%istageout, 'ISTAGEOUT', this%memoryPath)
    call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
    call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
    call mem_allocate(this%ipakcsv, 'IPAKCSV', this%memoryPath)
    call mem_allocate(this%idiversions, 'IDIVERSIONS', this%memoryPath)
    call mem_allocate(this%maxsfrpicard, 'MAXSFRPICARD', this%memoryPath)
    call mem_allocate(this%maxsfrit, 'MAXSFRIT', this%memoryPath)
    call mem_allocate(this%bditems, 'BDITEMS', this%memoryPath)
    call mem_allocate(this%cbcauxitems, 'CBCAUXITEMS', this%memoryPath)
    call mem_allocate(this%unitconv, 'UNITCONV', this%memoryPath)
    call mem_allocate(this%lengthconv, 'LENGTHCONV', this%memoryPath)
    call mem_allocate(this%timeconv, 'TIMECONV', this%memoryPath)
    call mem_allocate(this%dmaxchg, 'DMAXCHG', this%memoryPath)
    call mem_allocate(this%deps, 'DEPS', this%memoryPath)
    call mem_allocate(this%nconn, 'NCONN', this%memoryPath)
    call mem_allocate(this%icheck, 'ICHECK', this%memoryPath)
    call mem_allocate(this%iconvchk, 'ICONVCHK', this%memoryPath)
    call mem_allocate(this%idense, 'IDENSE', this%memoryPath)
    call mem_allocate(this%ianynone, 'IANYNONE', this%memoryPath)
    call mem_allocate(this%ncrossptstot, 'NCROSSPTSTOT', this%memoryPath)
    !
    ! -- set pointer to gwf iss
    call mem_setptr(this%gwfiss, 'ISS', create_mem_path(this%name_model))
    !
    ! -- Set values
    this%iprhed = 0
    this%istageout = 0
    this%ibudgetout = 0
    this%ibudcsv = 0
    this%ipakcsv = 0
    this%idiversions = 0
    this%maxsfrpicard = 100
    this%maxsfrit = maxadpit
    this%bditems = 8
    this%cbcauxitems = 1
    this%unitconv = done
    this%lengthconv = dnodata
    this%timeconv = dnodata
    this%dmaxchg = dem5
    this%deps = dp999 * this%dmaxchg
    this%nconn = 0
    this%icheck = 1
    this%iconvchk = 1
    this%idense = 0
    this%ivsc = 0
    this%ianynone = 0
    this%ncrossptstot = 0
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_ar()

subroutine sfrmodule::sfr_ar

(

class(sfrtype), intent(inout)

this

)

Method to read and prepare period data for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 780 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: ierr
    !
    ! -- allocate and read observations
    call this%obs%obs_ar()
    !
    ! -- call standard BndType allocate scalars
    call this%BndType%allocate_arrays()
    !
    ! -- set boundname for each connection
    if (this%inamedbound /= 0) then
      do n = 1, this%maxbound
        this%boundname(n) = this%sfrname(n)
      end do
    end if
    !
    ! -- copy boundname into boundname_cst
    call this%copy_boundname()
    !
    ! -- copy igwfnode into nodelist
    do n = 1, this%maxbound
      this%nodelist(n) = this%igwfnode(n)
    end do
    !
    ! -- check the sfr unit conversion data
    call this%sfr_check_conversion()
    !
    ! -- check the sfr reach data
    call this%sfr_check_reaches()
 
    ! -- check the connection data
    call this%sfr_check_connections()
 
    ! -- check the diversion data
    if (this%idiversions /= 0) then
      call this%sfr_check_diversions()
    end if
    !
    ! -- terminate if errors were detected in any of the static sfr data
    ierr = count_errors()
    if (ierr > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- setup pakmvrobj
    if (this%imover /= 0) then
      allocate (this%pakmvrobj)
      call this%pakmvrobj%ar(this%maxbound, this%maxbound, this%memoryPath)
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_bd_obs()

subroutine sfrmodule::sfr_bd_obs ( class(sfrtype)  this )

private

Method to save simulated values for the SFR package.

Parameters

this	SfrType object

Definition at line 2926 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: n
    real(DP) :: v
    character(len=100) :: msg
    type(ObserveType), pointer :: obsrv => null()
    !
    ! Write simulated values for all sfr observations
    if (this%obs%npakobs > 0) then
      call this%obs%obs_bd_clear()
      do i = 1, this%obs%npakobs
        obsrv => this%obs%pakobs(i)%obsrv
        do j = 1, obsrv%indxbnds_count
          n = obsrv%indxbnds(j)
          v = dzero
          select case (obsrv%ObsTypeId)
          case ('STAGE')
            v = this%stage(n)
          case ('TO-MVR')
            v = dnodata
            if (this%imover == 1) then
              v = this%pakmvrobj%get_qtomvr(n)
              if (v > dzero) then
                v = -v
              end if
            end if
          case ('FROM-MVR')
            v = dnodata
            if (this%imover == 1) then
              v = this%pakmvrobj%get_qfrommvr(n)
            end if
          case ('EXT-INFLOW')
            v = this%inflow(n)
          case ('INFLOW')
            v = this%usflow(n)
          case ('OUTFLOW')
            v = this%qoutflow(n)
          case ('EXT-OUTFLOW')
            v = this%qextoutflow(n)
          case ('RAINFALL')
            if (this%iboundpak(n) /= 0) then
              v = this%rain(n)
            else
              v = dzero
            end if
          case ('RUNOFF')
            v = this%simrunoff(n)
          case ('EVAPORATION')
            v = this%simevap(n)
          case ('SFR')
            v = this%gwflow(n)
          case ('UPSTREAM-FLOW')
            v = this%usflow(n)
            if (this%imover == 1) then
              v = v + this%pakmvrobj%get_qfrommvr(n)
            end if
          case ('DOWNSTREAM-FLOW')
            v = this%dsflow(n)
            if (v > dzero) then
              v = -v
            end if
          case ('DEPTH')
            v = this%depth(n)
          case ('WET-PERIMETER')
            v = this%calc_perimeter_wet(n, this%depth(n))
          case ('WET-AREA')
            v = this%calc_area_wet(n, this%depth(n))
          case ('WET-WIDTH')
            v = this%calc_top_width_wet(n, this%depth(n))
          case default
            msg = 'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
            call store_error(msg)
          end select
          call this%obs%SaveOneSimval(obsrv, v)
        end do
      end do
      !
      ! -- write summary of package error messages
      if (count_errors() > 0) then
        call this%parser%StoreErrorUnit()
      end if
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_calc_cond()

subroutine sfrmodule::sfr_calc_cond ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth, real(dp), intent(inout)  cond, real(dp), intent(in), optional  hsfr, real(dp), intent(in), optional  h_temp  )

private

Method to calculate the reach-aquifer conductance for a SFR package reach.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth
[in,out]	cond	reach-aquifer conductance
[in]	hsfr	stream stage
[in]	h_temp	head in gw cell

Definition at line 3921 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    real(DP), intent(inout) :: cond !< reach-aquifer conductance
    real(DP), intent(in), optional :: hsfr !< stream stage
    real(DP), intent(in), optional :: h_temp !< head in gw cell
    ! -- local variables
    integer(I4B) :: node
    real(DP) :: wp
    real(DP) :: vscratio
    !
    ! -- initialize conductance
    cond = dzero
    !
    ! -- initial viscosity ratio to 1
    vscratio = done
    !
    ! -- calculate conductance if GWF cell is active
    !    rch-gwf flow will not occur if reach connected to an constant head cell
    node = this%igwfnode(n)
    if (node > 0) then
      if (this%ibound(node) > 0) then
        !
        ! -- direction of gradient across streambed determines which vsc ratio
        if (this%ivsc == 1) then
          if (hsfr > h_temp) then
            ! strm stg > gw head
            vscratio = this%viscratios(1, n)
          else
            vscratio = this%viscratios(2, n)
          end if
        end if
        wp = this%calc_perimeter_wet(n, depth)
        cond = this%hk(n) * vscratio * this%length(n) * wp / this%bthick(n)
      end if
    end if
    !
    ! -- return
    return

sfr_calc_constant()

subroutine sfrmodule::sfr_calc_constant ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(inout)  d1, real(dp), intent(in)  hgwf, real(dp), intent(inout)  qgwf, real(dp), intent(inout)  qd  )

private

Parameters

	this	SfrType object
[in]	n	reach number
[in,out]	d1	current reach depth estimate
[in]	hgwf	head in gw cell
[in,out]	qgwf	reach-aquifer exchange
[in,out]	qd	reach outflow

Definition at line 3964 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(inout) :: d1 !< current reach depth estimate
    real(DP), intent(in) :: hgwf !< head in gw cell
    real(DP), intent(inout) :: qgwf !< reach-aquifer exchange
    real(DP), intent(inout) :: qd !< reach outflow
    ! -- local variables
    real(DP) :: qsrc
 
    this%stage(n) = this%sstage(n)
    d1 = max(dzero, this%stage(n) - this%strtop(n))
 
    call this%sfr_calc_qsource(n, d1, qsrc)
 
    if (this%sfr_gwf_conn(n) == 1) then
      call this%sfr_calc_qgwf(n, d1, hgwf, qgwf)
      qgwf = -qgwf
    else
      qgwf = dzero
    end if
 
    ! -- leakage exceeds inflow
    if (qgwf > qsrc) then
      d1 = dzero
      call this%sfr_calc_qsource(n, d1, qsrc)
      qgwf = qsrc
    end if
 
    ! -- set qd
    qd = qsrc - qgwf
 

sfr_calc_div()

subroutine sfrmodule::sfr_calc_div ( class(sfrtype)  this, integer(i4b), intent(in)  n, integer(i4b), intent(in)  i, real(dp), intent(inout)  qd, real(dp), intent(inout)  qdiv  )

private

Method to calculate the diversion flow for a diversion connected to a SFR package reach. The downstream flow for a reach is passed in and adjusted by the diversion flow amount calculated in this method.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	i	diversion number in reach
[in,out]	qd	remaining downstream flow for reach
[in,out]	qdiv	diversion flow for diversion i

Definition at line 4303 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    integer(I4B), intent(in) :: i !< diversion number in reach
    real(DP), intent(inout) :: qd !< remaining downstream flow for reach
    real(DP), intent(inout) :: qdiv !< diversion flow for diversion i
    ! -- local variables
    character(len=10) :: cp
    integer(I4B) :: jpos
    integer(I4B) :: n2
    real(DP) :: v
    !
    ! -- set local variables
    jpos = this%iadiv(n) + i - 1
    n2 = this%divreach(jpos)
    cp = this%divcprior(jpos)
    v = this%divflow(jpos)
    !
    ! -- calculate diversion
    select case (cp)
      ! -- flood diversion
    case ('EXCESS')
      if (qd < v) then
        v = dzero
      else
        v = qd - v
      end if
      ! -- diversion percentage
    case ('FRACTION')
      v = qd * v
      ! -- STR priority algorithm
    case ('THRESHOLD')
      if (qd < v) then
        v = dzero
      end if
      ! -- specified diversion
    case ('UPTO')
      if (v > qd) then
        v = qd
      end if
    case default
      v = dzero
    end select
    !
    ! -- update upstream from for downstream reaches
    qd = qd - v
    qdiv = v
    !
    ! -- return
    return

sfr_calc_qd()

subroutine sfrmodule::sfr_calc_qd ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth, real(dp), intent(in)  hgwf, real(dp), intent(inout)  qgwf, real(dp), intent(inout)  qd  )

private

Method to calculate downstream flow for a SFR package reach.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth
[in]	hgwf	groundwater head in connected GWF cell
[in,out]	qgwf	groundwater leakage for reach
[in,out]	qd	residual

Definition at line 3666 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    real(DP), intent(in) :: hgwf !< groundwater head in connected GWF cell
    real(DP), intent(inout) :: qgwf !< groundwater leakage for reach
    real(DP), intent(inout) :: qd !< residual
    ! -- local variables
    real(DP) :: qsrc
    !
    ! -- initialize residual
    qd = dzero
    !
    ! -- calculate total water sources excluding groundwater leakage
    call this%sfr_calc_qsource(n, depth, qsrc)
    !
    ! -- estimate groundwater leakage
    call this%sfr_calc_qgwf(n, depth, hgwf, qgwf)
    if (-qgwf > qsrc) qgwf = -qsrc
    !
    ! -- calculate down stream flow
    qd = qsrc + qgwf
    !
    ! -- limit downstream flow to a positive value
    if (qd < dem30) qd = dzero
    !
    ! -- return
    return

sfr_calc_qgwf()

subroutine sfrmodule::sfr_calc_qgwf ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth, real(dp), intent(in)  hgwf, real(dp), intent(inout)  qgwf, real(dp), intent(inout), optional  gwfhcof, real(dp), intent(inout), optional  gwfrhs  )

private

Method to calculate the reach-aquifer exchange for a SFR package reach. The reach-aquifer exchange is relative to the reach. Calculated flow is positive if flow is from the aquifer to the reach.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth
[in]	hgwf	head in GWF cell connected to reach
[in,out]	qgwf	reach-aquifer exchange
[in,out]	gwfhcof	diagonal coefficient term for reach
[in,out]	gwfrhs	right-hand side term for reach

Definition at line 3829 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    real(DP), intent(in) :: hgwf !< head in GWF cell connected to reach
    real(DP), intent(inout) :: qgwf !< reach-aquifer exchange
    real(DP), intent(inout), optional :: gwfhcof !< diagonal coefficient term for reach
    real(DP), intent(inout), optional :: gwfrhs !< right-hand side term for reach
    ! -- local variables
    integer(I4B) :: node
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: hsfr
    real(DP) :: h_temp
    real(DP) :: cond
    real(DP) :: sat
    real(DP) :: derv
    real(DP) :: gwfhcof0
    real(DP) :: gwfrhs0
    !
    ! -- initialize qgwf
    qgwf = dzero
    !
    ! -- skip sfr-aquifer exchange in external cells
    node = this%igwfnode(n)
    if (node < 1) return
    !
    ! -- skip sfr-aquifer exchange in inactive cells
    if (this%ibound(node) == 0) return
    !
    ! -- calculate saturation
    call schsmooth(depth, sat, derv)
    !
    ! -- terms for calculating direction of gradient across streambed
    tp = this%strtop(n)
    bt = tp - this%bthick(n)
    hsfr = tp + depth
    h_temp = hgwf
    if (h_temp < bt) then
      h_temp = bt
    end if
    !
    ! -- calculate conductance
    call this%sfr_calc_cond(n, depth, cond, hsfr, h_temp)
    !
    ! -- calculate groundwater leakage
    qgwf = sat * cond * (h_temp - hsfr)
    gwfrhs0 = -sat * cond * hsfr
    gwfhcof0 = -sat * cond
    !
    ! Add density contributions, if active
    if (this%idense /= 0) then
      call this%sfr_calculate_density_exchange(n, hsfr, hgwf, cond, tp, &
                                               qgwf, gwfhcof0, gwfrhs0)
    end if
    !
    ! -- Set gwfhcof and gwfrhs if present
    if (present(gwfhcof)) gwfhcof = gwfhcof0
    if (present(gwfrhs)) gwfrhs = gwfrhs0
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_calc_qman()

subroutine sfrmodule::sfr_calc_qman ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth, real(dp), intent(inout)  qman  )

private

Method to calculate the streamflow using Manning's equation for a SFR package reach.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth
[in,out]	qman	streamflow

Definition at line 3755 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    real(DP), intent(inout) :: qman !< streamflow
    ! -- local variables
    integer(I4B) :: npts
    integer(I4B) :: i0
    integer(I4B) :: i1
    real(DP) :: sat
    real(DP) :: derv
    real(DP) :: s
    real(DP) :: r
    real(DP) :: aw
    real(DP) :: wp
    real(DP) :: rh
    !
    ! -- initialize variables
    qman = dzero
    !
    ! -- calculate Manning's discharge for non-zero depths
    if (depth > dzero) then
      npts = this%ncrosspts(n)
      !
      ! -- set constant terms for Manning's equation
      call schsmooth(depth, sat, derv)
      s = this%slope(n)
      !
      ! -- calculate the mannings coefficient that is a
      !    function of depth
      if (npts > 1) then
        !
        ! -- get the location of the cross-section data for the reach
        i0 = this%iacross(n)
        i1 = this%iacross(n + 1) - 1
        !
        ! -- get the Manning's sum of the Manning's discharge
        !    for each section
        qman = get_mannings_section(npts, &
                                    this%station(i0:i1), &
                                    this%xsheight(i0:i1), &
                                    this%xsrough(i0:i1), &
                                    this%rough(n), &
                                    this%unitconv, &
                                    s, &
                                    depth)
      else
        r = this%rough(n)
        aw = this%calc_area_wet(n, depth)
        wp = this%calc_perimeter_wet(n, depth)
        if (wp > dzero) then
          rh = aw / wp
        else
          rh = dzero
        end if
        qman = this%unitconv * aw * (rh**dtwothirds) * sqrt(s) / r
      end if
      !
      ! -- calculate stream flow
      qman = sat * qman
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_calc_qsource()

subroutine sfrmodule::sfr_calc_qsource ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  depth, real(dp), intent(inout)  qsrc  )

private

Method to calculate the sum of sources for reach, excluding reach leakage, for a SFR package reach.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	depth	reach depth
[in,out]	qsrc	sum of sources for reach

Definition at line 3703 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    real(DP), intent(inout) :: qsrc !< sum of sources for reach
    ! -- local variables
    real(DP) :: qu
    real(DP) :: qi
    real(DP) :: qr
    real(DP) :: qe
    real(DP) :: qro
    real(DP) :: qfrommvr
    real(DP) :: a
    real(DP) :: ae
    !
    ! -- initialize residual
    qsrc = dzero
    !
    ! -- calculate flow terms
    qu = this%usflow(n)
    qi = this%inflow(n)
    qro = this%runoff(n)
    !
    ! -- calculate rainfall and evap
    a = this%calc_surface_area(n)
    ae = this%calc_surface_area_wet(n, depth)
    qr = this%rain(n) * a
    qe = this%evap(n) * ae
    !
    ! -- calculate mover term
    qfrommvr = dzero
    if (this%imover == 1) then
      qfrommvr = this%pakmvrobj%get_qfrommvr(n)
    end if
    !
    ! -- calculate down stream flow
    qsrc = qu + qi + qr - qe + qro + qfrommvr
    !
    ! -- adjust runoff or evaporation if sum of sources is negative
    call this%sfr_adjust_ro_ev(qsrc, qu, qi, qr, qro, qe, qfrommvr)
    !
    ! -- return
    return

sfr_calc_reach_depth()

subroutine sfrmodule::sfr_calc_reach_depth ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  q1, real(dp), intent(inout)  d1  )

private

Method to calculate the depth at the midpoint of a reach.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	q1	streamflow
[in,out]	d1	stream depth at midpoint of reach

Definition at line 4361 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: q1 !< streamflow
    real(DP), intent(inout) :: d1 !< stream depth at midpoint of reach
    ! -- local variables
    real(DP) :: w
    real(DP) :: s
    real(DP) :: r
    real(DP) :: qconst
    !
    ! -- initialize slope and roughness
    s = this%slope(n)
    r = this%rough(n)
    !
    ! -- calculate stream depth at the midpoint
    if (q1 > dzero) then
      if (this%ncrosspts(n) > 1) then
        call this%sfr_calc_xs_depth(n, q1, d1)
      else
        w = this%station(this%iacross(n))
        qconst = this%unitconv * w * sqrt(s) / r
        d1 = (q1 / qconst)**dp6
      end if
    else
      d1 = dzero
    end if
    !
    ! -- return
    return

sfr_calc_steady()

subroutine sfrmodule::sfr_calc_steady ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(inout)  d1, real(dp), intent(in)  hgwf, real(dp), intent(in)  qu, real(dp), intent(in)  qi, real(dp), intent(in)  qfrommvr, real(dp), intent(in)  qr, real(dp), intent(in)  qe, real(dp), intent(in)  qro, real(dp), intent(inout)  qgwf, real(dp), intent(inout)  qd  )

private

Parameters

	this	SfrType object
[in]	n	reach number
[in,out]	d1	current reach depth estimate
[in]	hgwf	head in gw cell
[in]	qu	reach upstream flow
[in]	qi	reach specified inflow
[in]	qfrommvr	reach flow from mover
[in]	qr	reach rainfall
[in]	qe	reach evaporation
[in]	qro	reach runoff flow
[in,out]	qgwf	reach-aquifer exchange
[in,out]	qd	reach outflow

Definition at line 3999 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(inout) :: d1 !< current reach depth estimate
    real(DP), intent(in) :: hgwf !< head in gw cell
    real(DP), intent(in) :: qu !< reach upstream flow
    real(DP), intent(in) :: qi !< reach specified inflow
    real(DP), intent(in) :: qfrommvr !< reach flow from mover
    real(DP), intent(in) :: qr !< reach rainfall
    real(DP), intent(in) :: qe !< reach evaporation
    real(DP), intent(in) :: qro !< reach runoff flow
    real(DP), intent(inout) :: qgwf !< reach-aquifer exchange
    real(DP), intent(inout) :: qd !< reach outflow
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: isolve
    integer(I4B) :: iic
    integer(I4B) :: iic2
    integer(I4B) :: iic3
    integer(I4B) :: iic4
    integer(I4B) :: ibflg
    real(DP) :: qmp
    real(DP) :: qsrc
    real(DP) :: qsrcmp
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: hsfr
    real(DP) :: qc
    real(DP) :: en1
    real(DP) :: en2
    real(DP) :: qen1
    real(DP) :: f1
    real(DP) :: f2
    real(DP) :: qgwf1
    real(DP) :: qgwf2
    real(DP) :: qgwfp
    real(DP) :: qgwfold
    real(DP) :: fhstr1
    real(DP) :: fhstr2
    real(DP) :: d2
    real(DP) :: dpp
    real(DP) :: dx
    real(DP) :: q1
    real(DP) :: q2
    real(DP) :: derv
    real(DP) :: dlh
    real(DP) :: dlhold
    real(DP) :: fp
    real(DP) :: err
    real(DP) :: errold
    !
    ! -- initialize local variables
    d2 = dzero
    q1 = dzero
    q2 = dzero
    qsrc = dzero
    qgwf = dzero
    qgwfold = dzero
    !
    ! -- calculate the flow at end of the reach
    !    excluding groundwater leakage
    qc = qu + qi + qr - qe + qro + qfrommvr
    !
    ! -- calculate flow at the middle of the reach
    !    excluding groundwater leakage
    qsrcmp = qu + qi + qfrommvr + dhalf * (qr - qe + qro)
    qmp = qsrcmp ! initial estimate flow at the midpoint
    !
    ! -- calculate stream depth at the midpoint
    call this%sfr_calc_reach_depth(n, qmp, d1)
    !
    ! -- calculate sources/sinks for reach
    !    excluding groundwater leakage
    call this%sfr_calc_qsource(n, d1, qsrc)
    !
    ! -- calculate initial reach stage, downstream flow,
    !    and groundwater leakage
    tp = this%strtop(n)
    bt = tp - this%bthick(n)
    hsfr = d1 + tp
    qd = max(qsrc, dzero)
    qgwf = dzero
    !
    ! -- set flag to skip iterations
    isolve = 1
    if (hsfr <= tp .and. hgwf <= tp) isolve = 0
    if (hgwf <= tp .and. qc < dem30) isolve = 0
    if (this%sfr_gwf_conn(n) == 0) isolve = 0
    !
    ! -- iterate to achieve solution
    calc_solution: if (isolve /= 0) then
      !
      ! -- estimate initial end points
      en1 = dzero
      if (d1 > dem30) then
        if ((tp - hgwf) > dem30) then
          en2 = dp9 * d1
        else
          en2 = d1p1 * d1 - (tp - hgwf)
        end if
      else if ((tp - hgwf) > dem30) then
        en2 = done
      else
        en2 = dp99 * (hgwf - tp)
      end if
      !
      ! -- estimate flow at end points
      ! -- end point 1
      if (hgwf > tp) then
        call this%sfr_calc_qgwf(n, dzero, hgwf, qgwf1)
        qgwf1 = -qgwf1
        qen1 = qmp - dhalf * qgwf1
      else
        qgwf1 = dzero
        qen1 = qsrcmp
      end if
      if (hgwf > bt) then
        call this%sfr_calc_qgwf(n, en2, hgwf, qgwf2)
        qgwf2 = -qgwf2
      else
        call this%sfr_calc_qgwf(n, en2, bt, qgwf2)
        qgwf2 = -qgwf2
      end if
      if (qgwf2 > qsrc) qgwf2 = qsrc
      ! -- calculate two depths
      call this%sfr_calc_reach_depth(n, (qsrcmp - dhalf * qgwf1), d1)
      call this%sfr_calc_reach_depth(n, (qsrcmp - dhalf * qgwf2), d2)
      ! -- determine roots
      if (d1 > dem30) then
        f1 = en1 - d1
      else
        en1 = dzero
        f1 = en1 - dzero
      end if
      if (d2 > dem30) then
        f2 = en2 - d2
        if (f2 < dem30) en2 = d2
      else
        d2 = dzero
        f2 = en2 - dzero
      end if
      !
      ! -- iterate to find a solution
      dpp = dhalf * (en1 + en2)
      dx = dpp
      iic = 0
      iic2 = 0
      iic3 = 0
      fhstr1 = dzero
      fhstr2 = dzero
      qgwfp = dzero
      dlhold = dzero
      do i = 1, this%maxsfrit
        ibflg = 0
        d1 = dpp
        d2 = d1 + dtwo * this%deps
        ! -- calculate q at midpoint at both end points
        call this%sfr_calc_qman(n, d1, q1)
        call this%sfr_calc_qman(n, d2, q2)
        ! -- calculate groundwater leakage at both end points
        call this%sfr_calc_qgwf(n, d1, hgwf, qgwf1)
        qgwf1 = -qgwf1
        call this%sfr_calc_qgwf(n, d2, hgwf, qgwf2)
        qgwf2 = -qgwf2
        !
        if (qgwf1 >= qsrc) then
          en2 = dpp
          dpp = dhalf * (en1 + en2)
          call this%sfr_calc_qgwf(n, dpp, hgwf, qgwfp)
          qgwfp = -qgwfp
          if (qgwfp > qsrc) qgwfp = qsrc
          call this%sfr_calc_reach_depth(n, (qsrcmp - dhalf * qgwfp), dx)
          ibflg = 1
        else
          fhstr1 = (qsrcmp - dhalf * qgwf1) - q1
          fhstr2 = (qsrcmp - dhalf * qgwf2) - q2
        end if
        !
        if (ibflg == 0) then
          derv = dzero
          if (abs(d1 - d2) > dzero) then
            derv = (fhstr1 - fhstr2) / (d1 - d2)
          end if
          if (abs(derv) > dem30) then
            dlh = -fhstr1 / derv
          else
            dlh = dzero
          end if
          dpp = d1 + dlh
          !
          ! -- updated depth outside of endpoints - use bisection instead
          if ((dpp >= en2) .or. (dpp <= en1)) then
            if (abs(dlh) > abs(dlhold) .or. dpp < dem30) then
              ibflg = 1
              dpp = dhalf * (en1 + en2)
            end if
          end if
          !
          ! -- check for slow convergence
          ! -- set flags to determine if the Newton-Raphson method oscillates
          !    or if convergence is slow
          if (qgwf1 * qgwfold < dem30) then
            iic2 = iic2 + 1
          else
            iic2 = 0
          end if
          if (qgwf1 < dem30) then
            iic3 = iic3 + 1
          else
            iic3 = 0
          end if
          if (dlh * dlhold < dem30 .or. abs(dlh) > abs(dlhold)) then
            iic = iic + 1
          end if
          iic4 = 0
          if (iic3 > 7 .and. iic > 12) then
            iic4 = 1
          end if
          !
          ! -- switch to bisection when the Newton-Raphson method oscillates
          !    or when convergence is slow
          if (iic2 > 7 .or. iic > 12 .or. iic4 == 1) then
            ibflg = 1
            dpp = dhalf * (en1 + en2)
          end if
          !
          ! -- Calculate perturbed gwf flow
          call this%sfr_calc_qgwf(n, dpp, hgwf, qgwfp)
          qgwfp = -qgwfp
          if (qgwfp > qsrc) then
            qgwfp = qsrc
            if (abs(en1 - en2) < this%dmaxchg * dem6) then
              call this%sfr_calc_reach_depth(n, (qsrcmp - dhalf * qgwfp), dpp)
            end if
          end if
          call this%sfr_calc_reach_depth(n, (qsrcmp - dhalf * qgwfp), dx)
        end if
        !
        ! -- bisection to update end points
        fp = dpp - dx
        if (ibflg == 1) then
          dlh = fp
          ! -- change end points
          ! -- root is between f1 and fp
          if (f1 * fp < dzero) then
            en2 = dpp
            f2 = fp
            ! -- root is between fp and f2
          else
            en1 = dpp
            f1 = fp
          end if
          err = min(abs(fp), abs(en2 - en1))
        else
          err = abs(dlh)
        end if
        !
        ! -- check for convergence and exit if converged
        if (err < this%dmaxchg) then
          d1 = dpp
          qgwf = qgwfp
          qd = qsrc - qgwf
          exit
        end if
        !
        ! -- save iterates
        errold = err
        dlhold = dlh
        if (ibflg == 1) then
          qgwfold = qgwfp
        else
          qgwfold = qgwf1
        end if
        !
        ! -- end of iteration
      end do
      ! -- depth = 0 and hgwf < bt
    else
      call this%sfr_calc_qgwf(n, d1, hgwf, qgwf)
      qgwf = -qgwf
      !
      ! -- leakage exceeds inflow
      if (qgwf > qsrc) then
        d1 = dzero
        call this%sfr_calc_qsource(n, d1, qsrc)
        qgwf = qsrc
      end if
      ! -- set qd
      qd = qsrc - qgwf
    end if calc_solution
 

sfr_calc_xs_depth()

subroutine sfrmodule::sfr_calc_xs_depth ( class(sfrtype)  this, integer(i4b), intent(in)  n, real(dp), intent(in)  qrch, real(dp), intent(inout)  d  )

private

Method to calculate the depth at the midpoint of a reach with a irregular cross-section using Newton-Raphson.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	qrch	streamflow
[in,out]	d	stream depth at midpoint of reach

Definition at line 4400 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: qrch !< streamflow
    real(DP), intent(inout) :: d !< stream depth at midpoint of reach
    ! -- local variables
    integer(I4B) :: iter
    real(DP) :: perturbation
    real(DP) :: q0
    real(DP) :: q1
    real(DP) :: dq
    real(DP) :: derv
    real(DP) :: dd
    real(DP) :: residual
    !
    ! -- initialize variables
    perturbation = this%deps * dtwo
    d = dzero
    q0 = dzero
    residual = q0 - qrch
    !
    ! -- Newton-Raphson iteration
    nriter: do iter = 1, this%maxsfrit
      call this%sfr_calc_qman(n, d + perturbation, q1)
      dq = (q1 - q0)
      if (dq /= dzero) then
        derv = perturbation / (q1 - q0)
      else
        derv = dzero
      end if
      dd = derv * residual
      d = d - dd
      call this%sfr_calc_qman(n, d, q0)
      residual = q0 - qrch
      !
      ! -- check for convergence
      if (abs(dd) < this%dmaxchg) then
        exit nriter
      end if
    end do nriter
    !
    ! -- return
    return

sfr_calculate_density_exchange()

subroutine sfrmodule::sfr_calculate_density_exchange	(	class(sfrtype), intent(inout)	this,
		integer(i4b), intent(in)	n,
		real(dp), intent(in)	stage,
		real(dp), intent(in)	head,
		real(dp), intent(in)	cond,
		real(dp), intent(in)	bots,
		real(dp), intent(inout)	flow,
		real(dp), intent(inout)	gwfhcof,
		real(dp), intent(inout)	gwfrhs
	)

Method to galculate groundwater-reach density exchange terms for a SFR package reach.

Member variable used here denseterms : shape (3, MAXBOUND), filled by buoyancy package col 1 is relative density of sfr (densesfr / denseref) col 2 is relative density of gwf cell (densegwf / denseref) col 3 is elevation of gwf cell

Parameters

[in,out]	this	SfrType object
[in]	n	reach number
[in]	stage	reach stage
[in]	head	head in connected GWF cell
[in]	cond	reach conductance
[in]	bots	bottom elevation of reach
[in,out]	flow	calculated flow, updated here with density terms
[in,out]	gwfhcof	GWF diagonal coefficient, updated here with density terms
[in,out]	gwfrhs	GWF right-hand-side value, updated here with density terms

Definition at line 5970 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: stage !< reach stage
    real(DP), intent(in) :: head !< head in connected GWF cell
    real(DP), intent(in) :: cond !< reach conductance
    real(DP), intent(in) :: bots !< bottom elevation of reach
    real(DP), intent(inout) :: flow !< calculated flow, updated here with density terms
    real(DP), intent(inout) :: gwfhcof !< GWF diagonal coefficient, updated here with density terms
    real(DP), intent(inout) :: gwfrhs !< GWF right-hand-side value, updated here with density terms
    ! -- local variables
    real(DP) :: ss
    real(DP) :: hh
    real(DP) :: havg
    real(DP) :: rdensesfr
    real(DP) :: rdensegwf
    real(DP) :: rdenseavg
    real(DP) :: elevsfr
    real(DP) :: elevgwf
    real(DP) :: elevavg
    real(DP) :: d1
    real(DP) :: d2
    logical(LGP) :: stage_below_bot
    logical(LGP) :: head_below_bot
    !
    ! -- Set sfr density to sfr density or gwf density
    if (stage >= bots) then
      ss = stage
      stage_below_bot = .false.
      rdensesfr = this%denseterms(1, n) ! sfr rel density
    else
      ss = bots
      stage_below_bot = .true.
      rdensesfr = this%denseterms(2, n) ! gwf rel density
    end if
    !
    ! -- set hh to head or bots
    if (head >= bots) then
      hh = head
      head_below_bot = .false.
      rdensegwf = this%denseterms(2, n) ! gwf rel density
    else
      hh = bots
      head_below_bot = .true.
      rdensegwf = this%denseterms(1, n) ! sfr rel density
    end if
    !
    ! -- todo: hack because denseterms not updated in a cf calculation
    if (rdensegwf == dzero) return
    !
    ! -- Update flow
    if (stage_below_bot .and. head_below_bot) then
      !
      ! -- flow is zero, so no terms are updated
      !
    else
      !
      ! -- calculate average relative density
      rdenseavg = dhalf * (rdensesfr + rdensegwf)
      !
      ! -- Add contribution of first density term:
      !      cond * (denseavg/denseref - 1) * (hgwf - hsfr)
      d1 = cond * (rdenseavg - done)
      gwfhcof = gwfhcof - d1
      gwfrhs = gwfrhs - d1 * ss
      d1 = d1 * (hh - ss)
      flow = flow + d1
      !
      ! -- Add second density term if stage and head not below bottom
      if (.not. stage_below_bot .and. .not. head_below_bot) then
        !
        ! -- Add contribution of second density term:
        !      cond * (havg - elevavg) * (densegwf - densesfr) / denseref
        elevgwf = this%denseterms(3, n)
        elevsfr = bots
        elevavg = dhalf * (elevsfr + elevgwf)
        havg = dhalf * (hh + ss)
        d2 = cond * (havg - elevavg) * (rdensegwf - rdensesfr)
        gwfrhs = gwfrhs + d2
        flow = flow + d2
      end if
    end if
    !
    ! -- return
    return

sfr_cc()

subroutine sfrmodule::sfr_cc ( class(sfrtype), intent(inout)  this, integer(i4b), intent(in)  innertot, integer(i4b), intent(in)  kiter, integer(i4b), intent(in)  iend, integer(i4b), intent(in)  icnvgmod, character(len=lenpakloc), intent(inout)  cpak, integer(i4b), intent(inout)  ipak, real(dp), intent(inout)  dpak  )

private

Perform additional convergence checks on the flow between the SFR package and the model it is attached to.

Parameters

[in,out]	this	SfrType object
[in]	innertot	total number of inner iterations
[in]	kiter	Picard iteration number
[in]	iend	flag indicating if this is the last Picard iteration
[in]	icnvgmod	flag inficating if the model has met specific convergence criteria
[in,out]	cpak	string for user node
[in,out]	ipak	location of the maximum dependent variable change
[in,out]	dpak	maximum dependent variable change

Definition at line 2162 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: totim, kstp, kper, delt
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    integer(I4B), intent(in) :: innertot !< total number of inner iterations
    integer(I4B), intent(in) :: kiter !< Picard iteration number
    integer(I4B), intent(in) :: iend !< flag indicating if this is the last Picard iteration
    integer(I4B), intent(in) :: icnvgmod !< flag inficating if the model has met specific convergence criteria
    character(len=LENPAKLOC), intent(inout) :: cpak !< string for user node
    integer(I4B), intent(inout) :: ipak !< location of the maximum dependent variable change
    real(DP), intent(inout) :: dpak !< maximum dependent variable change
    ! -- local variables
    character(len=LENPAKLOC) :: cloc
    character(len=LINELENGTH) :: tag
    integer(I4B) :: icheck
    integer(I4B) :: ipakfail
    integer(I4B) :: locdhmax
    integer(I4B) :: locrmax
    integer(I4B) :: locdqfrommvrmax
    integer(I4B) :: ntabrows
    integer(I4B) :: ntabcols
    integer(I4B) :: n
    real(DP) :: q
    real(DP) :: q0
    real(DP) :: qtolfact
    real(DP) :: dh
    real(DP) :: r
    real(DP) :: dhmax
    real(DP) :: rmax
    real(DP) :: dqfrommvr
    real(DP) :: dqfrommvrmax
    !
    ! -- initialize local variables
    icheck = this%iconvchk
    ipakfail = 0
    locdhmax = 0
    locrmax = 0
    r = dzero
    dhmax = dzero
    rmax = dzero
    locdqfrommvrmax = 0
    dqfrommvrmax = dzero
    !
    ! -- if not saving package convergence data on check convergence if
    !    the model is considered converged
    if (this%ipakcsv == 0) then
      if (icnvgmod == 0) then
        icheck = 0
      end if
      !
      ! -- saving package convergence data
    else
      !
      ! -- header for package csv
      if (.not. associated(this%pakcsvtab)) then
        !
        ! -- determine the number of columns and rows
        ntabrows = 1
        ntabcols = 9
        if (this%imover == 1) then
          ntabcols = ntabcols + 2
        end if
        !
        ! -- setup table
        call table_cr(this%pakcsvtab, this%packName, '')
        call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
                                     lineseparator=.false., separator=',', &
                                     finalize=.false.)
        !
        ! -- add columns to package csv
        tag = 'total_inner_iterations'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'totim'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'kper'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'kstp'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'nouter'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'dvmax'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'dvmax_loc'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'dinflowmax'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'dinflowmax_loc'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        if (this%imover == 1) then
          tag = 'dqfrommvrmax'
          call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
          tag = 'dqfrommvrmax_loc'
          call this%pakcsvtab%initialize_column(tag, 16, alignment=tableft)
        end if
      end if
    end if
    !
    ! -- perform package convergence check
    if (icheck /= 0) then
      final_check: do n = 1, this%maxbound
        if (this%iboundpak(n) == 0) cycle
        !
        ! -- set the Q to length factor
        qtolfact = delt / this%calc_surface_area(n)
        !
        ! -- calculate stage change
        dh = this%stage0(n) - this%stage(n)
        !
        ! -- evaluate flow difference if the time step is transient
        if (this%gwfiss == 0) then
          r = this%usflow0(n) - this%usflow(n)
          !
          ! -- normalize flow difference and convert to a depth
          r = r * qtolfact
        end if
        !
        ! -- q from mvr
        dqfrommvr = dzero
        if (this%imover == 1) then
          q = this%pakmvrobj%get_qfrommvr(n)
          q0 = this%pakmvrobj%get_qfrommvr0(n)
          dqfrommvr = qtolfact * (q0 - q)
        end if
        !
        ! -- evaluate magnitude of differences
        if (n == 1) then
          locdhmax = n
          dhmax = dh
          locrmax = n
          rmax = r
          dqfrommvrmax = dqfrommvr
          locdqfrommvrmax = n
        else
          if (abs(dh) > abs(dhmax)) then
            locdhmax = n
            dhmax = dh
          end if
          if (abs(r) > abs(rmax)) then
            locrmax = n
            rmax = r
          end if
          if (abs(dqfrommvr) > abs(dqfrommvrmax)) then
            dqfrommvrmax = dqfrommvr
            locdqfrommvrmax = n
          end if
        end if
      end do final_check
      !
      ! -- set dpak and cpak
      if (abs(dhmax) > abs(dpak)) then
        ipak = locdhmax
        dpak = dhmax
        write (cloc, "(a,'-',a)") trim(this%packName), 'stage'
        cpak = trim(cloc)
      end if
      if (abs(rmax) > abs(dpak)) then
        ipak = locrmax
        dpak = rmax
        write (cloc, "(a,'-',a)") trim(this%packName), 'inflow'
        cpak = trim(cloc)
      end if
      if (this%imover == 1) then
        if (abs(dqfrommvrmax) > abs(dpak)) then
          ipak = locdqfrommvrmax
          dpak = dqfrommvrmax
          write (cloc, "(a,'-',a)") trim(this%packName), 'qfrommvr'
          cpak = trim(cloc)
        end if
      end if
      !
      ! -- write convergence data to package csv
      if (this%ipakcsv /= 0) then
        !
        ! -- write the data
        call this%pakcsvtab%add_term(innertot)
        call this%pakcsvtab%add_term(totim)
        call this%pakcsvtab%add_term(kper)
        call this%pakcsvtab%add_term(kstp)
        call this%pakcsvtab%add_term(kiter)
        call this%pakcsvtab%add_term(dhmax)
        call this%pakcsvtab%add_term(locdhmax)
        call this%pakcsvtab%add_term(rmax)
        call this%pakcsvtab%add_term(locrmax)
        if (this%imover == 1) then
          call this%pakcsvtab%add_term(dqfrommvrmax)
          call this%pakcsvtab%add_term(locdqfrommvrmax)
        end if
        !
        ! -- finalize the package csv
        if (iend == 1) then
          call this%pakcsvtab%finalize_table()
        end if
      end if
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_cf()

subroutine sfrmodule::sfr_cf

(

class(sfrtype)

this

)

Formulate the hcof and rhs terms for the WEL package that will be added to the coefficient matrix and right-hand side vector.

Parameters

this	SfrType object

Definition at line 1972 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: igwfnode
    !
    ! -- return if no sfr reaches
    if (this%nbound == 0) return
    !
    ! -- find highest active cell
    do n = 1, this%nbound
      igwfnode = this%igwftopnode(n)
      if (igwfnode > 0) then
        if (this%ibound(igwfnode) == 0) then
          call this%dis%highest_active(igwfnode, this%ibound)
        end if
      end if
      this%igwfnode(n) = igwfnode
      this%nodelist(n) = igwfnode
    end do
    !
    ! -- return
    return

sfr_check_connections()

subroutine sfrmodule::sfr_check_connections ( class(sfrtype)  this )

private

Method to check connection data for a SFR package. This method also creates the tables used to print input data, if this option in enabled in the SFR package.

Parameters

this	SfrType object

Definition at line 4617 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    logical(LGP) :: lreorder
    character(len=5) :: crch
    character(len=5) :: crch2
    character(len=LINELENGTH) :: text
    character(len=LINELENGTH) :: title
    integer(I4B) :: n
    integer(I4B) :: nn
    integer(I4B) :: nc
    integer(I4B) :: i
    integer(I4B) :: ii
    integer(I4B) :: j
    integer(I4B) :: ifound
    integer(I4B) :: ierr
    integer(I4B) :: maxconn
    integer(I4B) :: ntabcol
    !
    ! -- determine if the reaches have been reordered
    lreorder = .false.
    do j = 1, this%MAXBOUND
      n = this%isfrorder(j)
      if (n /= j) then
        lreorder = .true.
        exit
      end if
    end do
    !
    ! -- write message that the solution order h
    if (lreorder) then
      write (this%iout, '(/,1x,a)') &
        trim(adjustl(this%text))//' PACKAGE ('// &
        trim(adjustl(this%packName))//') REACH SOLUTION HAS BEEN '// &
        'REORDERED USING A DAG'
      !
      ! -- print table
      if (this%iprpak /= 0) then
        !
        ! -- reset the input table object
        ntabcol = 2
        title = trim(adjustl(this%text))//' PACKAGE ('// &
                trim(adjustl(this%packName))//') REACH SOLUTION ORDER'
        call table_cr(this%inputtab, this%packName, title)
        call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
        text = 'ORDER'
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
        text = 'REACH'
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
        !
        ! -- upstream connection data
        do j = 1, this%maxbound
          n = this%isfrorder(j)
          call this%inputtab%add_term(j)
          call this%inputtab%add_term(n)
        end do
      end if
    end if
    !
    ! -- create input table for reach connections data
    if (this%iprpak /= 0) then
      !
      ! -- calculate the maximum number of connections
      maxconn = 0
      do n = 1, this%maxbound
        maxconn = max(maxconn, this%nconnreach(n))
      end do
      ntabcol = 1 + maxconn
      !
      ! -- reset the input table object
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STATIC REACH CONNECTION DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
      text = 'REACH'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      do n = 1, maxconn
        write (text, '(a,1x,i6)') 'CONN', n
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      end do
    end if
    !
    ! -- check the reach connections for simple errors
    ! -- connection check
    do n = 1, this%MAXBOUND
      write (crch, '(i5)') n
      eachconn: do i = this%ia(n) + 1, this%ia(n + 1) - 1
        nn = this%ja(i)
        write (crch2, '(i5)') nn
        ifound = 0
        connreach: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
          nc = this%ja(ii)
          if (nc == n) then
            ifound = 1
            exit connreach
          end if
        end do connreach
        if (ifound /= 1) then
          errmsg = 'Reach '//crch//' is connected to '// &
                   'reach '//crch2//' but reach '//crch2// &
                   ' is not connected to reach '//crch//'.'
          call store_error(errmsg)
        end if
      end do eachconn
      !
      ! -- write connection data to the table
      if (this%iprpak /= 0) then
        call this%inputtab%add_term(n)
        do i = this%ia(n) + 1, this%ia(n + 1) - 1
          call this%inputtab%add_term(this%ja(i))
        end do
        nn = maxconn - this%nconnreach(n)
        do i = 1, nn
          call this%inputtab%add_term(' ')
        end do
      end if
    end do
    !
    ! -- check for incorrect connections between upstream connections
    !
    ! -- check upstream connections for each reach
    ierr = 0
    do n = 1, this%maxbound
      write (crch, '(i5)') n
      eachconnv: do i = this%ia(n) + 1, this%ia(n + 1) - 1
        !
        ! -- skip downstream connections
        if (this%idir(i) < 0) cycle eachconnv
        nn = this%ja(i)
        write (crch2, '(i5)') nn
        connreachv: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
          ! -- skip downstream connections
          if (this%idir(ii) < 0) cycle connreachv
          nc = this%ja(ii)
          !
          ! -- if nc == n then that means reach n is an upstream connection for
          !    reach nn and reach nn is an upstream connection for reach n
          if (nc == n) then
            ierr = ierr + 1
            errmsg = 'Reach '//crch//' is connected to '// &
                     'reach '//crch2//' but streamflow from reach '// &
                     crch//' to reach '//crch2//' is not permitted.'
            call store_error(errmsg)
            exit connreachv
          end if
        end do connreachv
      end do eachconnv
    end do
    !
    ! -- terminate if connectivity errors
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- check that downstream reaches for a reach are
    !    the upstream reaches for the reach
    do n = 1, this%maxbound
      write (crch, '(i5)') n
      eachconnds: do i = this%ia(n) + 1, this%ia(n + 1) - 1
        nn = this%ja(i)
        if (this%idir(i) > 0) cycle eachconnds
        write (crch2, '(i5)') nn
        ifound = 0
        connreachds: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
          nc = this%ja(ii)
          if (nc == n) then
            if (this%idir(i) /= this%idir(ii)) then
              ifound = 1
            end if
            exit connreachds
          end if
        end do connreachds
        if (ifound /= 1) then
          errmsg = 'Reach '//crch//' downstream connected reach '// &
                   'is reach '//crch2//' but reach '//crch//' is not'// &
                   ' the upstream connected reach for reach '//crch2//'.'
          call store_error(errmsg)
        end if
      end do eachconnds
    end do
    !
    ! -- create input table for upstream and downstream connections
    if (this%iprpak /= 0) then
      !
      ! -- calculate the maximum number of upstream connections
      maxconn = 0
      do n = 1, this%maxbound
        ii = 0
        do i = this%ia(n) + 1, this%ia(n + 1) - 1
          if (this%idir(i) > 0) then
            ii = ii + 1
          end if
        end do
        maxconn = max(maxconn, ii)
      end do
      ntabcol = 1 + maxconn
      !
      ! -- reset the input table object
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STATIC UPSTREAM REACH '// &
              'CONNECTION DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
      text = 'REACH'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      do n = 1, maxconn
        write (text, '(a,1x,i6)') 'UPSTREAM CONN', n
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      end do
      !
      ! -- upstream connection data
      do n = 1, this%maxbound
        call this%inputtab%add_term(n)
        ii = 0
        do i = this%ia(n) + 1, this%ia(n + 1) - 1
          if (this%idir(i) > 0) then
            call this%inputtab%add_term(this%ja(i))
            ii = ii + 1
          end if
        end do
        nn = maxconn - ii
        do i = 1, nn
          call this%inputtab%add_term(' ')
        end do
      end do
      !
      ! -- calculate the maximum number of downstream connections
      maxconn = 0
      do n = 1, this%maxbound
        ii = 0
        do i = this%ia(n) + 1, this%ia(n + 1) - 1
          if (this%idir(i) < 0) then
            ii = ii + 1
          end if
        end do
        maxconn = max(maxconn, ii)
      end do
      ntabcol = 1 + maxconn
      !
      ! -- reset the input table object
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STATIC DOWNSTREAM '// &
              'REACH CONNECTION DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
      text = 'REACH'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      do n = 1, maxconn
        write (text, '(a,1x,i6)') 'DOWNSTREAM CONN', n
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      end do
      !
      ! -- downstream connection data
      do n = 1, this%maxbound
        call this%inputtab%add_term(n)
        ii = 0
        do i = this%ia(n) + 1, this%ia(n + 1) - 1
          if (this%idir(i) < 0) then
            call this%inputtab%add_term(this%ja(i))
            ii = ii + 1
          end if
        end do
        nn = maxconn - ii
        do i = 1, nn
          call this%inputtab%add_term(' ')
        end do
      end do
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_check_conversion()

subroutine sfrmodule::sfr_check_conversion ( class(sfrtype)  this )

private

Method to check unit conversion data for a SFR package. This method also calculates unitconv that is used in the Manning's equation.

Parameters

this	SfrType object

Definition at line 4452 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    ! -- formats
    character(len=*), parameter :: fmtunitconv_error = &
      &"('SFR (',a,') UNIT_CONVERSION SPECIFIED VALUE (',g0,') AND', &
      &1x,'LENGTH_CONVERSION OR TIME_CONVERSION SPECIFIED.')"
    character(len=*), parameter :: fmtunitconv = &
      &"(1x,'SFR PACKAGE (',a,') CONVERSION DATA',&
      &/4x,'UNIT CONVERSION VALUE (',g0,').',/)"
    !
    ! -- check the reach data for simple errors
    if (this%lengthconv /= dnodata .or. this%timeconv /= dnodata) then
      if (this%unitconv /= done) then
        write (errmsg, fmtunitconv_error) &
          trim(adjustl(this%packName)), this%unitconv
        call store_error(errmsg)
      else
        if (this%lengthconv /= dnodata) then
          this%unitconv = this%unitconv * this%lengthconv**donethird
        end if
        if (this%timeconv /= dnodata) then
          this%unitconv = this%unitconv * this%timeconv
        end if
        write (this%iout, fmtunitconv) &
          trim(adjustl(this%packName)), this%unitconv
      end if
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_check_diversions()

subroutine sfrmodule::sfr_check_diversions ( class(sfrtype)  this )

private

Method to check diversion data for a SFR package. This method also creates the tables used to print input data, if this option in enabled in the SFR package.

Parameters

this	SfrType object

Definition at line 4898 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    character(len=5) :: crch
    character(len=5) :: cdiv
    character(len=5) :: crch2
    character(len=10) :: cprior
    integer(I4B) :: maxdiv
    integer(I4B) :: n
    integer(I4B) :: nn
    integer(I4B) :: nc
    integer(I4B) :: ii
    integer(I4B) :: idiv
    integer(I4B) :: ifound
    integer(I4B) :: jpos
    ! -- format
10  format('Diversion ', i0, ' of reach ', i0, &
           ' is invalid or has not been defined.')
    !
    ! -- write header
    if (this%iprpak /= 0) then
      !
      ! -- determine the maximum number of diversions
      maxdiv = 0
      do n = 1, this%maxbound
        maxdiv = maxdiv + this%ndiv(n)
      end do
      !
      ! -- reset the input table object
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') REACH DIVERSION DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(maxdiv, 4, this%iout)
      text = 'REACH'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'DIVERSION'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'REACH 2'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'CPRIOR'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
    end if
    !
    ! -- check that diversion data are correct
    do n = 1, this%maxbound
      if (this%ndiv(n) < 1) cycle
      write (crch, '(i5)') n
 
      do idiv = 1, this%ndiv(n)
        !
        ! -- determine diversion index
        jpos = this%iadiv(n) + idiv - 1
        !
        ! -- write idiv to cdiv
        write (cdiv, '(i5)') idiv
        !
        !
        nn = this%divreach(jpos)
        write (crch2, '(i5)') nn
        !
        ! -- make sure diversion reach is connected to current reach
        ifound = 0
        if (nn < 1 .or. nn > this%maxbound) then
          write (errmsg, 10) idiv, n
          call store_error(errmsg)
          cycle
        end if
        connreach: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
          nc = this%ja(ii)
          if (nc == n) then
            if (this%idir(ii) > 0) then
              ifound = 1
            end if
            exit connreach
          end if
        end do connreach
        if (ifound /= 1) then
          errmsg = 'Reach '//crch//' is not a upstream reach for '// &
                   'reach '//crch2//' as a result diversion '//cdiv// &
                   ' from reach '//crch//' to reach '//crch2// &
                   ' is not possible. Check reach connectivity.'
          call store_error(errmsg)
        end if
        ! -- iprior
        cprior = this%divcprior(jpos)
        !
        ! -- add terms to the table
        if (this%iprpak /= 0) then
          call this%inputtab%add_term(n)
          call this%inputtab%add_term(idiv)
          call this%inputtab%add_term(nn)
          call this%inputtab%add_term(cprior)
        end if
      end do
    end do
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_check_reaches()

subroutine sfrmodule::sfr_check_reaches ( class(sfrtype)  this )

private

Method to check specified data for a SFR package. This method also creates the tables used to print input data, if this option in enabled in the SFR package.

Parameters

this	SfrType object

Definition at line 4493 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    character(len=5) :: crch
    character(len=10) :: cval
    character(len=30) :: nodestr
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    integer(I4B) :: n
    integer(I4B) :: nn
    real(DP) :: btgwf
    real(DP) :: bt
    !
    ! -- setup inputtab tableobj
    if (this%iprpak /= 0) then
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STATIC REACH DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%maxbound, 10, this%iout)
      text = 'NUMBER'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'CELLID'
      call this%inputtab%initialize_column(text, 20, alignment=tableft)
      text = 'LENGTH'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      text = 'WIDTH'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      text = 'SLOPE'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      text = 'TOP'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      text = 'THICKNESS'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      text = 'HK'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      text = 'ROUGHNESS'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      text = 'UPSTREAM FRACTION'
      call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
    end if
    !
    ! --
    !
    ! -- check the reach data for simple errors
    do n = 1, this%maxbound
      write (crch, '(i5)') n
      nn = this%igwfnode(n)
      if (nn > 0) then
        btgwf = this%dis%bot(nn)
        call this%dis%noder_to_string(nn, nodestr)
      else
        nodestr = 'none'
      end if
      ! -- check reach length
      if (this%length(n) <= dzero) then
        errmsg = 'Reach '//crch//' length must be greater than 0.0.'
        call store_error(errmsg)
      end if
      ! -- check reach width
      if (this%width(n) <= dzero) then
        errmsg = 'Reach '//crch//' width must be greater than 0.0.'
        call store_error(errmsg)
      end if
      ! -- check reach slope
      if (this%slope(n) <= dzero) then
        errmsg = 'Reach '//crch//' slope must be greater than 0.0.'
        call store_error(errmsg)
      end if
      ! -- check bed thickness and bed hk for reaches connected to GWF
      if (nn > 0) then
        bt = this%strtop(n) - this%bthick(n)
        if (bt <= btgwf .and. this%icheck /= 0) then
          write (cval, '(f10.4)') bt
          errmsg = 'Reach '//crch//' bed bottom (rtp-rbth ='// &
                   cval//') must be greater than the bottom of cell ('// &
                   nodestr
          write (cval, '(f10.4)') btgwf
          errmsg = trim(adjustl(errmsg))//'='//cval//').'
          call store_error(errmsg)
        end if
        if (this%hk(n) < dzero) then
          errmsg = 'Reach '//crch//' hk must be greater than or equal to 0.0.'
          call store_error(errmsg)
        end if
      end if
      ! -- check reach roughness
      if (this%rough(n) <= dzero) then
        errmsg = 'Reach '//crch//" Manning's roughness "// &
                 'coefficient must be greater than 0.0.'
        call store_error(errmsg)
      end if
      ! -- check reach upstream fraction
      if (this%ustrf(n) < dzero) then
        errmsg = 'Reach '//crch//' upstream fraction must be greater '// &
                 'than or equal to 0.0.'
        call store_error(errmsg)
      end if
      ! -- write summary of reach information
      if (this%iprpak /= 0) then
        call this%inputtab%add_term(n)
        call this%inputtab%add_term(nodestr)
        call this%inputtab%add_term(this%length(n))
        call this%inputtab%add_term(this%width(n))
        call this%inputtab%add_term(this%slope(n))
        call this%inputtab%add_term(this%strtop(n))
        call this%inputtab%add_term(this%bthick(n))
        call this%inputtab%add_term(this%hk(n))
        call this%inputtab%add_term(this%rough(n))
        call this%inputtab%add_term(this%ustrf(n))
      end if
    end do
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_check_ustrf()

subroutine sfrmodule::sfr_check_ustrf ( class(sfrtype)  this )

private

Method to check upstream fraction data for a SFR package. This method also creates the tables used to print input data, if this option in enabled in the SFR package.

Parameters

this	SfrType object

Definition at line 5008 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    logical(LGP) :: lcycle
    logical(LGP) :: ladd
    character(len=5) :: crch
    character(len=5) :: crch2
    character(len=10) :: cval
    integer(I4B) :: maxcols
    integer(I4B) :: npairs
    integer(I4B) :: ipair
    integer(I4B) :: i
    integer(I4B) :: n
    integer(I4B) :: n2
    integer(I4B) :: idiv
    integer(I4B) :: i0
    integer(I4B) :: i1
    integer(I4B) :: jpos
    integer(I4B) :: ids
    real(DP) :: f
    real(DP) :: rval
    !
    ! -- write table header
    if (this%iprpak /= 0) then
      !
      ! -- determine the maximum number of columns
      npairs = 0
      do n = 1, this%maxbound
        ipair = 0
        ec: do i = this%ia(n) + 1, this%ia(n + 1) - 1
          !
          ! -- skip upstream connections
          if (this%idir(i) > 0) cycle ec
          n2 = this%ja(i)
          !
          ! -- skip inactive downstream reaches
          if (this%iboundpak(n2) == 0) cycle ec
          !
          ! -- increment ipair and see if it exceeds npairs
          ipair = ipair + 1
          npairs = max(npairs, ipair)
        end do ec
      end do
      maxcols = 1 + npairs * 2
      !
      ! -- reset the input table object
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') CONNECTED REACH UPSTREAM '// &
              'FRACTION DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%maxbound, maxcols, this%iout)
      text = 'REACH'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      do i = 1, npairs
        write (cval, '(i10)') i
        text = 'DOWNSTREAM REACH '//trim(adjustl(cval))
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
        text = 'FRACTION '//trim(adjustl(cval))
        call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
      end do
    end if
    !
    ! -- fill diversion number for each connection
    do n = 1, this%maxbound
      do idiv = 1, this%ndiv(n)
        i0 = this%iadiv(n)
        i1 = this%iadiv(n + 1) - 1
        do jpos = i0, i1
          do i = this%ia(n) + 1, this%ia(n + 1) - 1
            n2 = this%ja(i)
            if (this%divreach(jpos) == n2) then
              this%idiv(i) = jpos - i0 + 1
              exit
            end if
          end do
        end do
      end do
    end do
    !
    ! -- check that the upstream fraction for reach connected by
    !    a diversion is zero
    do n = 1, this%maxbound
      !
      ! -- determine the number of downstream reaches
      ids = 0
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(i) < 0) then
          ids = ids + 1
        end if
      end do
      !
      ! -- evaluate the diversions
      do idiv = 1, this%ndiv(n)
        jpos = this%iadiv(n) + idiv - 1
        n2 = this%divreach(jpos)
        f = this%ustrf(n2)
        if (f /= dzero) then
          write (errmsg, '(a,2(1x,i0,1x,a),1x,a,g0,a,2(1x,a))') &
            'Reach', n, 'is connected to reach', n2, 'by a diversion', &
            'but the upstream fraction is not equal to zero (', f, '). Check', &
            trim(this%packName), 'package diversion and package data.'
          if (ids > 1) then
            call store_error(errmsg)
          else
            write (warnmsg, '(a,3(1x,a))') &
              trim(warnmsg), &
              'A warning instead of an error is issued because', &
              'the reach is only connected to the diversion reach in the ', &
              'downstream direction.'
            call store_warning(warnmsg)
          end if
        end if
      end do
    end do
    !
    ! -- calculate the total fraction of connected reaches that are
    !    not diversions and check that the sum of upstream fractions
    !    is equal to 1 for each reach
    do n = 1, this%maxbound
      ids = 0
      rval = dzero
      f = dzero
      write (crch, '(i5)') n
      if (this%iprpak /= 0) then
        call this%inputtab%add_term(n)
      end if
      ipair = 0
      eachconn: do i = this%ia(n) + 1, this%ia(n + 1) - 1
        lcycle = .false.
        !
        ! -- initialize downstream connection q
        this%qconn(i) = dzero
        !
        ! -- skip upstream connections
        if (this%idir(i) > 0) then
          lcycle = .true.
        end if
        n2 = this%ja(i)
        !
        ! -- skip inactive downstream reaches
        if (this%iboundpak(n2) == 0) then
          lcycle = .true.
        end if
        if (lcycle) then
          cycle eachconn
        end if
        ipair = ipair + 1
        write (crch2, '(i5)') n2
        ids = ids + 1
        ladd = .true.
        f = f + this%ustrf(n2)
        write (cval, '(f10.4)') this%ustrf(n2)
        !
        ! -- write upstream fractions
        if (this%iprpak /= 0) then
          call this%inputtab%add_term(n2)
          call this%inputtab%add_term(this%ustrf(n2))
        end if
        eachdiv: do idiv = 1, this%ndiv(n)
          jpos = this%iadiv(n) + idiv - 1
          if (this%divreach(jpos) == n2) then
            ladd = .false.
            exit eachdiv
          end if
        end do eachdiv
        if (ladd) then
          rval = rval + this%ustrf(n2)
        end if
      end do eachconn
      this%ftotnd(n) = rval
      !
      ! -- write remaining table columns
      if (this%iprpak /= 0) then
        ipair = ipair + 1
        do i = ipair, npairs
          call this%inputtab%add_term('  ')
          call this%inputtab%add_term('  ')
        end do
      end if
      !
      ! -- evaluate if an error condition has occurred
      !    the sum of fractions is not equal to 1
      if (ids /= 0) then
        if (abs(f - done) > dem6) then
          write (errmsg, '(a,1x,i0,1x,a,g0,a,3(1x,a))') &
            'Upstream fractions for reach ', n, 'is not equal to one (', f, &
            '). Check', trim(this%packName), 'package reach connectivity and', &
            'package data.'
          call store_error(errmsg)
        end if
      end if
    end do
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_cq()

subroutine sfrmodule::sfr_cq	(	class(sfrtype), intent(inout)	this,
		real(dp), dimension(:), intent(in)	x,
		real(dp), dimension(:), intent(inout), contiguous	flowja,
		integer(i4b), intent(in), optional	iadv
	)

Calculate the flow between connected SFR package control volumes.

Parameters

[in,out]	this	SfrType object
[in]	x	current dependent-variable value
[in,out]	flowja	flow between two connected control volumes
[in]	iadv	flag that indicates if this is an advance package

Definition at line 2367 of file gwf-sfr.f90.

    ! -- modules
    use budgetmodule, only: budgettype
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    real(DP), dimension(:), intent(in) :: x !< current dependent-variable value
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
    integer(I4B), optional, intent(in) :: iadv !< flag that indicates if this is an advance package
    ! -- local variables
    integer(I4B) :: i
    real(DP) :: qext
    ! -- for budget
    integer(I4B) :: n
    integer(I4B) :: n2
    real(DP) :: qoutflow
    real(DP) :: qfrommvr
    real(DP) :: qtomvr
    !
    ! -- call base functionality in bnd_cq.  This will calculate sfr-gwf flows
    !    and put them into this%simvals
    call this%BndType%bnd_cq(x, flowja, iadv=1)
    !
    ! -- Calculate qextoutflow and qoutflow for subsequent budgets
    do n = 1, this%maxbound
      !
      ! -- mover
      qfrommvr = dzero
      qtomvr = dzero
      if (this%imover == 1) then
        qfrommvr = this%pakmvrobj%get_qfrommvr(n)
        qtomvr = this%pakmvrobj%get_qtomvr(n)
        if (qtomvr > dzero) then
          qtomvr = -qtomvr
        end if
      end if
      !
      ! -- external downstream stream flow
      qext = this%dsflow(n)
      qoutflow = dzero
      if (qext > dzero) then
        qext = -qext
      end if
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(i) > 0) cycle
        n2 = this%ja(i)
        if (this%iboundpak(n2) == 0) cycle
        qext = dzero
        exit
      end do
      !
      ! -- adjust external downstream stream flow using qtomvr
      if (qext < dzero) then
        if (qtomvr < dzero) then
          qext = qext - qtomvr
        end if
      else
        qoutflow = this%dsflow(n)
        if (qoutflow > dzero) then
          qoutflow = -qoutflow
        end if
      end if
      !
      ! -- set qextoutflow and qoutflow for cell by cell budget
      !    output and observations
      this%qextoutflow(n) = qext
      this%qoutflow(n) = qoutflow
      !
    end do
    !
    ! -- fill the budget object
    call this%sfr_fill_budobj()
    !
    ! -- return
    return

sfr_create()

subroutine, public sfrmodule::sfr_create	(	class(bndtype), pointer	packobj,
		integer(i4b), intent(in)	id,
		integer(i4b), intent(in)	ibcnum,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout,
		character(len=*), intent(in)	namemodel,
		character(len=*), intent(in)	pakname
	)

Create a new SFR Package object

Parameters

	packobj	pointer to default package type
[in]	id	package id
[in]	ibcnum	boundary condition number
[in]	inunit	unit number of SFR package input file
[in]	iout	unit number of model listing file
[in]	namemodel	model name
[in]	pakname	package name

Definition at line 235 of file gwf-sfr.f90.

    ! -- modules
    use memoryhelpermodule, only: create_mem_path
    ! -- dummy variables
    class(BndType), pointer :: packobj !< pointer to default package type
    integer(I4B), intent(in) :: id !< package id
    integer(I4B), intent(in) :: ibcnum !< boundary condition number
    integer(I4B), intent(in) :: inunit !< unit number of SFR package input file
    integer(I4B), intent(in) :: iout !< unit number of model listing file
    character(len=*), intent(in) :: namemodel !< model name
    character(len=*), intent(in) :: pakname !< package name
    ! -- local variables
    type(SfrType), pointer :: sfrobj
    !
    ! -- allocate the object and assign values to object variables
    allocate (sfrobj)
    packobj => sfrobj
    !
    ! -- create name and memory path
    call packobj%set_names(ibcnum, namemodel, pakname, ftype)
    packobj%text = text
    !
    ! -- allocate scalars
    call sfrobj%sfr_allocate_scalars()
    !
    ! -- initialize package
    call packobj%pack_initialize()
 
    packobj%inunit = inunit
    packobj%iout = iout
    packobj%id = id
    packobj%ibcnum = ibcnum
    packobj%ncolbnd = 4
    packobj%iscloc = 0 ! not supported
    packobj%isadvpak = 1
    packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
    !
    ! -- return
    return

Here is the call graph for this function:

Here is the caller graph for this function:

sfr_da()

subroutine sfrmodule::sfr_da

(

class(sfrtype)

this

)

Deallocate SFR package scalars and arrays.

Parameters

this	SfrType object

Definition at line 2642 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    !
    ! -- deallocate arrays
    call mem_deallocate(this%qoutflow)
    call mem_deallocate(this%qextoutflow)
    deallocate (this%csfrbudget)
    call mem_deallocate(this%sfrname, 'SFRNAME', this%memoryPath)
    call mem_deallocate(this%dbuff)
    deallocate (this%cauxcbc)
    call mem_deallocate(this%qauxcbc)
    call mem_deallocate(this%iboundpak)
    call mem_deallocate(this%igwfnode)
    call mem_deallocate(this%igwftopnode)
    call mem_deallocate(this%length)
    call mem_deallocate(this%width)
    call mem_deallocate(this%strtop)
    call mem_deallocate(this%bthick)
    call mem_deallocate(this%hk)
    call mem_deallocate(this%slope)
    call mem_deallocate(this%nconnreach)
    call mem_deallocate(this%ustrf)
    call mem_deallocate(this%ftotnd)
    call mem_deallocate(this%usflow)
    call mem_deallocate(this%dsflow)
    call mem_deallocate(this%depth)
    call mem_deallocate(this%stage)
    call mem_deallocate(this%gwflow)
    call mem_deallocate(this%simevap)
    call mem_deallocate(this%simrunoff)
    call mem_deallocate(this%stage0)
    call mem_deallocate(this%usflow0)
    call mem_deallocate(this%denseterms)
    call mem_deallocate(this%viscratios)
    !
    ! -- deallocate reach order and connection data
    call mem_deallocate(this%isfrorder)
    call mem_deallocate(this%ia)
    call mem_deallocate(this%ja)
    call mem_deallocate(this%idir)
    call mem_deallocate(this%idiv)
    call mem_deallocate(this%qconn)
    !
    ! -- deallocate boundary data
    call mem_deallocate(this%rough)
    call mem_deallocate(this%rain)
    call mem_deallocate(this%evap)
    call mem_deallocate(this%inflow)
    call mem_deallocate(this%runoff)
    call mem_deallocate(this%sstage)
    !
    ! -- deallocate aux variables
    call mem_deallocate(this%rauxvar)
    !
    ! -- deallocate diversion variables
    call mem_deallocate(this%iadiv)
    call mem_deallocate(this%divreach)
    if (associated(this%divcprior)) then
      deallocate (this%divcprior)
    end if
    call mem_deallocate(this%divflow)
    call mem_deallocate(this%divq)
    call mem_deallocate(this%ndiv)
    !
    ! -- deallocate cross-section data
    call mem_deallocate(this%ncrosspts)
    call mem_deallocate(this%iacross)
    call mem_deallocate(this%station)
    call mem_deallocate(this%xsheight)
    call mem_deallocate(this%xsrough)
    !
    ! -- deallocate budobj
    call this%budobj%budgetobject_da()
    deallocate (this%budobj)
    nullify (this%budobj)
    !
    ! -- deallocate stage table
    if (this%iprhed > 0) then
      call this%stagetab%table_da()
      deallocate (this%stagetab)
      nullify (this%stagetab)
    end if
    !
    ! -- deallocate package csv table
    if (this%ipakcsv > 0) then
      if (associated(this%pakcsvtab)) then
        call this%pakcsvtab%table_da()
        deallocate (this%pakcsvtab)
        nullify (this%pakcsvtab)
      end if
    end if
    !
    ! -- deallocate scalars
    call mem_deallocate(this%iprhed)
    call mem_deallocate(this%istageout)
    call mem_deallocate(this%ibudgetout)
    call mem_deallocate(this%ibudcsv)
    call mem_deallocate(this%ipakcsv)
    call mem_deallocate(this%idiversions)
    call mem_deallocate(this%maxsfrpicard)
    call mem_deallocate(this%maxsfrit)
    call mem_deallocate(this%bditems)
    call mem_deallocate(this%cbcauxitems)
    call mem_deallocate(this%unitconv)
    call mem_deallocate(this%lengthconv)
    call mem_deallocate(this%timeconv)
    call mem_deallocate(this%dmaxchg)
    call mem_deallocate(this%deps)
    call mem_deallocate(this%nconn)
    call mem_deallocate(this%icheck)
    call mem_deallocate(this%iconvchk)
    call mem_deallocate(this%idense)
    call mem_deallocate(this%ianynone)
    call mem_deallocate(this%ncrossptstot)
    nullify (this%gwfiss)
    !
    ! -- call base BndType deallocate
    call this%BndType%bnd_da()
    !
    ! -- return
    return

sfr_df_obs()

subroutine sfrmodule::sfr_df_obs ( class(sfrtype)  this )

private

Method to define the observation types available in the SFR package.

Parameters

this	SfrType object

Definition at line 2826 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: indx
    !
    ! -- Store obs type and assign procedure pointer
    !    for stage observation type.
    call this%obs%StoreObsType('stage', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for inflow observation type.
    call this%obs%StoreObsType('inflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for inflow observation type.
    call this%obs%StoreObsType('ext-inflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for rainfall observation type.
    call this%obs%StoreObsType('rainfall', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for runoff observation type.
    call this%obs%StoreObsType('runoff', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for evaporation observation type.
    call this%obs%StoreObsType('evaporation', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for outflow observation type.
    call this%obs%StoreObsType('outflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for ext-outflow observation type.
    call this%obs%StoreObsType('ext-outflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for to-mvr observation type.
    call this%obs%StoreObsType('to-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for sfr-frommvr observation type.
    call this%obs%StoreObsType('from-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for sfr observation type.
    call this%obs%StoreObsType('sfr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for upstream flow observation type.
    call this%obs%StoreObsType('upstream-flow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for downstream flow observation type.
    call this%obs%StoreObsType('downstream-flow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for depth observation type.
    call this%obs%StoreObsType('depth', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for wetted-perimeter observation type.
    call this%obs%StoreObsType('wet-perimeter', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for wetted-area observation type.
    call this%obs%StoreObsType('wet-area', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for wetted-width observation type.
    call this%obs%StoreObsType('wet-width', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_fc()

subroutine sfrmodule::sfr_fc ( class(sfrtype)  this, real(dp), dimension(:), intent(inout)  rhs, integer(i4b), dimension(:), intent(in)  ia, integer(i4b), dimension(:), intent(in)  idxglo, class(matrixbasetype), pointer  matrix_sln  )

private

Add the hcof and rhs terms for the SFR package to the coefficient matrix and right-hand side vector.

Parameters

	this	SfrType object
[in,out]	rhs	right-hand side vector for model
[in]	ia	solution CRS row pointers
[in]	idxglo	mapping vector for model (local) to solution (global)
	matrix_sln	solution coefficient matrix

Definition at line 2004 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model
    integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers
    integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global)
    class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: ipos
    integer(I4B) :: node
    real(DP) :: s0
    real(DP) :: ds
    real(DP) :: dsmax
    real(DP) :: hgwf
    real(DP) :: v
    real(DP) :: hhcof
    real(DP) :: rrhs
    !
    ! -- picard iterations for sfr to achieve good solution regardless
    !    of reach order
    sfrpicard: do i = 1, this%maxsfrpicard
      !
      ! -- initialize maximum stage change for iteration to zero
      dsmax = dzero
      !
      ! -- pakmvrobj fc - reset qformvr to zero
      if (this%imover == 1) then
        call this%pakmvrobj%fc()
      end if
      !
      ! -- solve for each sfr reach
      reachsolve: do j = 1, this%nbound
        n = this%isfrorder(j)
        node = this%igwfnode(n)
        if (node > 0) then
          hgwf = this%xnew(node)
        else
          hgwf = dep20
        end if
        !
        ! -- save previous stage and upstream flow
        if (i == 1) then
          this%stage0(n) = this%stage(n)
          this%usflow0(n) = this%usflow(n)
        end if
        !
        ! -- set initial stage to calculate stage change
        s0 = this%stage(n)
        !
        ! -- solve for flow in swr
        if (this%iboundpak(n) /= 0) then
          call this%sfr_solve(n, hgwf, hhcof, rrhs)
        else
          this%depth(n) = dzero
          this%stage(n) = this%strtop(n)
          v = dzero
          call this%sfr_update_flows(n, v, v)
          hhcof = dzero
          rrhs = dzero
        end if
        !
        ! -- set package hcof and rhs
        this%hcof(n) = hhcof
        this%rhs(n) = rrhs
        !
        ! -- calculate stage change
        ds = s0 - this%stage(n)
        !
        ! -- evaluate if stage change exceeds dsmax
        if (abs(ds) > abs(dsmax)) then
          dsmax = ds
        end if
 
      end do reachsolve
      !
      ! -- evaluate if the sfr picard iterations should be terminated
      if (abs(dsmax) <= this%dmaxchg) then
        exit sfrpicard
      end if
 
    end do sfrpicard
    !
    ! -- Copy package rhs and hcof into solution rhs and amat
    do n = 1, this%nbound
      node = this%nodelist(n)
      if (node < 1) cycle
      rhs(node) = rhs(node) + this%rhs(n)
      ipos = ia(node)
      call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(n))
    end do
    !
    ! -- return
    return

sfr_fill_budobj()

subroutine sfrmodule::sfr_fill_budobj ( class(sfrtype)  this )

private

Method to copy flows into the budget object that stores all the sfr flows The terms listed here must correspond in number and order to the ones added in the sfr_setup_budobj method.

Parameters

this	SfrType object

Definition at line 5434 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: naux
    integer(I4B) :: i
    integer(I4B) :: n
    integer(I4B) :: n1
    integer(I4B) :: n2
    integer(I4B) :: ii
    integer(I4B) :: idx
    integer(I4B) :: idiv
    integer(I4B) :: jpos
    real(DP) :: q
    real(DP) :: qt
    real(DP) :: d
    real(DP) :: ca
    real(DP) :: a
    real(DP) :: wp
    real(DP) :: l
    !
    ! -- initialize counter
    idx = 0
    !
    ! -- FLOW JA FACE
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nconn)
    do n = 1, this%maxbound
      n1 = n
      q = dzero
      ca = dzero
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        n2 = this%ja(i)
        if (this%iboundpak(n) /= 0) then
          ! flow to downstream reaches
          if (this%idir(i) < 0) then
            qt = this%dsflow(n)
            q = -this%qconn(i)
            ! flow from upstream reaches
          else
            qt = this%usflow(n)
            do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
              if (this%idir(ii) > 0) cycle
              if (this%ja(ii) /= n) cycle
              q = this%qconn(ii)
              exit
            end do
          end if
          ! calculate flow area
          call this%sfr_calc_reach_depth(n, qt, d)
          ca = this%calc_area_wet(n, d)
        else
          q = dzero
          ca = dzero
        end if
        this%qauxcbc(1) = ca
        call this%budobj%budterm(idx)%update_term(n1, n2, q, this%qauxcbc)
      end do
    end do
    !
    ! -- GWF (LEAKAGE)
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound - this%ianynone)
    do n = 1, this%maxbound
      n2 = this%igwfnode(n)
      if (n2 > 0) then
        if (this%iboundpak(n) /= 0) then
          ! -- calc_perimeter_wet() does not enforce depth dependence
          if (this%depth(n) > dzero) then
            wp = this%calc_perimeter_wet(n, this%depth(n))
          else
            wp = dzero
          end if
          l = this%length(n)
          a = wp * l
          this%qauxcbc(1) = a
          q = -this%gwflow(n)
        else
          this%qauxcbc(1) = dzero
          q = dzero
        end if
        call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
      end if
    end do
    !
    ! -- RAIN
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    do n = 1, this%maxbound
      if (this%iboundpak(n) /= 0) then
        a = this%calc_surface_area(n)
        q = this%rain(n) * a
      else
        q = dzero
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- EVAPORATION
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    do n = 1, this%maxbound
      if (this%iboundpak(n) /= 0) then
        q = -this%simevap(n)
      else
        q = dzero
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- RUNOFF
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    do n = 1, this%maxbound
      if (this%iboundpak(n) /= 0) then
        q = this%simrunoff(n)
      else
        q = dzero
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- INFLOW
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    do n = 1, this%maxbound
      if (this%iboundpak(n) /= 0) then
        q = this%inflow(n)
      else
        q = dzero
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- EXTERNAL OUTFLOW
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    do n = 1, this%maxbound
      q = dzero
      if (this%iboundpak(n) /= 0) then
        do i = this%ia(n) + 1, this%ia(n + 1) - 1
          if (this%idir(i) > 0) cycle
          idiv = this%idiv(i)
          if (idiv > 0) then
            jpos = this%iadiv(n) + idiv - 1
            q = q + this%divq(jpos)
          else
            q = q + this%qconn(i)
          end if
        end do
        q = q - this%dsflow(n)
        if (this%imover == 1) then
          q = q + this%pakmvrobj%get_qtomvr(n)
        end if
      else
        if (this%imover == 1) then
          q = this%pakmvrobj%get_qfrommvr(n)
        end if
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- STORAGE
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    do n = 1, this%maxbound
      q = dzero
      if (this%iboundpak(n) /= 0) then
        d = this%depth(n)
        a = this%calc_surface_area_wet(n, d)
        this%qauxcbc(1) = a * d
      else
        q = dzero
        this%qauxcbc(1) = dzero
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc)
    end do
    !
    ! -- MOVER
    if (this%imover == 1) then
      !
      ! -- FROM MOVER
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%maxbound)
      do n = 1, this%maxbound
        q = dzero
        if (this%iboundpak(n) /= 0) then
          q = this%pakmvrobj%get_qfrommvr(n)
        end if
        call this%budobj%budterm(idx)%update_term(n, n, q)
      end do
      !
      ! -- TO MOVER
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%maxbound)
      do n = 1, this%maxbound
        if (this%iboundpak(n) /= 0) then
          q = this%pakmvrobj%get_qtomvr(n)
          if (q > dzero) then
            q = -q
          end if
        else
          q = dzero
        end if
        call this%budobj%budterm(idx)%update_term(n, n, q)
      end do
    end if
    !
    ! -- AUXILIARY VARIABLES
    naux = this%naux
    if (naux > 0) then
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%maxbound)
      do n = 1, this%maxbound
        q = dzero
        call this%budobj%budterm(idx)%update_term(n, n, q, this%auxvar(:, n))
      end do
    end if
    !
    ! --Terms are filled, now accumulate them for this time step
    call this%budobj%accumulate_terms()
    !
    ! -- return
    return

sfr_fn()

subroutine sfrmodule::sfr_fn ( class(sfrtype)  this, real(dp), dimension(:), intent(inout)  rhs, integer(i4b), dimension(:), intent(in)  ia, integer(i4b), dimension(:), intent(in)  idxglo, class(matrixbasetype), pointer  matrix_sln  )

private

Calculate and add the Newton-Raphson terms for the SFR package to the coefficient matrix and right-hand side vector.

Parameters

	this	SfrType object
[in,out]	rhs	right-hand side vector for model
[in]	ia	solution CRS row pointers
[in]	idxglo	mapping vector for model (local) to solution (global)
	matrix_sln	solution coefficient matrix

Definition at line 2108 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model
    integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers
    integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global)
    class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: ipos
    real(DP) :: rterm
    real(DP) :: drterm
    real(DP) :: rhs1
    real(DP) :: hcof1
    real(DP) :: q1
    real(DP) :: q2
    real(DP) :: hgwf
    !
    ! -- Copy package rhs and hcof into solution rhs and amat
    do j = 1, this%nbound
      i = this%isfrorder(j)
      ! -- skip inactive reaches
      if (this%iboundpak(i) < 1) cycle
      ! -- skip if reach is not connected to gwf
      n = this%nodelist(i)
      if (n < 1) cycle
      ipos = ia(n)
      rterm = this%hcof(i) * this%xnew(n)
      ! -- calculate perturbed head
      hgwf = this%xnew(n) + dem4
      call this%sfr_solve(i, hgwf, hcof1, rhs1, update=.false.)
      q1 = rhs1 - hcof1 * hgwf
      ! -- calculate unperturbed head
      q2 = this%rhs(i) - this%hcof(i) * this%xnew(n)
      ! -- calculate derivative
      drterm = (q2 - q1) / dem4
      ! -- add terms to convert conductance formulation into
      !    newton-raphson formulation
      call matrix_sln%add_value_pos(idxglo(ipos), drterm - this%hcof(i))
      rhs(n) = rhs(n) - rterm + drterm * this%xnew(n)
    end do
    !
    ! -- return
    return

sfr_gwf_conn()

integer(i4b) function sfrmodule::sfr_gwf_conn ( class(sfrtype)  this, integer(i4b), intent(in)  n  )

private

Function to determine if a reach is connected to a gwf cell. If connected, the return value is 1. Otherwise, the return value is 0.

Returns

flag indicating if reach is connected to a gwf cell

Parameters

	this	SfrType object
[in]	n	reach number

Definition at line 3900 of file gwf-sfr.f90.

    ! -- return variable
    integer(I4B) :: sfr_gwf_conn !< flag indicating if reach is connected to a gwf cell
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    ! -- local variables
    integer(I4B) :: node
 
    sfr_gwf_conn = 0
    node = this%igwfnode(n)
    if (node > 0 .and. this%hk(n) > dzero) then
      sfr_gwf_conn = 1
    end if

sfr_obs_supported()

logical function sfrmodule::sfr_obs_supported ( class(sfrtype)  this )

private

Function to determine if observations are supported by the SFR package. Observations are supported by the SFR package.

Returns

sfr_obs_supported boolean indicating if observations are supported

Parameters

this	SfrType object

Definition at line 2810 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    !
    ! -- set boolean
    sfr_obs_supported = .true.
    !
    ! -- return
    return

sfr_options()

subroutine sfrmodule::sfr_options ( class(sfrtype), intent(inout)  this, character(len=*), intent(inout)  option, logical(lgp), intent(inout)  found  )

private

Read additional options for SFR package.

Parameters

[in,out]	this	SfrType object
[in,out]	option	option keyword string
[in,out]	found	boolean indicating if option found

Definition at line 629 of file gwf-sfr.f90.

    ! -- modules
    use openspecmodule, only: access, form
    use inputoutputmodule, only: getunit, openfile
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    character(len=*), intent(inout) :: option !< option keyword string
    logical(LGP), intent(inout) :: found !< boolean indicating if option found
    ! -- local variables
    real(DP) :: r
    character(len=MAXCHARLEN) :: fname
    character(len=MAXCHARLEN) :: keyword
    ! -- formats
    character(len=*), parameter :: fmttimeconv = &
      &"(4x, 'TIME CONVERSION VALUE (',g0,') SPECIFIED.')"
    character(len=*), parameter :: fmtlengthconv = &
      &"(4x, 'LENGTH CONVERSION VALUE (',g0,') SPECIFIED.')"
    character(len=*), parameter :: fmtpicard = &
      &"(4x, 'MAXIMUM SFR PICARD ITERATION VALUE (',i0,') SPECIFIED.')"
    character(len=*), parameter :: fmtiter = &
      &"(4x, 'MAXIMUM SFR ITERATION VALUE (',i0,') SPECIFIED.')"
    character(len=*), parameter :: fmtdmaxchg = &
      &"(4x, 'MAXIMUM DEPTH CHANGE VALUE (',g0,') SPECIFIED.')"
    character(len=*), parameter :: fmtsfrbin = &
      "(4x, 'SFR ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, /4x, &
    &'OPENED ON UNIT: ', I0)"
    !
    ! -- Check for SFR options
    found = .true.
    select case (option)
    case ('PRINT_STAGE')
      this%iprhed = 1
      write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
        ' STAGES WILL BE PRINTED TO LISTING FILE.'
    case ('STAGE')
      call this%parser%GetStringCaps(keyword)
      if (keyword == 'FILEOUT') then
        call this%parser%GetString(fname)
        this%istageout = getunit()
        call openfile(this%istageout, this%iout, fname, 'DATA(BINARY)', &
                      form, access, 'REPLACE', mnormal)
        write (this%iout, fmtsfrbin) &
          'STAGE', trim(adjustl(fname)), this%istageout
      else
        call store_error('Optional stage keyword must &
                         &be followed by fileout.')
      end if
    case ('BUDGET')
      call this%parser%GetStringCaps(keyword)
      if (keyword == 'FILEOUT') then
        call this%parser%GetString(fname)
        this%ibudgetout = getunit()
        call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
                      form, access, 'REPLACE', mnormal)
        write (this%iout, fmtsfrbin) &
          'BUDGET', trim(adjustl(fname)), this%ibudgetout
      else
        call store_error('Optional budget keyword must be '// &
                         'followed by fileout.')
      end if
    case ('BUDGETCSV')
      call this%parser%GetStringCaps(keyword)
      if (keyword == 'FILEOUT') then
        call this%parser%GetString(fname)
        this%ibudcsv = getunit()
        call openfile(this%ibudcsv, this%iout, fname, 'CSV', &
                      filstat_opt='REPLACE')
        write (this%iout, fmtsfrbin) &
          'BUDGET CSV', trim(adjustl(fname)), this%ibudcsv
      else
        call store_error('OPTIONAL BUDGETCSV KEYWORD MUST BE FOLLOWED BY &
          &FILEOUT')
      end if
    case ('PACKAGE_CONVERGENCE')
      call this%parser%GetStringCaps(keyword)
      if (keyword == 'FILEOUT') then
        call this%parser%GetString(fname)
        this%ipakcsv = getunit()
        call openfile(this%ipakcsv, this%iout, fname, 'CSV', &
                      filstat_opt='REPLACE', mode_opt=mnormal)
        write (this%iout, fmtsfrbin) &
          'PACKAGE_CONVERGENCE', trim(adjustl(fname)), this%ipakcsv
      else
        call store_error('Optional package_convergence keyword must be '// &
                         'followed by fileout.')
      end if
    case ('UNIT_CONVERSION')
      this%unitconv = this%parser%GetDouble()
      !
      ! -- create warning message
      write (warnmsg, '(a)') &
        'SETTING UNIT_CONVERSION DIRECTLY'
      !
      ! -- create deprecation warning
      call deprecation_warning('OPTIONS', 'UNIT_CONVERSION', '6.4.2', &
                               warnmsg, this%parser%GetUnit())
    case ('LENGTH_CONVERSION')
      this%lengthconv = this%parser%GetDouble()
      write (this%iout, fmtlengthconv) this%lengthconv
    case ('TIME_CONVERSION')
      this%timeconv = this%parser%GetDouble()
      write (this%iout, fmttimeconv) this%timeconv
    case ('MAXIMUM_PICARD_ITERATIONS')
      this%maxsfrpicard = this%parser%GetInteger()
      write (this%iout, fmtpicard) this%maxsfrpicard
    case ('MAXIMUM_ITERATIONS')
      this%maxsfrit = this%parser%GetInteger()
      write (this%iout, fmtiter) this%maxsfrit
    case ('MAXIMUM_DEPTH_CHANGE')
      r = this%parser%GetDouble()
      this%dmaxchg = r
      this%deps = dp999 * r
      write (this%iout, fmtdmaxchg) this%dmaxchg
    case ('MOVER')
      this%imover = 1
      write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
      !
      ! -- right now these are options that are only available in the
      !    development version and are not included in the documentation.
      !    These options are only available when IDEVELOPMODE in
      !    constants module is set to 1
    case ('DEV_NO_CHECK')
      call this%parser%DevOpt()
      this%icheck = 0
      write (this%iout, '(4x,A)') 'SFR CHECKS OF REACH GEOMETRY '// &
        'RELATIVE TO MODEL GRID AND '// &
        'REASONABLE PARAMETERS WILL NOT '// &
        'BE PERFORMED.'
    case ('DEV_NO_FINAL_CHECK')
      call this%parser%DevOpt()
      this%iconvchk = 0
      write (this%iout, '(4x,a)') &
        'A FINAL CONVERGENCE CHECK OF THE CHANGE IN STREAM FLOW ROUTING &
        &STAGES AND FLOWS WILL NOT BE MADE'
      !
      ! -- no valid options found
    case default
      !
      ! -- No options found
      found = .false.
    end select
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_ot_bdsummary()

subroutine sfrmodule::sfr_ot_bdsummary	(	class(sfrtype)	this,
		integer(i4b), intent(in)	kstp,
		integer(i4b), intent(in)	kper,
		integer(i4b), intent(in)	iout,
		integer(i4b), intent(in)	ibudfl
	)

Output SFR package budget summary.

Parameters

	this	SfrType object
[in]	kstp	time step number
[in]	kper	period number
[in]	iout	flag and unit number for the model listing file
[in]	ibudfl	flag indicating budget should be written

Definition at line 2621 of file gwf-sfr.f90.

    ! -- module
    use tdismodule, only: totim, delt
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: kstp !< time step number
    integer(I4B), intent(in) :: kper !< period number
    integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
    integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
    !
    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
    !
    ! -- return
    return

sfr_ot_dv()

subroutine sfrmodule::sfr_ot_dv	(	class(sfrtype)	this,
		integer(i4b), intent(in)	idvsave,
		integer(i4b), intent(in)	idvprint
	)

Output SFR boundary package dependent-variable terms.

Parameters

	this	SfrType object
[in]	idvsave	flag and unit number for dependent-variable output
[in]	idvprint	flag indicating if dependent-variable should be written to the model listing file

Definition at line 2505 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: kstp, kper, pertim, totim
    use inputoutputmodule, only: ulasav
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output
    integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file
    ! -- local variables
    character(len=20) :: cellid
    integer(I4B) :: ibinun
    integer(I4B) :: n
    integer(I4B) :: node
    real(DP) :: d
    real(DP) :: v
    real(DP) :: hgwf
    real(DP) :: sbot
    real(DP) :: depth
    real(DP) :: stage
    real(DP) :: w
    real(DP) :: cond
    real(DP) :: grad
    !
    ! -- set unit number for binary dependent variable output
    ibinun = 0
    if (this%istageout /= 0) then
      ibinun = this%istageout
    end if
    if (idvsave == 0) ibinun = 0
    !
    ! -- write sfr binary output
    if (ibinun > 0) then
      do n = 1, this%maxbound
        d = this%depth(n)
        v = this%stage(n)
        if (this%iboundpak(n) == 0) then
          v = dhnoflo
        else if (d == dzero) then
          v = dhdry
        end if
        this%dbuff(n) = v
      end do
      call ulasav(this%dbuff, '           STAGE', kstp, kper, pertim, totim, &
                  this%maxbound, 1, 1, ibinun)
    end if
    !
    ! -- print sfr stage and depth table
    if (idvprint /= 0 .and. this%iprhed /= 0) then
      !
      ! -- set table kstp and kper
      call this%stagetab%set_kstpkper(kstp, kper)
      !
      ! -- fill stage data
      do n = 1, this%maxbound
        node = this%igwfnode(n)
        if (node > 0) then
          call this%dis%noder_to_string(node, cellid)
          hgwf = this%xnew(node)
        else
          cellid = 'NONE'
        end if
        if (this%inamedbound == 1) then
          call this%stagetab%add_term(this%boundname(n))
        end if
        call this%stagetab%add_term(n)
        call this%stagetab%add_term(cellid)
        if (this%iboundpak(n) /= 0) then
          depth = this%depth(n)
          stage = this%stage(n)
          w = this%calc_top_width_wet(n, depth)
          call this%sfr_calc_cond(n, depth, cond, stage, hgwf)
        else
          depth = dhnoflo
          stage = dhnoflo
          w = dhnoflo
          cond = dhnoflo
        end if
        if (depth == dzero) then
          call this%stagetab%add_term(dhdry)
        else
          call this%stagetab%add_term(stage)
        end if
        call this%stagetab%add_term(depth)
        call this%stagetab%add_term(w)
        if (node > 0) then
          if (this%iboundpak(n) /= 0) then
            sbot = this%strtop(n) - this%bthick(n)
            if (hgwf < sbot) then
              grad = stage - sbot
            else
              grad = stage - hgwf
            end if
            grad = grad / this%bthick(n)
          else
            grad = dhnoflo
          end if
          call this%stagetab%add_term(hgwf)
          call this%stagetab%add_term(cond)
          call this%stagetab%add_term(grad)
        else
          call this%stagetab%add_term('--')
          call this%stagetab%add_term('--')
          call this%stagetab%add_term('--')
        end if
      end do
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_ot_package_flows()

subroutine sfrmodule::sfr_ot_package_flows	(	class(sfrtype)	this,
		integer(i4b), intent(in)	icbcfl,
		integer(i4b), intent(in)	ibudfl
	)

Output SFR package flow terms.

Parameters

	this	SfrType object
[in]	icbcfl	flag and unit number for cell-by-cell output
[in]	ibudfl	flag indication if cell-by-cell data should be saved

Definition at line 2448 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: kstp, kper, delt, pertim, totim
    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: icbcfl !< flag and unit number for cell-by-cell output
    integer(I4B), intent(in) :: ibudfl !< flag indication if cell-by-cell data should be saved
    ! -- local variables
    integer(I4B) :: ibinun
    character(len=20), dimension(:), allocatable :: cellidstr
    integer(I4B) :: n
    integer(I4B) :: node
    !
    ! -- write the flows from the budobj
    ibinun = 0
    if (this%ibudgetout /= 0) then
      ibinun = this%ibudgetout
    end if
    if (icbcfl == 0) ibinun = 0
    if (ibinun > 0) then
      call this%budobj%save_flows(this%dis, ibinun, kstp, kper, delt, &
                                  pertim, totim, this%iout)
    end if
    !
    ! -- Print sfr flows table
    if (ibudfl /= 0 .and. this%iprflow /= 0) then
      !
      ! -- If there are any 'none' gwf connections then need to calculate
      !    a vector of cellids and pass that in to the budget flow table because
      !    the table assumes that there are maxbound gwf entries, which is not
      !    the case if any 'none's are specified.
      if (this%ianynone > 0) then
        allocate (cellidstr(this%maxbound))
        do n = 1, this%maxbound
          node = this%igwfnode(n)
          if (node > 0) then
            call this%dis%noder_to_string(node, cellidstr(n))
          else
            cellidstr(n) = 'NONE'
          end if
        end do
        call this%budobj%write_flowtable(this%dis, kstp, kper, cellidstr)
        deallocate (cellidstr)
      else
        call this%budobj%write_flowtable(this%dis, kstp, kper)
      end if
    end if
    !
    ! -- return
    return

sfr_process_obsid()

subroutine sfrmodule::sfr_process_obsid	(	type(observetype), intent(inout)	obsrv,
		class(disbasetype), intent(in)	dis,
		integer(i4b), intent(in)	inunitobs,
		integer(i4b), intent(in)	iout
	)

Method to process observation ID strings for a SFR package.

Parameters

[in,out]	obsrv	Observation object
[in]	dis	Discretization object
[in]	inunitobs	file unit number for the package observation file
[in]	iout	model listing file unit number

Definition at line 3136 of file gwf-sfr.f90.

    ! -- dummy variables
    type(ObserveType), intent(inout) :: obsrv !< Observation object
    class(DisBaseType), intent(in) :: dis !< Discretization object
    integer(I4B), intent(in) :: inunitobs !< file unit number for the package observation file
    integer(I4B), intent(in) :: iout !< model listing file unit number
    ! -- local variables
    integer(I4B) :: nn1
    integer(I4B) :: icol
    integer(I4B) :: istart
    integer(I4B) :: istop
    character(len=LINELENGTH) :: string
    character(len=LENBOUNDNAME) :: bndname
    !
    ! -- initialize local variables
    string = obsrv%IDstring
    !
    ! -- Extract reach number from string and store it.
    !    If 1st item is not an integer(I4B), it should be a
    !    boundary name--deal with it.
    icol = 1
    !
    ! -- get reach number or boundary name
    call extract_idnum_or_bndname(string, icol, istart, istop, nn1, bndname)
    if (nn1 == namedboundflag) then
      obsrv%FeatureName = bndname
    end if
    !
    ! -- store reach number (NodeNumber)
    obsrv%NodeNumber = nn1
    !
    ! -- return
    return

Here is the call graph for this function:

Here is the caller graph for this function:

sfr_read_connectiondata()

subroutine sfrmodule::sfr_read_connectiondata

(

class(sfrtype), intent(inout)

this

)

Method to read connectiondata for each reach for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 1239 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    use sparsemodule, only: sparsematrix
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    character(len=LINELENGTH) :: line
    logical(LGP) :: isfound
    logical(LGP) :: endOfBlock
    integer(I4B) :: n
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: jj
    integer(I4B) :: jcol
    integer(I4B) :: jcol2
    integer(I4B) :: nja
    integer(I4B) :: ival
    integer(I4B) :: idir
    integer(I4B) :: ierr
    integer(I4B) :: nconnmax
    integer(I4B) :: nup
    integer(I4B) :: ipos
    integer(I4B) :: istat
    integer(I4B), dimension(:), pointer, contiguous :: rowmaxnnz => null()
    integer, allocatable, dimension(:) :: nboundchk
    integer, allocatable, dimension(:, :) :: iconndata
    type(sparsematrix), pointer :: sparse => null()
    integer(I4B), dimension(:), allocatable :: iup
    integer(I4B), dimension(:), allocatable :: order
    type(dag) :: sfr_dag
    !
    ! -- allocate and initialize local variables for reach connections
    allocate (nboundchk(this%maxbound))
    do n = 1, this%maxbound
      nboundchk(n) = 0
    end do
    !
    ! -- calculate the number of non-zero entries (size of ja maxtrix)
    nja = 0
    nconnmax = 0
    allocate (rowmaxnnz(this%maxbound))
    do n = 1, this%maxbound
      ival = this%nconnreach(n)
      if (ival < 0) ival = 0
      rowmaxnnz(n) = ival + 1
      nja = nja + ival + 1
      if (ival > nconnmax) then
        nconnmax = ival
      end if
    end do
    !
    ! -- reallocate connection data for package
    call mem_reallocate(this%ja, nja, 'JA', this%memoryPath)
    call mem_reallocate(this%idir, nja, 'IDIR', this%memoryPath)
    call mem_reallocate(this%idiv, nja, 'IDIV', this%memoryPath)
    call mem_reallocate(this%qconn, nja, 'QCONN', this%memoryPath)
    !
    ! -- initialize connection data
    do n = 1, nja
      this%idir(n) = 0
      this%idiv(n) = 0
      this%qconn(n) = dzero
    end do
    !
    ! -- allocate space for iconndata
    allocate (iconndata(nconnmax, this%maxbound))
    !
    ! -- initialize iconndata
    do n = 1, this%maxbound
      do j = 1, nconnmax
        iconndata(j, n) = 0
      end do
    end do
    !
    ! -- allocate space for connectivity
    allocate (sparse)
    !
    ! -- set up sparse
    call sparse%init(this%maxbound, this%maxbound, rowmaxnnz)
    !
    ! -- read connection data
    call this%parser%GetBlock('CONNECTIONDATA', isfound, ierr, &
                              supportopenclose=.true.)
    !
    ! -- parse reach connectivity block if detected
    if (isfound) then
      write (this%iout, '(/1x,a)') &
        'PROCESSING '//trim(adjustl(this%text))//' CONNECTIONDATA'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        !
        ! -- get reach number
        n = this%parser%GetInteger()
        !
        ! -- check for error
        if (n < 1 .or. n > this%maxbound) then
          write (errmsg, '(a,1x,a,1x,i0)') &
            'SFR reach in connectiondata block is less than one or greater', &
            'than NREACHES:', n
          call store_error(errmsg)
          cycle
        end if
        !
        ! -- increment nboundchk
        nboundchk(n) = nboundchk(n) + 1
        !
        ! -- add diagonal connection for reach
        call sparse%addconnection(n, n, 1)
        !
        ! -- fill off diagonals
        do i = 1, this%nconnreach(n)
          !
          ! -- get connected reach
          ival = this%parser%GetInteger()
          !
          ! -- save connection data to temporary iconndata
          iconndata(i, n) = ival
          !
          ! -- determine idir
          if (ival < 0) then
            idir = -1
            ival = abs(ival)
          elseif (ival == 0) then
            call store_error('Missing or zero connection reach in line:')
            call store_error(line)
          else
            idir = 1
          end if
          if (ival > this%maxbound) then
            call store_error('Reach number exceeds NREACHES in line:')
            call store_error(line)
          end if
          !
          ! -- add connection to sparse
          call sparse%addconnection(n, ival, 1)
        end do
      end do
 
      write (this%iout, '(1x,a)') &
        'END OF '//trim(adjustl(this%text))//' CONNECTIONDATA'
 
      do n = 1, this%maxbound
        !
        ! -- check for missing or duplicate sfr connections
        if (nboundchk(n) == 0) then
          write (errmsg, '(a,1x,i0)') &
            'No connection data specified for reach', n
          call store_error(errmsg)
        else if (nboundchk(n) > 1) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
            'Connection data for reach', n, &
            'specified', nboundchk(n), 'times.'
          call store_error(errmsg)
        end if
      end do
    else
      call store_error('Required connectiondata block not found.')
    end if
    !
    ! -- terminate if errors encountered in connectiondata block
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- create ia and ja from sparse
    call sparse%filliaja(this%ia, this%ja, ierr, sort=.true.)
    !
    ! -- test for error condition
    if (ierr /= 0) then
      write (errmsg, '(a,3(1x,a))') &
        'Could not fill', trim(this%packName), &
        'package IA and JA connection data.', &
        'Check connectivity data in connectiondata block.'
      call store_error(errmsg)
    end if
    !
    ! -- fill flat connection storage
    do n = 1, this%maxbound
      do j = this%ia(n) + 1, this%ia(n + 1) - 1
        jcol = this%ja(j)
        do jj = 1, this%nconnreach(n)
          jcol2 = iconndata(jj, n)
          if (abs(jcol2) == jcol) then
            idir = 1
            if (jcol2 < 0) then
              idir = -1
            end if
            this%idir(j) = idir
            exit
          end if
        end do
      end do
    end do
    !
    ! -- deallocate temporary local storage for reach connections
    deallocate (rowmaxnnz)
    deallocate (nboundchk)
    deallocate (iconndata)
    !
    ! -- destroy sparse
    call sparse%destroy()
    deallocate (sparse)
    !
    ! -- calculate reach order using DAG
    !
    ! -- initialize the DAG
    call sfr_dag%set_vertices(this%maxbound)
    !
    ! -- fill DAG
    fill_dag: do n = 1, this%maxbound
      !
      ! -- determine the number of upstream reaches
      nup = 0
      do j = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(j) > 0) then
          nup = nup + 1
        end if
      end do
      !
      ! -- cycle if nu upstream reacches
      if (nup == 0) cycle fill_dag
      !
      ! -- allocate local storage
      allocate (iup(nup))
      !
      ! -- fill local storage
      ipos = 1
      do j = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(j) > 0) then
          iup(ipos) = this%ja(j)
          ipos = ipos + 1
        end if
      end do
      !
      ! -- add upstream connections to DAG
      call sfr_dag%set_edges(n, iup)
      !
      ! -- clean up local storage
      deallocate (iup)
    end do fill_dag
    !
    ! -- perform toposort on DAG
    call sfr_dag%toposort(order, istat)
    !
    ! -- write warning if circular dependency
    if (istat == -1) then
      write (warnmsg, '(a)') &
        trim(adjustl(this%text))//' PACKAGE ('// &
        trim(adjustl(this%packName))//') cannot calculate a '// &
        'Directed Asyclic Graph for reach connectivity because '// &
        'of circular dependency. Using the reach number for '// &
        'solution ordering.'
      call store_warning(warnmsg)
    end if
    !
    ! -- fill isfrorder
    do n = 1, this%maxbound
      if (istat == 0) then
        this%isfrorder(n) = order(n)
      else
        this%isfrorder(n) = n
      end if
    end do
    !
    ! -- clean up DAG and remaining local storage
    call sfr_dag%destroy()
    if (istat == 0) then
      deallocate (order)
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_read_crossection()

subroutine sfrmodule::sfr_read_crossection

(

class(sfrtype), intent(inout)

this

)

Method to read crosssection data for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 1097 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    use sfrcrosssectionmanager, only: cross_section_cr, sfrcrosssection
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    character(len=LINELENGTH) :: keyword
    character(len=LINELENGTH) :: line
    logical(LGP) :: isfound
    logical(LGP) :: endOfBlock
    integer(I4B) :: n
    integer(I4B) :: ierr
    integer(I4B) :: ncrossptstot
    integer, allocatable, dimension(:) :: nboundchk
    type(SfrCrossSection), pointer :: cross_data => null()
    !
    ! -- read cross-section data
    call this%parser%GetBlock('CROSSSECTIONS', isfound, ierr, &
                              supportopenclose=.true., &
                              blockrequired=.false.)
    !
    ! -- parse reach connectivity block if detected
    if (isfound) then
      write (this%iout, '(/1x,a)') &
        'PROCESSING '//trim(adjustl(this%text))//' CROSSSECTIONS'
      !
      ! -- allocate and initialize local variables for reach cross-sections
      allocate (nboundchk(this%maxbound))
      do n = 1, this%maxbound
        nboundchk(n) = 0
      end do
      !
      ! -- create and initialize cross-section data
      call cross_section_cr(cross_data, this%iout, this%iprpak, this%maxbound)
      call cross_data%initialize(this%ncrossptstot, this%ncrosspts, &
                                 this%iacross, &
                                 this%station, this%xsheight, &
                                 this%xsrough)
      !
      ! -- read all of the entries in the block
      readtable: do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        !
        ! -- get reach number
        n = this%parser%GetInteger()
        !
        ! -- check for reach number error
        if (n < 1 .or. n > this%maxbound) then
          write (errmsg, '(a,1x,a,1x,i0)') &
            'SFR reach in crosssections block is less than one or greater', &
            'than NREACHES:', n
          call store_error(errmsg)
          cycle readtable
        end if
        !
        ! -- increment nboundchk
        nboundchk(n) = nboundchk(n) + 1
        !
        ! -- read FILE keyword
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('TAB6')
          call this%parser%GetStringCaps(keyword)
          if (trim(adjustl(keyword)) /= 'FILEIN') then
            errmsg = 'TAB6 keyword must be followed by "FILEIN" '// &
                     'then by filename.'
            call store_error(errmsg)
            cycle readtable
          end if
          call this%parser%GetString(line)
          call cross_data%read_table(n, this%width(n), &
                                     trim(adjustl(line)))
        case default
          write (errmsg, '(a,1x,i4,1x,a)') &
            'CROSS-SECTION TABLE ENTRY for REACH ', n, &
            'MUST INCLUDE TAB6 KEYWORD'
          call store_error(errmsg)
          cycle readtable
        end select
      end do readtable
 
      write (this%iout, '(1x,a)') &
        'END OF '//trim(adjustl(this%text))//' CROSSSECTIONS'
 
      !
      ! -- check for duplicate sfr crosssections
      do n = 1, this%maxbound
        if (nboundchk(n) > 1) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
            'Cross-section data for reach', n, &
            'specified', nboundchk(n), 'times.'
          call store_error(errmsg)
        end if
      end do
      !
      ! -- terminate if errors encountered in cross-sections block
      if (count_errors() > 0) then
        call this%parser%StoreErrorUnit()
      end if
      !
      ! -- determine the current size of cross-section data
      ncrossptstot = cross_data%get_ncrossptstot()
      !
      ! -- reallocate sfr package cross-section data
      if (ncrossptstot /= this%ncrossptstot) then
        this%ncrossptstot = ncrossptstot
        call mem_reallocate(this%station, this%ncrossptstot, 'STATION', &
                            this%memoryPath)
        call mem_reallocate(this%xsheight, this%ncrossptstot, 'XSHEIGHT', &
                            this%memoryPath)
        call mem_reallocate(this%xsrough, this%ncrossptstot, 'XSROUGH', &
                            this%memoryPath)
      end if
      !
      ! -- write cross-section data to the model listing file
      call cross_data%output(this%width, this%rough)
      !
      ! -- pack cross-section data
      call cross_data%pack(this%ncrossptstot, this%ncrosspts, &
                           this%iacross, &
                           this%station, &
                           this%xsheight, &
                           this%xsrough)
      !
      ! -- deallocate temporary local storage for reach cross-sections
      deallocate (nboundchk)
      call cross_data%destroy()
      deallocate (cross_data)
      nullify (cross_data)
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_read_dimensions()

subroutine sfrmodule::sfr_read_dimensions

(

class(sfrtype), intent(inout)

this

)

Read dimensions for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 540 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    character(len=LINELENGTH) :: keyword
    integer(I4B) :: ierr
    logical(LGP) :: isfound
    logical(LGP) :: endOfBlock
    !
    ! -- initialize dimensions to 0
    this%maxbound = 0
    !
    ! -- get dimensions block
    call this%parser%GetBlock('DIMENSIONS', isfound, ierr, &
                              supportopenclose=.true.)
    !
    ! -- parse dimensions block if detected
    if (isfound) then
      write (this%iout, '(/1x,a)') &
        'PROCESSING '//trim(adjustl(this%text))//' DIMENSIONS'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('NREACHES')
          this%maxbound = this%parser%GetInteger()
          write (this%iout, '(4x,a,i0)') 'NREACHES = ', this%maxbound
        case default
          write (errmsg, '(2a)') &
            'Unknown '//trim(this%text)//' dimension: ', trim(keyword)
          call store_error(errmsg)
        end select
      end do
      write (this%iout, '(1x,a)') &
        'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
    else
      call store_error('Required dimensions block not found.')
    end if
    !
    ! -- verify dimensions were set
    if (this%maxbound < 1) then
      write (errmsg, '(a)') &
        'NREACHES was not specified or was specified incorrectly.'
      call store_error(errmsg)
    end if
    !
    ! -- write summary of error messages for block
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- Call define_listlabel to construct the list label that is written
    !    when PRINT_INPUT option is used.
    call this%define_listlabel()
    !
    ! -- Define default cross-section data size
    this%ncrossptstot = this%maxbound
    !
    ! -- Allocate arrays in package superclass
    call this%sfr_allocate_arrays()
    !
    ! -- read package data
    call this%sfr_read_packagedata()
    !
    ! -- read cross-section data
    call this%sfr_read_crossection()
    !
    ! -- read connection data
    call this%sfr_read_connectiondata()
    !
    ! -- read diversion data
    call this%sfr_read_diversions()
    !
    ! -- setup the budget object
    call this%sfr_setup_budobj()
    !
    ! -- setup the stage table object
    call this%sfr_setup_tableobj()
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_read_diversions()

subroutine sfrmodule::sfr_read_diversions

(

class(sfrtype), intent(inout)

this

)

Method to read diversions for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 1520 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    character(len=10) :: cnum
    character(len=10) :: cval
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: ierr
    integer(I4B) :: ival
    integer(I4B) :: i0
    integer(I4B) :: ipos
    integer(I4B) :: jpos
    integer(I4B) :: ndiv
    integer(I4B) :: ndiversions
    integer(I4B) :: idivreach
    logical(LGP) :: isfound
    logical(LGP) :: endOfBlock
    integer(I4B) :: idiv
    integer, allocatable, dimension(:) :: iachk
    integer, allocatable, dimension(:) :: nboundchk
    !
    ! -- determine the total number of diversions and fill iadiv
    ndiversions = 0
    i0 = 1
    this%iadiv(1) = i0
    do n = 1, this%maxbound
      ndiversions = ndiversions + this%ndiv(n)
      i0 = i0 + this%ndiv(n)
      this%iadiv(n + 1) = i0
    end do
    !
    ! -- reallocate memory for diversions
    if (ndiversions > 0) then
      call mem_reallocate(this%divreach, ndiversions, 'DIVREACH', &
                          this%memoryPath)
      allocate (this%divcprior(ndiversions))
      call mem_reallocate(this%divflow, ndiversions, 'DIVFLOW', this%memoryPath)
      call mem_reallocate(this%divq, ndiversions, 'DIVQ', this%memoryPath)
    end if
    !
    ! -- initialize diversion flow
    do n = 1, ndiversions
      this%divflow(n) = dzero
      this%divq(n) = dzero
    end do
    !
    ! -- read diversions
    call this%parser%GetBlock('DIVERSIONS', isfound, ierr, &
                              supportopenclose=.true., &
                              blockrequired=.false.)
    !
    ! -- parse reach connectivity block if detected
    if (isfound) then
      if (this%idiversions /= 0) then
        write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
          ' DIVERSIONS'
        !
        ! -- allocate and initialize local variables for diversions
        ndiv = 0
        do n = 1, this%maxbound
          ndiv = ndiv + this%ndiv(n)
        end do
        allocate (iachk(this%maxbound + 1))
        allocate (nboundchk(ndiv))
        iachk(1) = 1
        do n = 1, this%maxbound
          iachk(n + 1) = iachk(n) + this%ndiv(n)
        end do
        do n = 1, ndiv
          nboundchk(n) = 0
        end do
        !
        ! -- read diversion data
        do
          call this%parser%GetNextLine(endofblock)
          if (endofblock) exit
          !
          ! -- get reach number
          n = this%parser%GetInteger()
          if (n < 1 .or. n > this%maxbound) then
            write (cnum, '(i0)') n
            errmsg = 'Reach number should be between 1 and '// &
                     trim(cnum)//'.'
            call store_error(errmsg)
            cycle
          end if
          !
          ! -- make sure reach has at least one diversion
          if (this%ndiv(n) < 1) then
            write (cnum, '(i0)') n
            errmsg = 'Diversions cannot be specified '// &
                     'for reach '//trim(cnum)
            call store_error(errmsg)
            cycle
          end if
          !
          ! -- read diversion number
          ival = this%parser%GetInteger()
          if (ival < 1 .or. ival > this%ndiv(n)) then
            write (cnum, '(i0)') n
            errmsg = 'Reach  '//trim(cnum)
            write (cnum, '(i0)') this%ndiv(n)
            errmsg = trim(errmsg)//' diversion number should be between '// &
                     '1 and '//trim(cnum)//'.'
            call store_error(errmsg)
            cycle
          end if
 
          ! -- increment nboundchk
          ipos = iachk(n) + ival - 1
          nboundchk(ipos) = nboundchk(ipos) + 1
 
          idiv = ival
          !
          ! -- get target reach for diversion
          ival = this%parser%GetInteger()
          if (ival < 1 .or. ival > this%maxbound) then
            write (cnum, '(i0)') ival
            errmsg = 'Diversion target reach number should be '// &
                     'between 1 and '//trim(cnum)//'.'
            call store_error(errmsg)
            cycle
          end if
          idivreach = ival
          jpos = this%iadiv(n) + idiv - 1
          this%divreach(jpos) = idivreach
          !
          ! -- get cprior
          call this%parser%GetStringCaps(cval)
          ival = -1
          select case (cval)
          case ('UPTO')
            ival = 0
          case ('THRESHOLD')
            ival = -1
          case ('FRACTION')
            ival = -2
          case ('EXCESS')
            ival = -3
          case default
            errmsg = 'Invalid cprior type '//trim(cval)//'.'
            call store_error(errmsg)
          end select
          !
          ! -- set cprior for diversion
          this%divcprior(jpos) = cval
        end do
 
        write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))// &
          ' DIVERSIONS'
 
        do n = 1, this%maxbound
          do j = 1, this%ndiv(n)
            ipos = iachk(n) + j - 1
            !
            ! -- check for missing or duplicate reach diversions
            if (nboundchk(ipos) == 0) then
              write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
                'No data specified for reach', n, 'diversion', j
              call store_error(errmsg)
            else if (nboundchk(ipos) > 1) then
              write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
                'Data for reach', n, 'diversion', j, &
                'specified', nboundchk(ipos), 'times'
              call store_error(errmsg)
            end if
          end do
        end do
        !
        ! -- deallocate local variables
        deallocate (iachk)
        deallocate (nboundchk)
      else
        !
        ! -- error condition
        write (errmsg, '(a,1x,a)') &
          'A diversions block should not be', &
          'specified if diversions are not specified.'
        call store_error(errmsg)
      end if
    else
      if (this%idiversions /= 0) then
        call store_error('REQUIRED DIVERSIONS BLOCK NOT FOUND.')
      end if
    end if
    !
    ! -- write summary of diversion error messages
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_read_packagedata()

subroutine sfrmodule::sfr_read_packagedata ( class(sfrtype), intent(inout)  this )

private

Method to read packagedata for each reach for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 843 of file gwf-sfr.f90.

    ! -- modules
    use timeseriesmanagermodule, only: read_value_or_time_series_adv
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    character(len=LINELENGTH) :: text
    character(len=LINELENGTH) :: cellid
    character(len=10) :: cnum
    character(len=LENBOUNDNAME) :: bndName
    character(len=LENBOUNDNAME) :: bndNameTemp
    character(len=LENBOUNDNAME) :: hkname
    character(len=LENBOUNDNAME) :: manningname
    character(len=LENBOUNDNAME) :: ustrfname
    character(len=50), dimension(:), allocatable :: caux
    integer(I4B) :: n, ierr, ival
    logical(LGP) :: isfound
    logical(LGP) :: endOfBlock
    integer(I4B) :: i
    integer(I4B) :: ii
    integer(I4B) :: jj
    integer(I4B) :: iaux
    integer(I4B) :: nconzero
    integer(I4B) :: ipos
    integer, allocatable, dimension(:) :: nboundchk
    real(DP), pointer :: bndElem => null()
    !
    ! -- allocate space for checking sfr reach data
    allocate (nboundchk(this%maxbound))
    do i = 1, this%maxbound
      nboundchk(i) = 0
    end do
    nconzero = 0
    !
    ! -- allocate local storage for aux variables
    if (this%naux > 0) then
      allocate (caux(this%naux))
    end if
    !
    ! -- read reach data
    call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
                              supportopenclose=.true.)
    !
    ! -- parse reaches block if detected
    if (isfound) then
      write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
        ' PACKAGEDATA'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        ! -- read reach number
        n = this%parser%GetInteger()
 
        if (n < 1 .or. n > this%maxbound) then
          write (errmsg, '(a,1x,a,1x,i0)') &
            'Reach number (rno) must be greater than 0 and less', &
            'than or equal to', this%maxbound
          call store_error(errmsg)
          cycle
        end if
 
        ! -- increment nboundchk
        nboundchk(n) = nboundchk(n) + 1
        !
        ! -- get model node number
        call this%parser%GetCellid(this%dis%ndim, cellid, flag_string=.true.)
        this%igwfnode(n) = this%dis%noder_from_cellid(cellid, this%inunit, &
                                                      this%iout, &
                                                      flag_string=.true., &
                                                      allow_zero=.true.)
        this%igwftopnode(n) = this%igwfnode(n)
        !
        ! -- read the cellid string and determine if 'none' is specified
        if (this%igwfnode(n) < 1) then
          this%ianynone = this%ianynone + 1
          call upcase(cellid)
          if (cellid == 'NONE') then
            call this%parser%GetStringCaps(cellid)
            !
            ! -- create warning message
            write (cnum, '(i0)') n
            warnmsg = 'CELLID for unconnected reach '//trim(cnum)// &
                      ' specified to be NONE. Unconnected reaches '// &
                      'should be specified with a zero for each grid '// &
                      'dimension. For example, for a DIS grid a CELLID '// &
                      'of 0 0 0 should be specified for unconnected reaches'
            !
            ! -- create deprecation warning
            call deprecation_warning('PACKAGEDATA', 'CELLID=NONE', '6.4.3', &
                                     warnmsg, this%parser%GetUnit())
          else
 
          end if
        end if
        ! -- get reach length
        this%length(n) = this%parser%GetDouble()
        ! -- get reach width
        this%width(n) = this%parser%GetDouble()
        ! -- get reach slope
        this%slope(n) = this%parser%GetDouble()
        ! -- get reach stream bottom
        this%strtop(n) = this%parser%GetDouble()
        ! -- get reach bed thickness
        this%bthick(n) = this%parser%GetDouble()
        ! -- get reach bed hk
        call this%parser%GetStringCaps(hkname)
        ! -- get reach roughness
        call this%parser%GetStringCaps(manningname)
        ! -- get number of connections for reach
        ival = this%parser%GetInteger()
        this%nconnreach(n) = ival
        this%nconn = this%nconn + ival
        if (ival < 0) then
          write (errmsg, '(a,1x,i0,1x,a,i0,a)') &
            'NCON for reach', n, &
            'must be greater than or equal to 0 (', ival, ').'
          call store_error(errmsg)
        else if (ival == 0) then
          nconzero = nconzero + 1
        end if
        ! -- get upstream fraction for reach
        call this%parser%GetString(ustrfname)
        ! -- get number of diversions for reach
        ival = this%parser%GetInteger()
        this%ndiv(n) = ival
        if (ival > 0) then
          this%idiversions = 1
        else if (ival < 0) then
          ival = 0
        end if
 
        ! -- get aux data
        do iaux = 1, this%naux
          call this%parser%GetString(caux(iaux))
        end do
 
        ! -- set default bndName
        write (cnum, '(i10.10)') n
        bndname = 'Reach'//cnum
 
        ! -- get reach name
        if (this%inamedbound /= 0) then
          call this%parser%GetStringCaps(bndnametemp)
          if (bndnametemp /= '') then
            bndname = bndnametemp
          end if
          !this%boundname(n) = bndName
        end if
        this%sfrname(n) = bndname
        !
        ! -- set reach hydraulic conductivity
        text = hkname
        jj = 1 !for 'BEDK'
        bndelem => this%hk(n)
        call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                           this%packName, 'BND', &
                                           this%tsManager, this%iprpak, &
                                           'BEDK')
        !
        ! -- set Mannings
        text = manningname
        jj = 1 !for 'MANNING'
        bndelem => this%rough(n)
        call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                           this%packName, 'BND', &
                                           this%tsManager, this%iprpak, &
                                           'MANNING')
        !
        ! -- set upstream fraction
        text = ustrfname
        jj = 1 ! For 'USTRF'
        bndelem => this%ustrf(n)
        call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                           this%packName, 'BND', &
                                           this%tsManager, this%iprpak, 'USTRF')
        !
        ! -- get aux data
        do jj = 1, this%naux
          text = caux(jj)
          ii = n
          bndelem => this%rauxvar(jj, ii)
          call read_value_or_time_series_adv(text, ii, jj, bndelem, &
                                             this%packName, 'AUX', &
                                             this%tsManager, this%iprpak, &
                                             this%auxname(jj))
        end do
        !
        ! -- initialize sstage to the top of the reach
        !    this value would be used by simple routing reaches
        !    on kper = 1 and kstp = 1 if a stage is not specified
        !    on the status line for the reach
        this%sstage(n) = this%strtop(n)
 
      end do
      write (this%iout, '(1x,a)') &
        'END OF '//trim(adjustl(this%text))//' PACKAGEDATA'
    else
      call store_error('REQUIRED PACKAGEDATA BLOCK NOT FOUND.')
    end if
    !
    ! -- Check to make sure that every reach is specified and that no reach
    !    is specified more than once.
    do i = 1, this%maxbound
      if (nboundchk(i) == 0) then
        write (errmsg, '(a,i0,1x,a)') &
          'Information for reach ', i, 'not specified in packagedata block.'
        call store_error(errmsg)
      else if (nboundchk(i) > 1) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
          'Reach information specified', nboundchk(i), 'times for reach', i
        call store_error(errmsg)
      end if
    end do
    deallocate (nboundchk)
    !
    ! -- Submit warning message if any reach has zero connections
    if (nconzero > 0) then
      write (warnmsg, '(a,1x,a,1x,a,1x,i0,1x, a)') &
        'SFR Package', trim(this%packName), &
        'has', nconzero, 'reach(es) with zero connections.'
      call store_warning(warnmsg)
    end if
    !
    ! -- terminate if errors encountered in reach block
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- initialize the cross-section data
    ipos = 1
    this%iacross(1) = ipos
    do i = 1, this%maxbound
      this%ncrosspts(i) = 1
      this%station(ipos) = this%width(i)
      this%xsheight(ipos) = dzero
      this%xsrough(ipos) = done
      ipos = ipos + 1
      this%iacross(i + 1) = ipos
    end do
    !
    ! -- deallocate local storage for aux variables
    if (this%naux > 0) then
      deallocate (caux)
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_rp()

subroutine sfrmodule::sfr_rp

(

class(sfrtype), intent(inout)

this

)

Method to read and prepare period data for the SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 1723 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: kper, nper
    use memorymanagermodule, only: mem_reallocate
    use sfrcrosssectionmanager, only: cross_section_cr, sfrcrosssection
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: line
    character(len=LINELENGTH) :: crossfile
    integer(I4B) :: ierr
    integer(I4B) :: n
    integer(I4B) :: ichkustrm
    integer(I4B) :: ichkcross
    integer(I4B) :: ncrossptstot
    logical(LGP) :: isfound
    logical(LGP) :: endOfBlock
    type(SfrCrossSection), pointer :: cross_data => null()
    ! -- formats
    character(len=*), parameter :: fmtblkerr = &
      &"('Looking for BEGIN PERIOD iper.  Found ', a, ' instead.')"
    character(len=*), parameter :: fmtlsp = &
      &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
    character(len=*), parameter :: fmtnbd = &
      "(1X,/1X,'The number of active ',A,'S (',I6, &
      &') is greater than maximum (',I6,')')"
    !
    ! -- initialize flags
    ichkustrm = 0
    ichkcross = 0
    if (kper == 1) then
      ichkustrm = 1
    end if
    !
    ! -- set nbound to maxbound
    this%nbound = this%maxbound
    !
    ! -- Set ionper to the stress period number for which a new block of data
    !    will be read.
    if (this%ionper < kper) then
      !
      ! -- get period block
      call this%parser%GetBlock('PERIOD', isfound, ierr, &
                                supportopenclose=.true., &
                                blockrequired=.false.)
      if (isfound) then
        !
        ! -- read ionper and check for increasing period numbers
        call this%read_check_ionper()
      else
        !
        ! -- PERIOD block not found
        if (ierr < 0) then
          ! -- End of file found; data applies for remainder of simulation.
          this%ionper = nper + 1
        else
          ! -- Found invalid block
          call this%parser%GetCurrentLine(line)
          write (errmsg, fmtblkerr) adjustl(trim(line))
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
      end if
    end if
    !
    ! -- Read data if ionper == kper
    if (this%ionper == kper) then
      !
      ! -- create and initialize cross-section data
      call cross_section_cr(cross_data, this%iout, this%iprpak, this%maxbound)
      call cross_data%initialize(this%ncrossptstot, this%ncrosspts, &
                                 this%iacross, &
                                 this%station, this%xsheight, &
                                 this%xsrough)
      !
      ! -- setup table for period data
      if (this%iprpak /= 0) then
        !
        ! -- reset the input table object
        title = trim(adjustl(this%text))//' PACKAGE ('// &
                trim(adjustl(this%packName))//') DATA FOR PERIOD'
        write (title, '(a,1x,i6)') trim(adjustl(title)), kper
        call table_cr(this%inputtab, this%packName, title)
        call this%inputtab%table_df(1, 4, this%iout, finalize=.false.)
        text = 'NUMBER'
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
        text = 'KEYWORD'
        call this%inputtab%initialize_column(text, 20, alignment=tableft)
        do n = 1, 2
          write (text, '(a,1x,i6)') 'VALUE', n
          call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
        end do
      end if
      !
      ! -- read data
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        n = this%parser%GetInteger()
        if (n < 1 .or. n > this%maxbound) then
          write (errmsg, '(a,1x,a,1x,i0,a)') &
            'Reach number (RNO) must be greater than 0 and', &
            'less than or equal to', this%maxbound, '.'
          call store_error(errmsg)
          cycle
        end if
        !
        ! -- read data from the rest of the line
        call this%sfr_set_stressperiod(n, ichkustrm, crossfile)
        !
        ! -- write line to table
        if (this%iprpak /= 0) then
          call this%parser%GetCurrentLine(line)
          call this%inputtab%line_to_columns(line)
        end if
        !
        ! -- process cross-section file
        if (trim(adjustl(crossfile)) /= 'NONE') then
          call cross_data%read_table(n, this%width(n), &
                                     trim(adjustl(crossfile)))
        end if
      end do
      !
      ! -- write raw period data
      if (this%iprpak /= 0) then
        call this%inputtab%finalize_table()
      end if
      !
      ! -- finalize cross-sections
 
      !
      ! -- determine the current size of cross-section data
      ncrossptstot = cross_data%get_ncrossptstot()
      !
      ! -- reallocate sfr package cross-section data
      if (ncrossptstot /= this%ncrossptstot) then
        this%ncrossptstot = ncrossptstot
        call mem_reallocate(this%station, this%ncrossptstot, 'STATION', &
                            this%memoryPath)
        call mem_reallocate(this%xsheight, this%ncrossptstot, 'XSHEIGHT', &
                            this%memoryPath)
        call mem_reallocate(this%xsrough, this%ncrossptstot, 'XSROUGH', &
                            this%memoryPath)
      end if
      !
      ! -- write cross-section data to the model listing file
      call cross_data%output(this%width, this%rough, kstp=1, kper=kper)
      !
      ! -- pack cross-section data
      call cross_data%pack(this%ncrossptstot, this%ncrosspts, &
                           this%iacross, &
                           this%station, &
                           this%xsheight, &
                           this%xsrough)
      !
      ! -- deallocate temporary local storage for reach cross-sections
      call cross_data%destroy()
      deallocate (cross_data)
      nullify (cross_data)
      !
      ! -- Reuse data from last stress period
    else
      write (this%iout, fmtlsp) trim(this%filtyp)
    end if
    !
    ! -- check upstream fraction values
    if (ichkustrm /= 0) then
      call this%sfr_check_ustrf()
    end if
    !
    ! -- write summary of package block error messages
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_rp_obs()

subroutine sfrmodule::sfr_rp_obs ( class(sfrtype), intent(inout)  this )

private

Method to read and prepare observations for a SFR package.

Parameters

[in,out]

this

SfrType object

Definition at line 3022 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: kper
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: nn1
    character(len=LENBOUNDNAME) :: bname
    logical(LGP) :: jfound
    class(ObserveType), pointer :: obsrv => null()
    ! -- formats
10  format('Boundary "', a, '" for observation "', a, &
           '" is invalid in package "', a, '"')
30  format('Boundary name not provided for observation "', a, &
           '" in package "', a, '"')
    !
    ! -- process each package observation
    !    only done the first stress period since boundaries are fixed
    !    for the simulation
    if (kper == 1) then
      do i = 1, this%obs%npakobs
        obsrv => this%obs%pakobs(i)%obsrv
        !
        ! -- get node number 1
        nn1 = obsrv%NodeNumber
        if (nn1 == namedboundflag) then
          bname = obsrv%FeatureName
          if (bname /= '') then
            ! -- Observation location(s) is(are) based on a boundary name.
            !    Iterate through all boundaries to identify and store
            !    corresponding index(indices) in bound array.
            jfound = .false.
            do j = 1, this%maxbound
              if (this%boundname(j) == bname) then
                jfound = .true.
                call obsrv%AddObsIndex(j)
              end if
            end do
            if (.not. jfound) then
              write (errmsg, 10) &
                trim(bname), trim(obsrv%name), trim(this%packName)
              call store_error(errmsg)
            end if
          else
            write (errmsg, 30) trim(obsrv%name), trim(this%packName)
            call store_error(errmsg)
          end if
        else if (nn1 < 1 .or. nn1 > this%maxbound) then
          write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
            trim(adjustl(obsrv%ObsTypeId)), &
            'reach must be greater than 0 and less than or equal to', &
            this%maxbound, '(specified value is ', nn1, ')'
          call store_error(errmsg)
        else
          if (obsrv%indxbnds_count == 0) then
            call obsrv%AddObsIndex(nn1)
          else
            errmsg = 'Programming error in sfr_rp_obs'
            call store_error(errmsg)
          end if
        end if
        !
        ! -- catch non-cumulative observation assigned to observation defined
        !    by a boundname that is assigned to more than one element
        if (obsrv%ObsTypeId == 'STAGE' .or. &
            obsrv%ObsTypeId == 'DEPTH' .or. &
            obsrv%ObsTypeId == 'WET-PERIMETER' .or. &
            obsrv%ObsTypeId == 'WET-AREA' .or. &
            obsrv%ObsTypeId == 'WET-WIDTH') then
          nn1 = obsrv%NodeNumber
          if (nn1 == namedboundflag) then
            if (obsrv%indxbnds_count > 1) then
              write (errmsg, '(a,3(1x,a))') &
                trim(adjustl(obsrv%ObsTypeId)), &
                'for observation', trim(adjustl(obsrv%Name)), &
                ' must be assigned to a reach with a unique boundname.'
              call store_error(errmsg)
            end if
          end if
        end if
        !
        ! -- check that node number 1 is valid; call store_error if not
        do j = 1, obsrv%indxbnds_count
          nn1 = obsrv%indxbnds(j)
          if (nn1 < 1 .or. nn1 > this%maxbound) then
            write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
              trim(adjustl(obsrv%ObsTypeId)), &
              'reach must be greater than 0 and less than or equal to', &
              this%maxbound, '(specified value is ', nn1, ')'
            call store_error(errmsg)
          end if
        end do
      end do
      !
      ! -- evaluate if there are any observation errors
      if (count_errors() > 0) then
        call this%parser%StoreErrorUnit()
      end if
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_set_stressperiod()

subroutine sfrmodule::sfr_set_stressperiod ( class(sfrtype), intent(inout)  this, integer(i4b), intent(in)  n, integer(i4b), intent(inout)  ichkustrm, character(len=linelength), intent(inout)  crossfile  )

private

Method to read and set period data for a SFR package reach.

Parameters

[in,out]	this	SfrType object
[in]	n	reach number
[in,out]	ichkustrm	flag indicating if upstream fraction data specified
[in,out]	crossfile	cross-section file name

Definition at line 3180 of file gwf-sfr.f90.

    ! -- modules
    use timeseriesmanagermodule, only: read_value_or_time_series_adv
    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    integer(I4B), intent(inout) :: ichkustrm !< flag indicating if upstream fraction data specified
    character(len=LINELENGTH), intent(inout) :: crossfile !< cross-section file name
    ! -- local variables
    character(len=10) :: cnum
    character(len=LINELENGTH) :: text
    character(len=LINELENGTH) :: caux
    character(len=LINELENGTH) :: keyword
    integer(I4B) :: ival
    integer(I4B) :: ii
    integer(I4B) :: jj
    integer(I4B) :: idiv
    integer(I4B) :: ixserror
    character(len=10) :: cp
    real(DP) :: divq
    real(DP), pointer :: bndElem => null()
    !
    ! -- initialize variables
    crossfile = 'NONE'
    !
    ! -- read line
    call this%parser%GetStringCaps(keyword)
    select case (keyword)
    case ('STATUS')
      ichkustrm = 1
      call this%parser%GetStringCaps(text)
      if (text == 'INACTIVE') then
        this%iboundpak(n) = 0
      else if (text == 'ACTIVE') then
        this%iboundpak(n) = 1
      else if (text == 'SIMPLE') then
        this%iboundpak(n) = -1
      else
        write (errmsg, '(2a)') &
          'Unknown '//trim(this%text)//' sfr status keyword: ', trim(text)
        call store_error(errmsg)
      end if
    case ('BEDK')
      call this%parser%GetString(text)
      jj = 1 ! For 'BEDK'
      bndelem => this%hk(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, &
                                         'BEDK')
    case ('MANNING')
      call this%parser%GetString(text)
      jj = 1 ! For 'MANNING'
      bndelem => this%rough(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, &
                                         'MANNING')
    case ('STAGE')
      call this%parser%GetString(text)
      jj = 1 ! For 'STAGE'
      bndelem => this%sstage(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, 'STAGE')
    case ('RAINFALL')
      call this%parser%GetString(text)
      jj = 1 ! For 'RAIN'
      bndelem => this%rain(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, 'RAIN')
    case ('EVAPORATION')
      call this%parser%GetString(text)
      jj = 1 ! For 'EVAP'
      bndelem => this%evap(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, &
                                         'EVAP')
    case ('RUNOFF')
      call this%parser%GetString(text)
      jj = 1 ! For 'RUNOFF'
      bndelem => this%runoff(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, &
                                         'RUNOFF')
    case ('INFLOW')
      call this%parser%GetString(text)
      jj = 1 ! For 'INFLOW'
      bndelem => this%inflow(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, &
                                         'INFLOW')
    case ('DIVERSION')
      !
      ! -- make sure reach has at least one diversion
      if (this%ndiv(n) < 1) then
        write (cnum, '(i0)') n
        errmsg = 'diversions cannot be specified for reach '//trim(cnum)
        call store_error(errmsg)
      end if
      !
      ! -- read diversion number
      ival = this%parser%GetInteger()
      if (ival < 1 .or. ival > this%ndiv(n)) then
        write (cnum, '(i0)') n
        errmsg = 'Reach  '//trim(cnum)
        write (cnum, '(i0)') this%ndiv(n)
        errmsg = trim(errmsg)//' diversion number should be between 1 '// &
                 'and '//trim(cnum)//'.'
        call store_error(errmsg)
      end if
      idiv = ival
      !
      ! -- read value
      call this%parser%GetString(text)
      ii = this%iadiv(n) + idiv - 1
      jj = 1 ! For 'DIVERSION'
      bndelem => this%divflow(ii)
      call read_value_or_time_series_adv(text, ii, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, &
                                         'DIVFLOW')
      !
      ! -- if diversion cprior is 'fraction', ensure that 0.0 <= fraction <= 1.0
      cp = this%divcprior(ii)
      divq = this%divflow(ii)
      if (cp == 'FRACTION' .and. (divq < dzero .or. divq > done)) then
        write (errmsg, '(a,1x,i0,a)') &
          'cprior is type FRACTION for diversion no.', ii, &
          ', but divflow not within the range 0.0 to 1.0'
        call store_error(errmsg)
      end if
    case ('UPSTREAM_FRACTION')
      ichkustrm = 1
      call this%parser%GetString(text)
      jj = 1 ! For 'USTRF'
      bndelem => this%ustrf(n)
      call read_value_or_time_series_adv(text, n, jj, bndelem, &
                                         this%packName, 'BND', &
                                         this%tsManager, this%iprpak, 'USTRF')
 
    case ('CROSS_SECTION')
      ixserror = 0
      !
      ! -- read FILE keyword
      call this%parser%GetStringCaps(keyword)
      select case (keyword)
      case ('TAB6')
        call this%parser%GetStringCaps(keyword)
        if (trim(adjustl(keyword)) /= 'FILEIN') then
          errmsg = 'TAB6 keyword must be followed by "FILEIN" '// &
                   'then by filename.'
          call store_error(errmsg)
          ixserror = 1
        end if
        if (ixserror == 0) then
          call this%parser%GetString(crossfile)
        end if
      case default
        write (errmsg, '(a,1x,i4,1x,a)') &
          'CROSS-SECTION TABLE ENTRY for REACH ', n, &
          'MUST INCLUDE TAB6 KEYWORD'
        call store_error(errmsg)
      end select
 
    case ('AUXILIARY')
      call this%parser%GetStringCaps(caux)
      do jj = 1, this%naux
        if (trim(adjustl(caux)) /= trim(adjustl(this%auxname(jj)))) cycle
        call this%parser%GetString(text)
        ii = n
        bndelem => this%rauxvar(jj, ii)
        call read_value_or_time_series_adv(text, ii, jj, bndelem, &
                                           this%packName, 'AUX', &
                                           this%tsManager, this%iprpak, &
                                           this%auxname(jj))
        exit
      end do
 
    case default
      write (errmsg, '(a,a)') &
        'Unknown '//trim(this%text)//' sfr data keyword: ', &
        trim(keyword)//'.'
      call store_error(errmsg)
    end select
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_setup_budobj()

subroutine sfrmodule::sfr_setup_budobj ( class(sfrtype)  this )

private

Method to set up the budget object that stores all the sfr flows The terms listed here must correspond in number and order to the ones listed in the sfr_fill_budobj method.

Parameters

this	SfrType object

Definition at line 5215 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: nbudterm
    integer(I4B) :: i
    integer(I4B) :: n
    integer(I4B) :: n1
    integer(I4B) :: n2
    integer(I4B) :: maxlist
    integer(I4B) :: naux
    integer(I4B) :: idx
    real(DP) :: q
    character(len=LENBUDTXT) :: text
    character(len=LENBUDTXT), dimension(1) :: auxtxt
    !
    ! -- Determine the number of sfr budget terms. These are fixed for
    !    the simulation and cannot change.  This includes FLOW-JA-FACE
    !    so they can be written to the binary budget files, but these internal
    !    flows are not included as part of the budget table.
    nbudterm = 8
    if (this%imover == 1) nbudterm = nbudterm + 2
    if (this%naux > 0) nbudterm = nbudterm + 1
    !
    ! -- set up budobj
    call budgetobject_cr(this%budobj, this%packName)
    call this%budobj%budgetobject_df(this%maxbound, nbudterm, 0, 0, &
                                     ibudcsv=this%ibudcsv)
    idx = 0
    !
    ! -- Go through and set up each budget term
    text = '    FLOW-JA-FACE'
    idx = idx + 1
    maxlist = this%nconn
    naux = 1
    auxtxt(1) = '       FLOW-AREA'
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux, auxtxt)
    !
    ! -- store connectivity
    call this%budobj%budterm(idx)%reset(this%nconn)
    q = dzero
    do n = 1, this%maxbound
      n1 = n
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        n2 = this%ja(i)
        call this%budobj%budterm(idx)%update_term(n1, n2, q)
      end do
    end do
    !
    ! --
    text = '             GWF'
    idx = idx + 1
    maxlist = this%maxbound - this%ianynone
    naux = 1
    auxtxt(1) = '       FLOW-AREA'
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%name_model, &
                                             maxlist, .false., .true., &
                                             naux, auxtxt)
    call this%budobj%budterm(idx)%reset(this%maxbound)
    q = dzero
    do n = 1, this%maxbound
      n2 = this%igwfnode(n)
      if (n2 > 0) then
        call this%budobj%budterm(idx)%update_term(n, n2, q)
      end if
    end do
    !
    ! --
    text = '        RAINFALL'
    idx = idx + 1
    maxlist = this%maxbound
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '     EVAPORATION'
    idx = idx + 1
    maxlist = this%maxbound
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '          RUNOFF'
    idx = idx + 1
    maxlist = this%maxbound
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '      EXT-INFLOW'
    idx = idx + 1
    maxlist = this%maxbound
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '     EXT-OUTFLOW'
    idx = idx + 1
    maxlist = this%maxbound
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '         STORAGE'
    idx = idx + 1
    maxlist = this%maxbound
    naux = 1
    auxtxt(1) = '          VOLUME'
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux, auxtxt)
    !
    ! --
    if (this%imover == 1) then
      !
      ! --
      text = '        FROM-MVR'
      idx = idx + 1
      maxlist = this%maxbound
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux)
      !
      ! --
      text = '          TO-MVR'
      idx = idx + 1
      maxlist = this%maxbound
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux)
    end if
    !
    ! --
    naux = this%naux
    if (naux > 0) then
      !
      ! --
      text = '       AUXILIARY'
      idx = idx + 1
      maxlist = this%maxbound
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux, this%auxname)
    end if
    !
    ! -- if sfr flow for each reach are written to the listing file
    if (this%iprflow /= 0) then
      call this%budobj%flowtable_df(this%iout, cellids='GWF')
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_setup_tableobj()

subroutine sfrmodule::sfr_setup_tableobj ( class(sfrtype)  this )

private

Method to set up the table object that is used to write the sfr stage data. The terms listed here must correspond in number and order to the ones written to the stage table in the sfr_ot method.

Parameters

this	SfrType object

Definition at line 5667 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    ! -- local variables
    integer(I4B) :: nterms
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    !
    ! -- setup stage table
    if (this%iprhed > 0) then
      !
      ! -- Determine the number of sfr budget terms. These are fixed for
      !    the simulation and cannot change.  This includes FLOW-JA-FACE
      !    so they can be written to the binary budget files, but these internal
      !    flows are not included as part of the budget table.
      nterms = 8
      if (this%inamedbound == 1) then
        nterms = nterms + 1
      end if
      !
      ! -- set up table title
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STAGES FOR EACH CONTROL VOLUME'
      !
      ! -- set up stage tableobj
      call table_cr(this%stagetab, this%packName, title)
      call this%stagetab%table_df(this%maxbound, nterms, this%iout, &
                                  transient=.true.)
      !
      ! -- Go through and set up table budget term
      if (this%inamedbound == 1) then
        text = 'NAME'
        call this%stagetab%initialize_column(text, lenboundname, &
                                             alignment=tableft)
      end if
      !
      ! -- reach number
      text = 'NUMBER'
      call this%stagetab%initialize_column(text, 10, alignment=tabcenter)
      !
      ! -- cellids
      text = 'CELLID'
      call this%stagetab%initialize_column(text, 20, alignment=tableft)
      !
      ! -- reach stage
      text = 'STAGE'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- reach depth
      text = 'DEPTH'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- reach width
      text = 'WIDTH'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- gwf head
      text = 'GWF HEAD'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- streambed conductance
      text = 'STREAMBED CONDUCTANCE'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- streambed gradient
      text = 'STREAMBED GRADIENT'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sfr_solve()

subroutine sfrmodule::sfr_solve	(	class(sfrtype)	this,
		integer(i4b), intent(in)	n,
		real(dp), intent(in)	h,
		real(dp), intent(inout)	hcof,
		real(dp), intent(inout)	rhs,
		logical(lgp), intent(in), optional	update
	)

Method to solve the continuity equation for a SFR package reach.

Parameters

	this	SfrType object
[in]	n	reach number
[in]	h	groundwater head in cell connected to reach
[in,out]	hcof	coefficient term added to the diagonal
[in,out]	rhs	right-hand side term
[in]	update	boolean indicating if the reach depth and stage variables should be updated to current iterate

Definition at line 3379 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: h !< groundwater head in cell connected to reach
    real(DP), intent(inout) :: hcof !< coefficient term added to the diagonal
    real(DP), intent(inout) :: rhs !< right-hand side term
    logical(LGP), intent(in), optional :: update !< boolean indicating if the reach depth and stage variables should be updated to current iterate
    ! -- local variables
    logical(LGP) :: lupdate
    integer(I4B) :: i
    integer(I4B) :: ii
    integer(I4B) :: n2
    real(DP) :: hgwf
    real(DP) :: sa
    real(DP) :: sa_wet
    real(DP) :: qu
    real(DP) :: qi
    real(DP) :: qr
    real(DP) :: qe
    real(DP) :: qro
    real(DP) :: qsrc
    real(DP) :: qfrommvr
    real(DP) :: qgwf
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: hsfr
    real(DP) :: qd
    real(DP) :: d1
    real(DP) :: sumleak
    real(DP) :: sumrch
    real(DP) :: gwfhcof
    real(DP) :: gwfrhs
    !
    ! -- Process optional dummy variables
    if (present(update)) then
      lupdate = update
    else
      lupdate = .true.
    end if
 
    ! -- initialize variables
    hcof = dzero
    rhs = dzero
    !
    if (this%iboundpak(n) == 0) then
      this%depth(n) = dzero
      this%stage(n) = dhnoflo
      this%usflow(n) = dzero
      this%simevap(n) = dzero
      this%simrunoff(n) = dzero
      this%dsflow(n) = dzero
      this%gwflow(n) = dzero
    else
      hgwf = h
      d1 = dzero
      qsrc = dzero
      qgwf = dzero
 
      ! -- calculate initial depth assuming a wide cross-section and
      !    ignore groundwater leakage
      ! -- calculate upstream flow
      qu = dzero
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(i) < 0) cycle
        n2 = this%ja(i)
        do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
          if (this%idir(ii) > 0) cycle
          if (this%ja(ii) /= n) cycle
          qu = qu + this%qconn(ii)
        end do
      end do
      this%usflow(n) = qu
 
      ! -- calculate remaining terms
      sa = this%calc_surface_area(n)
      sa_wet = this%calc_surface_area_wet(n, this%depth(n))
      qi = this%inflow(n)
      qr = this%rain(n) * sa
      qe = this%evap(n) * sa_wet
      qro = this%runoff(n)
 
      ! -- Water mover term; assume that it goes in at the upstream end of the reach
      qfrommvr = dzero
      if (this%imover == 1) then
        qfrommvr = this%pakmvrobj%get_qfrommvr(n)
      end if
 
      ! -- calculate downstream flow ignoring groundwater leakage
      qsrc = qu + qi + qr - qe + qro + qfrommvr
 
      ! -- adjust runoff or evaporation if sum of sources is negative
      call this%sfr_adjust_ro_ev(qsrc, qu, qi, qr, qro, qe, qfrommvr)
 
      ! -- set simulated evaporation and runoff
      this%simevap(n) = qe
      this%simrunoff(n) = qro
 
      ! -- calculate reach flow using appropriate method
      if (this%iboundpak(n) < 0) then
        call this%sfr_calc_constant(n, d1, hgwf, qgwf, qd)
      else
        call this%sfr_calc_steady(n, d1, hgwf, qu, qi, &
                                  qfrommvr, qr, qe, qro, &
                                  qgwf, qd)
      end if
 
      ! -- update sfr stage
      tp = this%strtop(n)
      bt = tp - this%bthick(n)
      hsfr = tp + d1
 
      ! -- update stored values
      if (lupdate) then
        ! -- save depth and calculate stage
        this%depth(n) = d1
        this%stage(n) = hsfr
        ! -- update flows
        call this%sfr_update_flows(n, qd, qgwf)
      end if
 
      ! -- calculate sumleak and sumrch
      sumleak = dzero
      sumrch = dzero
      if (this%gwfiss == 0) then
        sumleak = qgwf
      else
        sumleak = qgwf
      end if
      if (hgwf < bt) then
        sumrch = qgwf
      end if
 
      ! -- make final qgwf calculation and obtain
      !    gwfhcof and gwfrhs values
      call this%sfr_calc_qgwf(n, d1, hgwf, qgwf, gwfhcof, gwfrhs)
 
      ! -- update hcof and rhs terms
      if (abs(sumleak) > dzero) then
        ! -- stream leakage is not head dependent
        if (hgwf < bt) then
          rhs = rhs - sumrch
          ! -- stream leakage is head dependent
        else if ((sumleak - qsrc) < -dem30) then
          if (this%gwfiss == 0) then
            rhs = rhs + gwfrhs - sumrch
          else
            rhs = rhs + gwfrhs
          end if
          hcof = gwfhcof
          ! -- place holder for UZF
        else
          if (this%gwfiss == 0) then
            rhs = rhs - sumleak - sumrch
          else
            rhs = rhs - sumleak
          end if
        end if
 
        ! -- add groundwater leakage
      else if (hgwf < bt) then
        rhs = rhs - sumrch
      end if
    end if

sfr_update_flows()

subroutine sfrmodule::sfr_update_flows ( class(sfrtype), intent(inout)  this, integer(i4b), intent(in)  n, real(dp), intent(inout)  qd, real(dp), intent(in)  qgwf  )

private

Method to update downstream flow and groundwater leakage terms for a SFR package reach.

Parameters

[in,out]	this	SfrType object
[in]	n	reach number
[in,out]	qd	downstream reach flow
[in]	qgwf	groundwater leakage for reach

Definition at line 3551 of file gwf-sfr.f90.

    ! -- dummy variables
    class(SfrType), intent(inout) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(inout) :: qd !< downstream reach flow
    real(DP), intent(in) :: qgwf !< groundwater leakage for reach
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: n2
    integer(I4B) :: idiv
    integer(I4B) :: jpos
    real(DP) :: qdiv
    real(DP) :: f
    !
    ! -- update reach terms
    !
    ! -- save final downstream stream flow
    this%dsflow(n) = qd
    !
    ! -- save groundwater leakage
    this%gwflow(n) = qgwf
    !
    ! -- route downstream flow
    if (qd > dzero) then
      !
      ! -- route water to diversions
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(i) > 0) cycle
        idiv = this%idiv(i)
        if (idiv == 0) cycle
        jpos = this%iadiv(n) + idiv - 1
        call this%sfr_calc_div(n, idiv, qd, qdiv)
        this%qconn(i) = qdiv
        this%divq(jpos) = qdiv
      end do
      !
      ! -- Mover terms: store outflow after diversion loss
      !    as qformvr and reduce outflow (qd)
      !    by how much was actually sent to the mover
      if (this%imover == 1) then
        call this%pakmvrobj%accumulate_qformvr(n, qd)
        qd = max(qd - this%pakmvrobj%get_qtomvr(n), dzero)
      end if
      !
      ! -- route remaining water to downstream reaches
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(i) > 0) cycle
        if (this%idiv(i) > 0) cycle
        n2 = this%ja(i)
        if (this%iboundpak(n2) == 0) cycle
        f = this%ustrf(n2) / this%ftotnd(n)
        this%qconn(i) = qd * f
      end do
    else
      do i = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(i) > 0) cycle
        this%qconn(i) = dzero
        idiv = this%idiv(i)
        if (idiv == 0) cycle
        jpos = this%iadiv(n) + idiv - 1
        this%divq(jpos) = dzero
      end do
    end if
    !
    ! -- return
    return

Variable Documentation

ftype

character(len=lenftype) sfrmodule::ftype = 'SFR'

Definition at line 46 of file gwf-sfr.f90.

  character(len=LENFTYPE) :: ftype = 'SFR' !< package ftype string

text

character(len=lenpackagename) sfrmodule::text = ' SFR'

Definition at line 47 of file gwf-sfr.f90.

  character(len=LENPACKAGENAME) :: text = '             SFR' !< package budget string