sfrmodule Module Reference

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

Data Types
type	sfrinputtype
	Input context array pointers for SFR. More...
type	sfrtype

Functions/Subroutines
subroutine, public	sfr_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname, mempath)
	@ 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_checkin_deprecated_names (this)
	Alias deprecated SFR variable names to their current storage. More...
subroutine	sfr_source_options (this)
	Source OPTIONS from input context. More...
subroutine	sfr_source_dimensions (this)
	Source DIMENSIONS block from input context and load static data. More...
subroutine	sfr_source_packagedata (this)
	Source PACKAGEDATA block from input context. More...
subroutine	sfr_source_crosssections (this)
	Source CROSSSECTIONS block from input context. More...
subroutine	sfr_source_connectiondata (this)
	Source CONNECTIONDATA block from input context. More...
subroutine	sfr_source_diversions (this)
	Source DIVERSIONS block from input context. More...
subroutine	sfr_source_initialstages (this)
	Source INITIALSTAGES block from input context. More...
subroutine	sfr_ar (this)
	@ brief Allocate and read method for 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_fp (this)
	Write the Courant-number summary table to the listing file. More...
subroutine	sfr_dt (this)
	Submit the ATS time step for the most Courant-constraining reach. More...
subroutine	sfr_calc_celerity (this, n, q, celerity)
	Kinematic-wave celerity from a flow perturbation (0 if dry) More...
subroutine	sfr_da (this)
	@ brief Deallocate package memory More...
subroutine	sfr_precompute_tvd (this)
	Pre-compute the single upstream reach index for the TVD limiter. 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_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_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_storage_weight (this)
	Check storage weight. 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_initialstages (this)
	Check initial stage data. More...
subroutine	sfr_check_ustrf (this)
	Check upstream fraction data. More...
subroutine	sfr_set_period_value (this, i, n, setting)
	Apply a PERIOD setting value to reach n. 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, tops, flow, gwfhcof, gwfrhs)
	Calculate density terms. More...
subroutine	sfr_input_init (this, mempath, naux)
	Bind SfrInputType PERIOD pointers to input context arrays. More...
subroutine	sfr_input_destroy (this)
	Nullify all SfrInputType pointers. 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 5556 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_area_wet !< wetted area
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local
    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

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 5584 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_perimeter_wet !< wetted perimeter
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local
    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

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 5612 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_surface_area !< surface area
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    ! -- local
    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)

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 5640 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_surface_area_wet !< wetted surface area
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local
    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)

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 5659 of file gwf-sfr.f90.

    ! -- return variable
    real(DP) :: calc_top_width_wet !< wetted top width
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: depth !< reach depth
    ! -- local
    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

Here is the call graph for this function:

define_listlabel()

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

private

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 3059 of file gwf-sfr.f90.

    ! -- dummy
    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

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 5691 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    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))

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 5718 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    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))

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 2033 of file gwf-sfr.f90.

    ! -- modules
    use constantsmodule, only: lenvarname
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: i
    integer(I4B) :: ii
    integer(I4B) :: idiv
    integer(I4B) :: iaux
    character(len=LENVARNAME) :: setting
 
    ! -- sync PACKAGEDATA AUX TS and advance observations;
    ! -- step auxvar period block overrides
    call this%BndExtType%bnd_ad()
    !
    ! -- update previous values
    if (this%istorage == 1) then
      do n = 1, this%maxbound
        this%stageold(n) = this%stage(n)
        this%usinflowold(n) = this%usinflow(n)
        this%dsflowold(n) = this%dsflow(n)
      end do
    end if
    !
    ! -- sync PACKAGEDATA and PERIOD TS fields from input context each
    !    timestep; skipped when no TS6 files are configured. PACKAGEDATA
    !    values are synced only for rows linked to a time series
    if (this%ts_active) then
      do i = 1, this%maxbound
        n = this%pkg_ifno(i)
        if (.not. this%bedk_set(n) .and. &
            this%input%rhk_tslink(i) /= 0) &
          this%hk(n) = this%input%rhk(i)
        if (.not. this%manning_set(n) .and. &
            this%input%man_tslink(i) /= 0) &
          this%rough(n) = this%input%man(i)
        if (.not. this%ustrf_set(n) .and. &
            this%input%ustrf_tslink(i) /= 0) &
          this%ustrf(n) = this%input%ustrf(i)
      end do
      !
      if (this%input%nbound > 0) then
        do i = 1, this%input%nbound
          n = this%input%ifno(i)
          if (n < 1 .or. n > this%maxbound) cycle
          ! -- BEDK, MANNING, UPSTREAM_FRAC, STAGE, INFLOW, RAINFALL, EVAPORATION, RUNOFF
          setting = this%input%setting(i)
          call this%sfr_set_period_value(i, n, setting)
          ! -- DIVFLOW
          if (trim(setting) == 'DIVERSION') then
            if (this%ndiv(n) >= 1) then
              idiv = this%input%idv(i)
              if (idiv >= 1 .and. idiv <= this%ndiv(n)) then
                ii = this%iadiv(n) + idiv - 1
              end if
            end if
          end if
        end do
      end if
    end if
    !
    ! -- check upstream fractions (only needed when USTRF may be TS-updated;
    !    sfr_rp already validates ustrf once per period via ichkustrm)
    if (this%ts_active) call this%sfr_check_ustrf()
    !
    ! -- 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%featureauxvar(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

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 3691 of file gwf-sfr.f90.

    ! -- dummy
    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
    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 483 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B) :: i
    !
    ! -- 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%slope, this%maxbound, 'SLOPE', this%memoryPath)
    call mem_allocate(this%nconnreach, this%maxbound, 'NCONNREACH', &
                      this%memoryPath)
    call mem_allocate(this%bedk_set, this%maxbound, 'BEDK_SET', this%memoryPath)
    call mem_allocate(this%manning_set, this%maxbound, 'MANNING_SET', &
                      this%memoryPath)
    call mem_allocate(this%ustrf_set, this%maxbound, 'USTRF_SET', &
                      this%memoryPath)
    this%bedk_set = .false.
    this%manning_set = .false.
    this%ustrf_set = .false.
    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)
    !
    ! -- stage, usflow, inflow, and dsflow for previous timestep
    if (this%istorage == 1) then
      call mem_allocate(this%stageold, this%maxbound, 'STAGEOLD', &
                        this%memoryPath)
      call mem_allocate(this%usinflow, this%maxbound, 'USINFLOW', &
                        this%memoryPath)
      call mem_allocate(this%usinflowold, this%maxbound, 'USINFLOWOLD', &
                        this%memoryPath)
      call mem_allocate(this%dsflowold, this%maxbound, 'DSFLOWOLD', &
                        this%memoryPath)
      call mem_allocate(this%storage, this%maxbound, 'STORAGE', &
                        this%memoryPath)
      call mem_allocate(this%crmin, this%maxbound, 'CRMIN', this%memoryPath)
      call mem_allocate(this%crmax, this%maxbound, 'CRMAX', this%memoryPath)
      call mem_allocate(this%crsum, this%maxbound, 'CRSUM', this%memoryPath)
      call mem_allocate(this%crcnt, this%maxbound, 'CRCNT', this%memoryPath)
      if (this%ats_courant /= dnodata) then
        call mem_allocate(this%itvd_upstream, this%maxbound, &
                          'ITVD_UPSTREAM', this%memoryPath)
      else
        call mem_allocate(this%itvd_upstream, 0, 'ITVD_UPSTREAM', this%memoryPath)
      end if
    end if
    !
    ! -- 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)
    !
    ! -- alias into the input context's permanent, feature-indexed arrays
    ! (allocated and DZERO-initialized by the loader)
    call mem_setptr(this%rain, 'RAINFALL', this%input_mempath)
    call mem_setptr(this%evap, 'EVAPORATION', this%input_mempath)
    call mem_setptr(this%inflow, 'INFLOW', this%input_mempath)
    call mem_setptr(this%runoff, 'RUNOFF', this%input_mempath)
    call mem_setptr(this%sstage, 'STAGE', this%input_mempath)
    call mem_setptr(this%hk, 'BEDK', this%input_mempath)
    call mem_setptr(this%rough, 'MANNING', this%input_mempath)
    call mem_setptr(this%ustrf, 'UPSTREAM_FRAC', this%input_mempath)
    !
    !
    ! -- 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%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
      !
      ! -- stage
      if (this%istorage == 1) then
        this%stageold(i) = dzero
        this%usinflow(i) = dzero
        this%usinflowold(i) = dzero
        this%dsflowold(i) = dzero
        this%storage(i) = dzero
        this%crmin(i) = dep20
        this%crmax(i) = -dep20
        this%crsum(i) = dzero
        this%crcnt(i) = 0
        if (this%ats_courant /= dnodata) then
          this%itvd_upstream(i) = 0
        end if
      end if
      !
      ! -- boundary data
      this%rough(i) = dzero
      !
      !
      ! -- 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)

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 412 of file gwf-sfr.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    use memoryhelpermodule, only: create_mem_path
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    !
    ! -- call standard BndExtType allocate scalars
    call this%BndExtType%allocate_scalars()
    !
    ! -- allocate the object and assign values to object variables
    call mem_allocate(this%istorage, 'ISTORAGE', this%memoryPath)
    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%storage_weight, 'STORAGE_WEIGHT', this%memoryPath)
    call mem_allocate(this%ats_courant, 'ATS_COURANT', 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%istorage = 0
    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%storage_weight = dnodata
    this%ats_courant = dnodata
    this%nconn = 0
    this%icheck = 1
    this%iconvchk = 1
    this%idense = 0
    this%ivsc = 0
    this%ianynone = 0
    this%ncrossptstot = 0

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 1682 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: inats
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: ierr
    !
    ! -- allocate and read observations
    call this%obs%obs_ar()
    !
    ! -- allocate BndType connection-level arrays (nodelist, bound, auxvar, etc.)
    call this%BndExtType%allocate_arrays()
    !
    ! -- allocate featureauxvar and set pkg_ifno/pkg_auxvar input context
    !    pointers
    call this%BndExtType%allocate_featureauxvar()
    !
    ! -- set persistent PERIOD pointers into input context (valid for
    !    simulation lifetime)
    call this%input%init(this%input_mempath, this%naux)
    !
    ! -- 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 storage_weight
    call this%sfr_check_storage_weight()
    !
    ! -- check that ATS_COURANT is only used with STORAGE
    if (this%ats_courant /= dnodata .and. this%istorage /= 1) then
      write (errmsg, '(a)') &
        'ATS_COURANT OPTION REQUIRES STORAGE OPTION TO BE ACTIVE'
      call store_error(errmsg)
    end if
    !
    ! -- warn when ATS_COURANT is specified but ATS is not active in TDIS
    if (this%ats_courant /= dnodata .and. inats == 0) then
      write (warnmsg, '(a)') &
        'ATS_COURANT IS SPECIFIED IN THE SFR OPTIONS BLOCK BUT THE '// &
        'ATS PACKAGE IS NOT ACTIVE IN TDIS, SO THE TIME STEP IS NOT '// &
        'ADAPTED TO THE COURANT NUMBER. IF THE COURANT NUMBER EXCEEDS 1 '// &
        'THE EXPLICIT KINEMATIC-WAVE ROUTING MAY OSCILLATE WHILE STILL '// &
        'CLOSING THE BUDGET. ACTIVATE THE ATS PACKAGE OR USE A SMALLER '// &
        'TIME STEP.'
      call store_warning(warnmsg)
    end if
    !
    ! -- pre-compute TVD upstream connectivity when ATS_COURANT is active
    if (this%ats_courant /= dnodata .and. this%istorage == 1) then
      call this%sfr_precompute_tvd()
    end if
    !
    ! -- 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
 
    ! -- check the diversion data
    if (this%istorage == 1) then
      call this%sfr_check_initialstages()
    end if
    !
    ! -- terminate if errors were detected in any of the static sfr data
    ierr = count_errors()
    if (ierr > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- setup pakmvrobj
    if (this%imover /= 0) then
      allocate (this%pakmvrobj)
      call this%pakmvrobj%ar(this%maxbound, this%maxbound, this%memoryPath)
    end if

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 3199 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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 store_error_filename(this%input_fname)
      end if
    end if

Here is the call graph for this function:

sfr_calc_celerity()

subroutine sfrmodule::sfr_calc_celerity	(	class(sfrtype)	this,
		integer(i4b), intent(in)	n,
		real(dp), intent(in)	q,
		real(dp), intent(out)	celerity
	)

Parameters

	this	SfrType object
[in]	n	reach number
[in]	q	reach flow
[out]	celerity	kinematic-wave celerity

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

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    real(DP), intent(in) :: q !< reach flow
    real(DP), intent(out) :: celerity !< kinematic-wave celerity
    ! -- local
    real(DP) :: d
    real(DP) :: a
    real(DP) :: a2
    !
    celerity = dzero
    call this%sfr_calc_reach_depth(n, q, d)
    if (d > dzero) then
      a = this%calc_area_wet(n, d)
      call this%sfr_calc_reach_depth(n, q + this%deps, d)
      a2 = this%calc_area_wet(n, d)
      if (a2 > a) then
        celerity = this%deps / (a2 - a)
      end if
    end if

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 3969 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 4016 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 3731 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 3882 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 3812 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 3764 of file gwf-sfr.f90.

    ! -- dummy
    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
    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)

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 4070 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 4105 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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)	tops,
		real(dp), intent(inout)	flow,
		real(dp), intent(inout)	gwfhcof,
		real(dp), intent(inout)	gwfrhs
	)

Method to calculate 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]	tops	top elevation of streambed
[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 5751 of file gwf-sfr.f90.

    ! -- dummy
    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) :: tops !< top elevation of streambed
    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
    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 >= tops) then
      ss = stage
      stage_below_bot = .false.
      rdensesfr = this%denseterms(1, n) ! sfr rel density
    else
      ss = tops
      stage_below_bot = .true.
      rdensesfr = this%denseterms(2, n) ! gwf rel density
    end if
    !
    ! -- set hh to head or tops (top elev of streambed)
    if (head >= tops) then
      hh = head
      head_below_bot = .false.
      rdensegwf = this%denseterms(2, n) ! gwf rel density
    else
      hh = tops
      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 = tops
        elevavg = dhalf * (elevsfr + elevgwf)
        havg = dhalf * (hh + ss)
        d2 = cond * (havg - elevavg) * (rdensegwf - rdensesfr)
        gwfrhs = gwfrhs + d2
        flow = flow + d2
      end if
    end if

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 2309 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: totim, kstp, kper, delt
    ! -- dummy
    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
    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

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 2131 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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

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 4332 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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
    !
    ! -- verify every listed connection is reciprocated: if reach n lists
    ! -- reach nn, reach nn must list reach n back
    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
    !
    ! -- verify upstream connections don't form a mutual pair: two reaches
    ! -- cannot each list the other as upstream
    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 store_error_filename(this%input_fname)
    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

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 4153 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    ! -- 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

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 4608 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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

Here is the call graph for this function:

sfr_check_initialstages()

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

private

Method to check initial data for a SFR package and calculates the initial upstream and downstream flows for the reach based on the initial staalso creates the tables used to print input data, if this option in enabled in the SFR package.

Parameters

this	SfrType object

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

    class(SfrType) :: this !< SfrType object
 
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    character(len=5) :: crch
    integer(I4B) :: n
    real(DP) :: qman
 
    ! skip check if storage is not activated
    if (this%istorage == 0) return
 
    ! write header
    if (this%iprpak /= 0) then
      !
      ! -- reset the input table object
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') REACH INITIAL STAGE DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%maxbound, 4, this%iout)
      text = 'REACH'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'INITIAL STAGE'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'INITIAL DEPTH'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'INITIAL FLOW'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
    end if
    !
    ! -- check that data are correct
    do n = 1, this%maxbound
      write (crch, '(i5)') n
 
      ! calculate the initial flows
      call this%sfr_calc_qman(n, this%depth(n), qman)
      this%usinflow(n) = qman
      this%dsflow(n) = qman
 
      ! add terms to the table
      if (this%iprpak /= 0) then
        call this%inputtab%add_term(n)
        call this%inputtab%add_term(this%stage(n))
        call this%inputtab%add_term(this%depth(n))
        call this%inputtab%add_term(qman)
      end if
    end do

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 4214 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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

Here is the call graph for this function:

sfr_check_storage_weight()

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

private

Method to check the kinematic storage weight for a SFR package. If the kinematic storage weight has not been set it is set to the default value.

Parameters

this	SfrType object

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

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- formats
    character(len=*), parameter :: fmtweight = &
      &"(1x,'SFR PACKAGE (',a,') SETTING DEFAULT',&
      &/4x,'STORAGE_WEIGHT VALUE (',g0,').',/)"
    !
    ! -- set storage weight if it has not been defined yet
    if (this%istorage == 1) then
      if (this%storage_weight == dnodata) then
        this%storage_weight = done
        write (this%iout, fmtweight) &
          trim(adjustl(this%packName)), this%storage_weight
      end if
    end if

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 4770 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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

Here is the call graph for this function:

sfr_checkin_deprecated_names()

subroutine sfrmodule::sfr_checkin_deprecated_names

(

class(sfrtype), intent(inout)

this

)

IDM integration moved several reach variables into the input context, renaming some. These aliases share memory with the relocated variables so that programs addressing them through the API by their old package-path names continue to resolve; each use is reported as deprecated. DIVFLOW is only set from the input context when diversions are present; otherwise an empty DIVFLOW is allocated at the package path, as before IDM integration.

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

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr, get_isize
    use memoryaliasmodule, only: mem_alias_deprecated
    ! -- dummy
    class(SfrType), intent(inout) :: this
    ! -- local
    real(DP), dimension(:, :), pointer, contiguous :: aux => null()
    integer(I4B) :: isize
    ! -- code
    !
    ! -- reach auxiliary values, present (possibly empty) in the input context
    call get_isize('AUX', this%input_mempath, isize)
    if (isize >= 0) then
      call mem_setptr(aux, 'AUX', this%input_mempath)
      call mem_alias_deprecated(aux, 'RAUXVAR', this%memoryPath, &
                                'AUX', this%input_mempath, '6.8.2')
    end if
    !
    ! -- diversion flows
    if (associated(this%divflow)) then
      call mem_alias_deprecated(this%divflow, 'DIVFLOW', this%memoryPath, &
                                'DIVFLOW', this%input_mempath, '6.8.2')
    else
      call mem_allocate(this%divflow, 0, 'DIVFLOW', this%memoryPath)
    end if
    !
    call mem_alias_deprecated(this%rain, 'RAIN', this%memoryPath, &
                              'RAINFALL', this%input_mempath, '6.8.2')
    call mem_alias_deprecated(this%evap, 'EVAP', this%memoryPath, &
                              'EVAPORATION', this%input_mempath, '6.8.2')
    call mem_alias_deprecated(this%inflow, 'INFLOW', this%memoryPath, &
                              'INFLOW', this%input_mempath, '6.8.2')
    call mem_alias_deprecated(this%runoff, 'RUNOFF', this%memoryPath, &
                              'RUNOFF', this%input_mempath, '6.8.2')
    call mem_alias_deprecated(this%sstage, 'SSTAGE', this%memoryPath, &
                              'STAGE', this%input_mempath, '6.8.2')
    call mem_alias_deprecated(this%hk, 'HK', this%memoryPath, &
                              'BEDK', this%input_mempath, '6.8.2')
    call mem_alias_deprecated(this%rough, 'ROUGH', this%memoryPath, &
                              'MANNING', this%input_mempath, '6.8.2')
    call mem_alias_deprecated(this%ustrf, 'USTRF', this%memoryPath, &
                              'UPSTREAM_FRAC', this%input_mempath, '6.8.2')

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 2510 of file gwf-sfr.f90.

    ! -- modules
    use budgetmodule, only: budgettype
    ! -- dummy
    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
    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()

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,
		character(len=*), intent(in)	mempath
	)

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
[in]	mempath	input mempath

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

    ! -- modules
    use memoryhelpermodule, only: create_mem_path
    ! -- dummy
    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
    character(len=*), intent(in) :: mempath !< input mempath
    ! -- local
    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, mempath)
    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')

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 )

private

Deallocate SFR package scalars and arrays.

Parameters

this	SfrType object

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

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    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%slope)
    call mem_deallocate(this%nconnreach)
    call mem_deallocate(this%bedk_set)
    call mem_deallocate(this%manning_set)
    call mem_deallocate(this%ustrf_set)
    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)
    !
    ! -- stage, usflow, and dsflow for previous timestep
    if (this%istorage == 1) then
      call mem_deallocate(this%stageold)
      call mem_deallocate(this%dsflowold)
      call mem_deallocate(this%storage)
      call mem_deallocate(this%usinflow)
      call mem_deallocate(this%usinflowold)
      call mem_deallocate(this%crmin)
      call mem_deallocate(this%crmax)
      call mem_deallocate(this%crsum)
      call mem_deallocate(this%crcnt)
      call mem_deallocate(this%itvd_upstream)
    end if
    !
    ! -- 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)
    !
    nullify (this%rain)
    nullify (this%evap)
    nullify (this%inflow)
    nullify (this%runoff)
    nullify (this%sstage)
    nullify (this%hk)
    nullify (this%rough)
    nullify (this%ustrf)
    !
    ! -- deallocate diversion variables
    call mem_deallocate(this%iadiv)
    call mem_deallocate(this%divreach)
    if (associated(this%divcprior)) then
      deallocate (this%divcprior)
    end if
    ! -- divflow is package-owned only when empty (no diversions); it is not
    !    associated if the simulation stopped before it was set
    if (associated(this%divflow)) then
      if (size(this%divflow) == 0) then
        call mem_deallocate(this%divflow)
      else
        nullify (this%divflow)
      end if
    end if
    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 Courant number table
    if (associated(this%couranttab)) then
      call this%couranttab%table_da()
      deallocate (this%couranttab)
      nullify (this%couranttab)
    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%ats_courant)
    call mem_deallocate(this%istorage)
    call mem_deallocate(this%storage_weight)
    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)
    !
    ! -- nullify SFR input context pointers (input context owns the memory)
    call this%input%destroy()
    !
    ! -- call base BndExtType deallocate
    call this%BndExtType%bnd_da()

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 3103 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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

sfr_dt()

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

private

Parameters

this	SfrType object

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

    ! -- modules
    use tdismodule, only: kstp, kper, ats
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: nrmin
    real(DP) :: celerity
    real(DP) :: dt_n
    real(DP) :: dtmin
    character(len=LINELENGTH) :: msg
    !
    if (this%ats_courant == dnodata) return
    if (this%istorage /= 1) return
    !
    dtmin = dnodata
    nrmin = 0
    !
    do n = 1, this%maxbound
      call this%sfr_calc_celerity(n, this%dsflow(n), celerity)
      if (celerity > dzero) then
        dt_n = this%ats_courant * this%length(n) / celerity
        if (dt_n < dtmin) then
          dtmin = dt_n
          nrmin = n
        end if
      end if
    end do
    !
    if (nrmin > 0) then
      write (msg, '(a,i0)') trim(this%packName)//'-REACH-', nrmin
      call ats%ats_submit_delt(kstp, kper, dtmin, trim(msg))
    end if

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 2159 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 5210 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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
        if (this%gwfiss == 0 .and. this%istorage == 1) then
          q = this%storage(n)
        end if
      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()

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 2259 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

sfr_fp()

subroutine sfrmodule::sfr_fp

(

class(sfrtype)

this

)

Parameters

this	SfrType object

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

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    integer(I4B) :: n
    real(DP) :: crmean
    !
    if (this%istorage == 1) then
      do n = 1, this%maxbound
        if (this%inamedbound == 1) then
          call this%couranttab%add_term(this%boundname(n))
        end if
        call this%couranttab%add_term(n)
        if (this%crmin(n) == dep20) then
          call this%couranttab%add_term('--')
        else
          call this%couranttab%add_term(this%crmin(n))
        end if
        if (this%crmax(n) < dzero) then
          call this%couranttab%add_term('--')
        else
          call this%couranttab%add_term(this%crmax(n))
        end if
        if (this%crcnt(n) > 0) then
          crmean = this%crsum(n) / real(this%crcnt(n), dp)
          call this%couranttab%add_term(crmean)
        else
          call this%couranttab%add_term('--')
        end if
      end do
    end if

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 3949 of file gwf-sfr.f90.

    ! -- return variable
    integer(I4B) :: sfr_gwf_conn !< flag indicating if reach is connected to a gwf cell
    ! -- dummy
    class(SfrType) :: this !< SfrType object
    integer(I4B), intent(in) :: n !< reach number
    ! -- local
    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_input_destroy()

subroutine sfrmodule::sfr_input_destroy

(

class(sfrinputtype), intent(inout)

this

)

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

    class(SfrInputType), intent(inout) :: this
    nullify (this%rhk)
    nullify (this%man)
    nullify (this%ustrf)
    nullify (this%rhk_tslink)
    nullify (this%man_tslink)
    nullify (this%ustrf_tslink)
    nullify (this%nbound)
    nullify (this%ifno)
    nullify (this%status)
    nullify (this%setting)
    nullify (this%idv)
    nullify (this%tab6_filename)
    nullify (this%auxname)
    nullify (this%auxval)

sfr_input_init()

subroutine sfrmodule::sfr_input_init ( class(sfrinputtype), intent(inout)  this, character(len=*), intent(in)  mempath, integer(i4b), intent(in)  naux  )

private

Called once from sfr_ar; PACKAGEDATA pointers are set separately in sfr_source_packagedata.

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

    use memorymanagermodule, only: mem_setptr
    class(SfrInputType), intent(inout) :: this
    character(len=*), intent(in) :: mempath
    integer(I4B), intent(in) :: naux
    call mem_setptr(this%nbound, 'NBOUND', mempath)
    call mem_setptr(this%ifno, 'IFNO', mempath)
    call mem_setptr(this%setting, 'SETTING', mempath)
    call mem_setptr(this%status, 'STATUS', mempath)
    call mem_setptr(this%idv, 'IDV', mempath)
    call mem_setptr(this%tab6_filename, 'TAB6_FILENAME', mempath)
    if (naux > 0) then
      call mem_setptr(this%auxname, 'AUXNAME', mempath)
      call mem_setptr(this%auxval, 'AUXVAL', mempath)
    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 3091 of file gwf-sfr.f90.

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

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 2752 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)

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 2640 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: kstp, kper, pertim, totim
    use inputoutputmodule, only: ulasav
    ! -- dummy
    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
    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

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 2587 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: kstp, kper, delt, pertim, totim
    ! -- dummy
    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
    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

sfr_precompute_tvd()

subroutine sfrmodule::sfr_precompute_tvd

(

class(sfrtype), intent(inout)

this

)

Index 0 for reaches with zero or multiple upstream connections.

Parameters

[in,out]

this

SfrType object

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

    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: j
    integer(I4B) :: iup_count
    integer(I4B) :: m
    !
    do n = 1, this%maxbound
      iup_count = 0
      m = 0
      do j = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(j) > 0) then
          iup_count = iup_count + 1
          m = this%ja(j)
        end if
      end do
      if (iup_count == 1) then
        this%itvd_upstream(n) = m
      else
        this%itvd_upstream(n) = 0
      end if
    end do

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 3401 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

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 1783 of file gwf-sfr.f90.

    ! -- modules
    use constantsmodule, only: lenvarname
    use tdismodule, only: kper
    use memorymanagermodule, only: mem_reallocate
    use sfrcrosssectionmanager, only: cross_section_cr, sfrcrosssection
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    character(len=LINELENGTH) :: crossfile
    character(len=LINELENGTH) :: str
    real(DP) :: auxvaldbl
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    character(len=LENVARNAME) :: setting
    character(len=10) :: cnum
    character(len=10) :: cp
    integer(I4B) :: n
    integer(I4B) :: i
    integer(I4B) :: ii
    integer(I4B) :: jj
    integer(I4B) :: idiv
    integer(I4B) :: ichkustrm
    integer(I4B) :: ncrossptstot
    type(SfrCrossSection), pointer :: cross_data => null()
    ! -- formats
    character(len=*), parameter :: fmtlsp = &
      &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
    !
    ! -- initialize flags
    ichkustrm = 0
    if (kper == 1) then
      ichkustrm = 1
    end if
    !
    ! -- set nbound to maxbound
    this%nbound = this%maxbound
    !
    ! -- check if input context has been updated for this stress period
    if (this%iper == 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 to echo period data
      if (this%iprpak /= 0) then
        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 jj = 1, 2
          write (text, '(a,1x,i6)') 'VALUE', jj
          call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
        end do
      end if
      !
      ! -- process each period data entry
      do i = 1, this%input%nbound
        n = this%input%ifno(i)
        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
        !
        ! -- dispatch key for this row
        setting = this%input%setting(i)
        !
        ! -- STATUS
        if (trim(setting) == 'STATUS') then
          str = this%input%status(i)
          ichkustrm = 1
          if (str == 'INACTIVE') then
            this%iboundpak(n) = 0
          else if (str == 'ACTIVE') then
            this%iboundpak(n) = 1
          else if (str == 'SIMPLE') then
            this%iboundpak(n) = -1
          else
            write (errmsg, '(2a)') &
              'Unknown '//trim(this%text)//' sfr status keyword: ', trim(str)
            call store_error(errmsg)
          end if
        end if
        !
        ! -- BEDK, MANNING, STAGE, INFLOW, RAINFALL, EVAPORATION, RUNOFF,
        !    UPSTREAM_FRAC; sfr_ad re-syncs these each timestep if TS-linked
        call this%sfr_set_period_value(i, n, setting)
        if (trim(setting) == 'UPSTREAM_FRAC') ichkustrm = 1
        !
        ! -- DIVERSION (compound group)
        if (trim(setting) == 'DIVERSION') then
          if (this%ndiv(n) < 1) then
            write (cnum, '(i0)') n
            errmsg = 'diversions cannot be specified for reach '//trim(cnum)
            call store_error(errmsg)
          else
            idiv = this%input%idv(i)
            if (idiv < 1 .or. idiv > 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)
            else
              ii = this%iadiv(n) + idiv - 1
              ! -- if cprior is 'FRACTION', check 0.0 <= divflow <= 1.0
              cp = this%divcprior(ii)
              if (cp == 'FRACTION' .and. &
                  (this%divflow(ii) < dzero .or. this%divflow(ii) > 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
            end if
          end if
        end if
        !
        ! -- CROSS_SECTION (compound group)
        if (trim(setting) == 'CROSS_SECTION') then
          crossfile = this%input%tab6_filename(i)
          call cross_data%read_table(n, this%width(n), trim(adjustl(crossfile)))
        end if
        !
        ! -- AUXILIARY (compound group): featureauxvar already holds the
        ! -- current value; only extract auxname for echoing this row
        if (this%naux > 0) then
          if (trim(setting) == 'PERIOD_AUXILIARY') then
            str = this%input%auxname(i)
            auxvaldbl = this%input%auxval(i)
          end if
        end if
        !
        ! -- echo row to period data table
        if (this%iprpak /= 0) then
          call this%inputtab%add_term(n)
          if (trim(setting) == 'PERIOD_AUXILIARY') then
            call this%inputtab%add_term('AUXILIARY')
          else
            call this%inputtab%add_term(trim(setting))
          end if
          select case (trim(setting))
          case ('STATUS')
            call this%inputtab%add_term(trim(str))
            call this%inputtab%add_term(' ')
          case ('BEDK')
            call this%inputtab%add_term(this%hk(n))
            call this%inputtab%add_term(' ')
          case ('MANNING')
            call this%inputtab%add_term(this%rough(n))
            call this%inputtab%add_term(' ')
          case ('STAGE')
            call this%inputtab%add_term(this%sstage(n))
            call this%inputtab%add_term(' ')
          case ('INFLOW')
            call this%inputtab%add_term(this%inflow(n))
            call this%inputtab%add_term(' ')
          case ('RAINFALL')
            call this%inputtab%add_term(this%rain(n))
            call this%inputtab%add_term(' ')
          case ('EVAPORATION')
            call this%inputtab%add_term(this%evap(n))
            call this%inputtab%add_term(' ')
          case ('RUNOFF')
            call this%inputtab%add_term(this%runoff(n))
            call this%inputtab%add_term(' ')
          case ('UPSTREAM_FRAC')
            call this%inputtab%add_term(this%ustrf(n))
            call this%inputtab%add_term(' ')
          case ('DIVERSION')
            call this%inputtab%add_term(idiv)
            call this%inputtab%add_term(this%divflow(ii))
          case ('CROSS_SECTION')
            call this%inputtab%add_term(trim(crossfile))
            call this%inputtab%add_term(' ')
          case ('PERIOD_AUXILIARY')
            call this%inputtab%add_term(trim(str))
            call this%inputtab%add_term(auxvaldbl)
          case default
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          end select
        end if
      end do
      if (this%iprpak /= 0) then
        call this%inputtab%finalize_table()
      end if
      !
      ! -- determine the current size of cross-section data
      ncrossptstot = cross_data%get_ncrossptstot()
      !
      ! -- reallocate sfr package cross-section data if size changed
      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)
    else
      write (this%iout, fmtlsp) trim(this%filtyp)
    end if
    !
    ! -- check upstream fraction values every period; ftotnd must be
    !    initialized for sfr_update_flows even when there is no period data
    if (ichkustrm /= 0) then
      call this%sfr_check_ustrf()
    end if
    !
    ! -- write summary of package block error messages
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if

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 3291 of file gwf-sfr.f90.

    ! -- modules
    use tdismodule, only: kper
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    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 store_error_filename(this%input_fname)
      end if
    end if

Here is the call graph for this function:

sfr_set_period_value()

subroutine sfrmodule::sfr_set_period_value ( class(sfrtype), intent(inout)  this, integer(i4b), intent(in)  i, integer(i4b), intent(in)  n, character(len=*), intent(in)  setting  )

private

Shared by sfr_rp and sfr_ad; DIVERSION is handled separately by both callers because its target index requires its own validation.

Parameters

[in]	i	row index into input context arrays
[in]	n	reach number
[in]	setting	uppercase mf6varname dispatch key

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

    ! -- dummy
    class(SfrType), intent(inout) :: this
    integer(I4B), intent(in) :: i !< row index into input context arrays
    integer(I4B), intent(in) :: n !< reach number
    character(len=*), intent(in) :: setting !< uppercase mf6varname dispatch key
    !
    select case (trim(setting))
    case ('BEDK')
      this%bedk_set(n) = .true.
    case ('MANNING')
      this%manning_set(n) = .true.
    case ('UPSTREAM_FRAC')
      this%ustrf_set(n) = .true.
    end select

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 4995 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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(maxlist)
    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

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 5442 of file gwf-sfr.f90.

    ! -- dummy
    class(SfrType) :: this !< SfrType object
    ! -- local
    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
    !
    ! -- setup Courant number table
    if (this%istorage == 1) then
      nterms = 4
      if (this%inamedbound == 1) then
        nterms = nterms + 1
      end if
      !
      ! -- set up table title
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') COURANT NUMBER FOR EACH REACH'
      !
      ! -- set up Courant tableobj
      call table_cr(this%couranttab, this%packName, title)
      call this%couranttab%table_df(this%maxbound, nterms, this%iout)
      !
      if (this%inamedbound == 1) then
        text = 'NAME'
        call this%couranttab%initialize_column(text, lenboundname, &
                                               alignment=tableft)
      end if
      !
      ! -- reach number
      text = 'NUMBER'
      call this%couranttab%initialize_column(text, 10, alignment=tabcenter)
      !
      ! -- minimum Courant number
      text = 'MINIMUM'
      call this%couranttab%initialize_column(text, 16, alignment=tabcenter)
      !
      ! -- maximum Courant number
      text = 'MAXIMUM'
      call this%couranttab%initialize_column(text, 16, alignment=tabcenter)
      !
      ! -- mean Courant number
      text = 'MEAN'
      call this%couranttab%initialize_column(text, 16, alignment=tabcenter)
    end if

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  )

private

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 3437 of file gwf-sfr.f90.

    ! -- dummy
    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
    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
      qd = this%dsflow(n)
 
      ! -- 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
        if (this%gwfiss == 0 .and. this%istorage == 1) then
          if (this%ats_courant /= dnodata) then
            call this%sfr_calc_tvd(n, d1, hgwf, qu, qi, &
                                   qfrommvr, qr, qe, qro, &
                                   qgwf, qd)
          else
            call this%sfr_calc_transient(n, d1, hgwf, qu, qi, &
                                         qfrommvr, qr, qe, qro, &
                                         qgwf, qd)
          end if
        else
          call this%sfr_calc_steady(n, d1, hgwf, qu, qi, &
                                    qfrommvr, qr, qe, qro, &
                                    qgwf, qd)
        end if
      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_source_connectiondata()

subroutine sfrmodule::sfr_source_connectiondata

(

class(sfrtype), intent(inout)

this

)

Parameters

[in,out]

this

SfrType object

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

    ! -- modules
    use memorymanagermodule, only: mem_setptr, mem_reallocate, get_isize
    use sparsemodule, only: sparsematrix
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B), dimension(:), pointer, contiguous :: cd_ifno => null()
    integer(I4B), dimension(:), pointer, contiguous :: ic_flat => null()
    integer(I4B), dimension(:), pointer, contiguous :: cd_ia => null()
    type(sparsematrix), pointer :: sparse => null()
    integer(I4B), dimension(:), pointer, contiguous :: rowmaxnnz => null()
    integer(I4B), allocatable :: nboundchk(:)
    integer(I4B), allocatable :: iconndata(:, :)
    integer(I4B), allocatable :: iup(:)
    integer(I4B), allocatable :: order(:)
    type(dag) :: sfr_dag
    integer(I4B) :: i, n, j, jj, jcol, jcol2, nja, nconnmax
    integer(I4B) :: ival, idir, nup, ipos, istat, ierr, jlocal, nread
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' CONNECTIONDATA'
    !
    ! -- get input context CONNECTIONDATA arrays
    call mem_setptr(cd_ifno, 'CONNDATA_IFNO', this%input_mempath)
    call mem_setptr(ic_flat, 'IC', this%input_mempath)
    call mem_setptr(cd_ia, 'IC_IA', this%input_mempath)
    !
    ! -- allocate and initialize check array
    allocate (nboundchk(this%maxbound))
    nboundchk = 0
    !
    ! -- build rowmaxnnz from nconnreach
    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) nconnmax = ival
    end do
    !
    ! -- reallocate connection arrays
    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)
    this%idir = 0
    this%idiv = 0
    this%qconn = dzero
    !
    ! -- allocate iconndata
    allocate (iconndata(nconnmax, this%maxbound))
    iconndata = 0
    !
    allocate (sparse)
    call sparse%init(this%maxbound, this%maxbound, rowmaxnnz)
    !
    ! -- process connection data rows
    do i = 1, this%maxbound
      n = this%validate_ifno(cd_ifno(i), this%maxbound, nboundchk, &
                             'reach', 'CONNECTIONDATA')
      if (n == 0) cycle
      nread = cd_ia(i + 1) - cd_ia(i)
      if (nread /= this%nconnreach(n)) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,a)') &
          'Reach', n, 'connectiondata specifies', nread, &
          'connections but NCON for this reach in packagedata is', &
          this%nconnreach(n), '.'
        call store_error(errmsg)
        cycle
      end if
      call sparse%addconnection(n, n, 1)
      jlocal = 0
      do j = cd_ia(i), cd_ia(i + 1) - 1
        jlocal = jlocal + 1
        ival = ic_flat(j)
        iconndata(jlocal, n) = ival
        if (ival < 0) then
          idir = -1
          ival = abs(ival)
        else if (ival == 0) then
          call store_error('Missing or zero connection reach in connectiondata.')
          cycle
        else
          idir = 1
        end if
        if (ival > this%maxbound) then
          call store_error('Reach number exceeds NREACHES in connectiondata.')
          cycle
        end if
        call sparse%addconnection(n, ival, 1)
      end do
    end do
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' CONNECTIONDATA'
    !
    call this%report_ifno_coverage(nboundchk, this%maxbound, 'reach', &
                                   'CONNECTIONDATA')
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- fill ia and ja from sparse
    call sparse%filliaja(this%ia, this%ja, ierr, sort=.true.)
    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 with direction
    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) idir = -1
            this%idir(j) = idir
            exit
          end if
        end do
      end do
    end do
    !
    ! -- deallocate temporaries
    deallocate (rowmaxnnz)
    deallocate (nboundchk)
    deallocate (iconndata)
    call sparse%destroy()
    deallocate (sparse)
    !
    ! -- compute reach order using DAG
    call sfr_dag%set_vertices(this%maxbound)
    fill_dag: do n = 1, this%maxbound
      nup = 0
      do j = this%ia(n) + 1, this%ia(n + 1) - 1
        if (this%idir(j) > 0) nup = nup + 1
      end do
      if (nup == 0) cycle fill_dag
      allocate (iup(nup))
      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
      call sfr_dag%set_edges(n, iup)
      deallocate (iup)
    end do fill_dag
    call sfr_dag%toposort(order, istat)
    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
    do n = 1, this%maxbound
      if (istat == 0) then
        this%isfrorder(n) = order(n)
      else
        this%isfrorder(n) = n
      end if
    end do
    call sfr_dag%destroy()
    if (istat == 0) deallocate (order)
    !
    ! -- release input context connectiondata memory
    call memorystore_release('CONNDATA_IFNO', this%input_mempath)
    call memorystore_release('IC', this%input_mempath)
    call memorystore_release('IC_IA', this%input_mempath)

Here is the call graph for this function:

sfr_source_crosssections()

subroutine sfrmodule::sfr_source_crosssections

(

class(sfrtype), intent(inout)

this

)

Parameters

[in,out]

this

SfrType object

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

    ! -- modules
    use constantsmodule, only: lenmempath
    use memorymanagermodule, only: mem_setptr, mem_reallocate, get_isize
    use characterstringmodule, only: characterstringtype
    use sfrcrosssectionmanager, only: cross_section_cr, sfrcrosssection
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    type(CharacterStringType), dimension(:), pointer, contiguous :: &
      tab6_filename => null()
    type(CharacterStringType), dimension(:), pointer, contiguous :: &
      sfrtab6_mempaths => null()
    type(SfrCrossSection), pointer :: cross_data => null()
    integer(I4B) :: i, n, ncrossptstot
    integer(I4B), allocatable :: nboundchk(:)
    integer(I4B) :: isize
    character(len=LINELENGTH) :: tabfname
    character(len=LENMEMPATH) :: mempath
    !
    ! -- check if CROSSSECTIONS block was loaded (optional block)
    call get_isize('CROSSSECT_IFNO', this%input_mempath, isize)
    if (isize < 1) return
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' CROSSSECTIONS'
    !
    ! -- get input context arrays
    call mem_setptr(ifno, 'CROSSSECT_IFNO', this%input_mempath)
    call mem_setptr(tab6_filename, 'SFRTAB6_FILENAME', this%input_mempath)
    call mem_setptr(sfrtab6_mempaths, 'SFRTAB6_MEMPATH', this%input_mempath)
    !
    ! -- create cross-section object
    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)
    !
    allocate (nboundchk(this%maxbound))
    nboundchk = 0
    !
    do i = 1, isize
      n = this%validate_ifno(ifno(i), this%maxbound, nboundchk, &
                             'reach', 'CROSSSECTIONS')
      if (n == 0) cycle
      tabfname = tab6_filename(i)
      mempath = sfrtab6_mempaths(i)
      call cross_data%source_table(n, this%width(n), mempath, &
                                   trim(adjustl(tabfname)))
    end do
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' CROSSSECTIONS'
    !
    ! -- CROSSSECTIONS is optional; only duplicates are reported
    call this%report_ifno_coverage(nboundchk, this%maxbound, 'reach', &
                                   'CROSSSECTIONS', check_missing=.false.)
    deallocate (nboundchk)
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- reallocate cross-section arrays if needed
    ncrossptstot = cross_data%get_ncrossptstot()
    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
    !
    call cross_data%output(this%width, this%rough)
    call cross_data%pack(this%ncrossptstot, this%ncrosspts, this%iacross, &
                         this%station, this%xsheight, this%xsrough)
    call cross_data%destroy()
    deallocate (cross_data)
    nullify (cross_data)
    !
    ! -- release input context crosssections memory
    call memorystore_release('CROSSSECT_IFNO', this%input_mempath)
    call memorystore_release('SFRTAB6_FILENAME', this%input_mempath)

Here is the call graph for this function:

sfr_source_dimensions()

subroutine sfrmodule::sfr_source_dimensions

(

class(sfrtype), intent(inout)

this

)

Parameters

[in,out]

this

SfrType object

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

    ! -- modules
    use memorymanagerextmodule, only: mem_set_value
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    logical(LGP) :: found
    !
    ! -- initialize dimensions to 0
    this%maxbound = 0
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' DIMENSIONS'
    !
    ! -- source NREACHES
    call mem_set_value(this%maxbound, 'MAXBOUND', this%input_mempath, found, &
                       release=.false.)
    if (found) then
      write (this%iout, '(4x,a,i0)') 'NREACHES = ', this%maxbound
    end if
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
    !
    if (this%maxbound < 1) then
      write (errmsg, '(a)') &
        'NREACHES was not specified or was specified incorrectly.'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- set npakeq (SFR adds no extra equations)
    this%npakeq = 0
    !
    ! -- define the default label
    call this%define_listlabel()
    !
    ! -- define default cross-section data size
    this%ncrossptstot = this%maxbound
    !
    ! -- allocate arrays
    call this%sfr_allocate_arrays()
    !
    ! -- source packagedata
    call this%sfr_source_packagedata()
    !
    ! -- source crosssections
    call this%sfr_source_crosssections()
    !
    ! -- source connectiondata
    call this%sfr_source_connectiondata()
    !
    ! -- source diversions
    call this%sfr_source_diversions()
    !
    ! -- alias deprecated variable names
    call this%sfr_checkin_deprecated_names()
    !
    ! -- source initialstages
    call this%sfr_source_initialstages()
    !
    ! -- setup the budget object
    call this%sfr_setup_budobj()
    !
    ! -- setup the stage table object
    call this%sfr_setup_tableobj()

Here is the call graph for this function:

sfr_source_diversions()

subroutine sfrmodule::sfr_source_diversions

(

class(sfrtype), intent(inout)

this

)

Parameters

[in,out]

this

SfrType object

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

    ! -- modules
    use memorymanagermodule, only: mem_setptr, mem_reallocate, get_isize
    use characterstringmodule, only: characterstringtype
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B), dimension(:), pointer, contiguous :: d_ifno => null()
    integer(I4B), dimension(:), pointer, contiguous :: d_idv => null()
    integer(I4B), dimension(:), pointer, contiguous :: d_iconr => null()
    type(CharacterStringType), dimension(:), pointer, contiguous :: &
      d_cprior => null()
    integer(I4B) :: i, n, idiv, ndiversions, i0, ipos, jpos
    integer(I4B) :: isize
    integer(I4B), allocatable :: iachk(:), nboundchk(:)
    character(len=16) :: cval
    !
    ! -- build iadiv and count total diversions
    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 diversion arrays
    if (ndiversions > 0) then
      call mem_reallocate(this%divreach, ndiversions, 'DIVREACH', &
                          this%memoryPath)
      allocate (this%divcprior(ndiversions))
      call mem_setptr(this%divflow, 'DIVFLOW', this%input_mempath)
      call mem_reallocate(this%divq, ndiversions, 'DIVQ', this%memoryPath)
      this%divq = dzero
    end if
    !
    ! -- check if DIVERSIONS block was loaded (optional block)
    call get_isize('DIVERSIONS_IFNO', this%input_mempath, isize)
    if (isize < 1) then
      if (this%idiversions /= 0) then
        call store_error('REQUIRED DIVERSIONS BLOCK NOT FOUND.')
        call store_error_filename(this%input_fname)
      end if
      return
    end if
    !
    if (this%idiversions == 0) then
      write (errmsg, '(a,1x,a)') &
        'A diversions block should not be', &
        'specified if diversions are not specified.'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
      return
    end if
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' DIVERSIONS'
    !
    ! -- get input context arrays
    call mem_setptr(d_ifno, 'DIVERSIONS_IFNO', this%input_mempath)
    call mem_setptr(d_idv, 'DIVERSIONS_IDV', this%input_mempath)
    call mem_setptr(d_iconr, 'ICONR', this%input_mempath)
    call mem_setptr(d_cprior, 'CPRIOR', this%input_mempath)
    !
    ! -- build check arrays
    allocate (iachk(this%maxbound + 1))
    allocate (nboundchk(ndiversions))
    nboundchk = 0
    iachk(1) = 1
    do n = 1, this%maxbound
      iachk(n + 1) = iachk(n) + this%ndiv(n)
    end do
    !
    do i = 1, isize
      n = d_ifno(i)
      if (n < 1 .or. n > this%maxbound) then
        write (errmsg, '(a,i0,a)') 'Reach number should be between 1 and ', &
          this%maxbound, '.'
        call store_error(errmsg)
        cycle
      end if
      if (this%ndiv(n) < 1) then
        write (errmsg, '(a,i0)') &
          'Diversions cannot be specified for reach ', n
        call store_error(errmsg)
        cycle
      end if
      idiv = d_idv(i)
      if (idiv < 1 .or. idiv > this%ndiv(n)) then
        write (errmsg, '(a,i0,a,i0,a)') &
          'Reach ', n, ' diversion number should be between 1 and ', &
          this%ndiv(n), '.'
        call store_error(errmsg)
        cycle
      end if
      ipos = iachk(n) + idiv - 1
      nboundchk(ipos) = nboundchk(ipos) + 1
      jpos = this%iadiv(n) + idiv - 1
      if (d_iconr(i) < 1 .or. d_iconr(i) > this%maxbound) then
        write (errmsg, '(a,i0,a)') &
          'Diversion target reach number should be between 1 and ', &
          this%maxbound, '.'
        call store_error(errmsg)
        cycle
      end if
      this%divreach(jpos) = d_iconr(i)
      cval = d_cprior(i)
      select case (trim(cval))
      case ('UPTO', 'THRESHOLD', 'FRACTION', 'EXCESS')
      case default
        errmsg = 'Invalid cprior type '//trim(cval)//'.'
        call store_error(errmsg)
      end select
      this%divcprior(jpos) = trim(cval)
    end do
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' DIVERSIONS'
    !
    do n = 1, this%maxbound
      do idiv = 1, this%ndiv(n)
        ipos = iachk(n) + idiv - 1
        if (nboundchk(ipos) == 0) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
            'No data specified for reach', n, 'diversion', idiv
          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', idiv, &
            'specified', nboundchk(ipos), 'times'
          call store_error(errmsg)
        end if
      end do
    end do
    deallocate (iachk)
    deallocate (nboundchk)
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- release input context diversions memory; DIVFLOW excluded
    call memorystore_release('DIVERSIONS_IFNO', this%input_mempath)
    call memorystore_release('DIVERSIONS_IDV', this%input_mempath)
    call memorystore_release('ICONR', this%input_mempath)
    call memorystore_release('CPRIOR', this%input_mempath)

Here is the call graph for this function:

sfr_source_initialstages()

subroutine sfrmodule::sfr_source_initialstages

(

class(sfrtype), intent(inout)

this

)

Parameters

[in,out]

this

SfrType object

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

    ! -- modules
    use memorymanagermodule, only: mem_setptr, get_isize
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B), dimension(:), pointer, contiguous :: is_ifno => null()
    real(DP), dimension(:), pointer, contiguous :: is_stage => null()
    integer(I4B) :: i, n
    integer(I4B) :: isize
    integer(I4B), allocatable :: nboundchk(:)
    real(DP) :: rval
    !
    ! -- check if INITIALSTAGES block was loaded (optional block)
    call get_isize('INITSTAGE_IFNO', this%input_mempath, isize)
    if (isize < 1) then
      ! -- set default initial stage to reach top if STORAGE active
      if (this%istorage == 1) then
        do n = 1, this%maxbound
          this%stage(n) = this%strtop(n)
        end do
      end if
      return
    end if
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' INITIALSTAGES'
    !
    call mem_setptr(is_ifno, 'INITSTAGE_IFNO', this%input_mempath)
    call mem_setptr(is_stage, 'INITIALSTAGE', this%input_mempath)
    !
    allocate (nboundchk(this%maxbound))
    nboundchk = 0
    !
    do i = 1, isize
      n = this%validate_ifno(is_ifno(i), this%maxbound, nboundchk, &
                             'reach', 'INITIALSTAGES')
      if (n == 0) cycle
      rval = is_stage(i)
      this%stage(n) = rval
      this%depth(n) = rval - this%strtop(n)
      if (rval < this%strtop(n)) then
        write (errmsg, '(a,g0,a,1x,i0,1x,a,g0,a)') &
          'Initial stage (', rval, ') for reach', n, &
          'is less than the reach top (', this%strtop(n), ').'
        call store_error(errmsg)
      end if
    end do
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' INITIALSTAGES'
    !
    call this%report_ifno_coverage(nboundchk, this%maxbound, 'reach', &
                                   'INITIALSTAGES')
    deallocate (nboundchk)
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- release input context initialstages memory
    call memorystore_release('INITSTAGE_IFNO', this%input_mempath)
    call memorystore_release('INITIALSTAGE', this%input_mempath)

Here is the call graph for this function:

sfr_source_options()

subroutine sfrmodule::sfr_source_options

(

class(sfrtype), intent(inout)

this

)

Parameters

[in,out]

this

SfrType object

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

    ! -- modules
    use memorymanagermodule, only: get_isize
    use memorymanagerextmodule, only: mem_set_value
    use inputoutputmodule, only: getunit, assign_iounit, openfile
    use openspecmodule, only: access, form
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    character(len=LINELENGTH) :: fname
    logical(LGP) :: found
    ! -- formats
    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 :: fmtlengthconv = &
      &"(4x, 'LENGTH CONVERSION VALUE (',g0,') SPECIFIED.')"
    character(len=*), parameter :: fmttimeconv = &
      &"(4x, 'TIME CONVERSION VALUE (',g0,') SPECIFIED.')"
    character(len=*), parameter :: fmtstoweight = &
      &"(4x, 'KINEMATIC STORAGE WEIGHT (',g0,') SPECIFIED.')"
    character(len=*), parameter :: fmtsfrbin = &
      "(4x, 'SFR ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, /4x, &
      &'OPENED ON UNIT: ', I0)"
    integer(I4B) :: isize
    !
    ! -- source base BndExtType options (NAUX, BOUNDNAMES, PRINT_INPUT, etc.)
    call this%BndExtType%source_options()
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' OPTIONS'
    !
    ! -- STORAGE
    call mem_set_value(this%istorage, 'STORAGE', this%input_mempath, found)
    if (found) then
      write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
        ' REACH STORAGE IS ACTIVE.'
    end if
    !
    ! -- ATS_COURANT
    call mem_set_value(this%ats_courant, 'ATS_COURANT', &
                       this%input_mempath, found)
    if (found) then
      if (this%ats_courant <= dzero) then
        write (errmsg, '(a,g0,a)') &
          "ATS_COURANT SPECIFIED TO BE '", this%ats_courant, &
          "' BUT MUST BE GREATER THAN ZERO"
        call store_error(errmsg)
        call store_error_filename(this%input_fname)
      else
        write (this%iout, '(4x,a,1pg15.6)') &
          'TARGET COURANT NUMBER FOR ADAPTIVE TIME STEPS: ', &
          this%ats_courant
      end if
    end if
    !
    ! -- PRINT_STAGE
    call mem_set_value(this%iprhed, 'PRINT_STAGE', this%input_mempath, found)
    if (found) then
      write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
        ' STAGES WILL BE PRINTED TO LISTING FILE.'
    end if
    !
    ! -- STAGEFILE
    call mem_set_value(fname, 'STAGEFILE', this%input_mempath, found)
    if (found) then
      this%istageout = getunit()
      call openfile(this%istageout, this%iout, fname, 'DATA(BINARY)', &
                    form, access, 'REPLACE', mnormal)
      write (this%iout, fmtsfrbin) 'STAGE', trim(fname), this%istageout
    end if
    !
    ! -- BUDGETFILE
    call mem_set_value(fname, 'BUDGETFILE', this%input_mempath, found)
    if (found) then
      call assign_iounit(this%ibudgetout, this%inunit, 'BUDGET fileout')
      call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
                    form, access, 'REPLACE', mnormal)
      write (this%iout, fmtsfrbin) 'BUDGET', trim(fname), this%ibudgetout
    end if
    !
    ! -- BUDGETCSVFILE
    call mem_set_value(fname, 'BUDGETCSVFILE', this%input_mempath, found)
    if (found) then
      call assign_iounit(this%ibudcsv, this%inunit, 'BUDGETCSV fileout')
      call openfile(this%ibudcsv, this%iout, fname, 'CSV', &
                    filstat_opt='REPLACE')
      write (this%iout, fmtsfrbin) 'BUDGET CSV', trim(fname), this%ibudcsv
    end if
    !
    ! -- PACKAGE_CONVERGENCE_FILENAME (mf6internal: pkgconvfname)
    call mem_set_value(fname, 'PKGCONVFNAME', this%input_mempath, found)
    if (found) then
      this%ipakcsv = getunit()
      call openfile(this%ipakcsv, this%iout, fname, 'CSV', &
                    filstat_opt='REPLACE', mode_opt=mnormal)
      write (this%iout, fmtsfrbin) 'PACKAGE_CONVERGENCE', trim(fname), &
        this%ipakcsv
    end if
    !
    ! -- MOVER
    call mem_set_value(this%imover, 'MOVER', this%input_mempath, found)
    if (found) then
      write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
    end if
    !
    ! -- MAXIMUM_PICARD_ITERATIONS (mf6internal: maxsfrpicard)
    call mem_set_value(this%maxsfrpicard, 'MAXSFRPICARD', &
                       this%input_mempath, found)
    if (found) write (this%iout, fmtpicard) this%maxsfrpicard
    !
    ! -- MAXIMUM_ITERATIONS (mf6internal: maxsfrit)
    call mem_set_value(this%maxsfrit, 'MAXSFRIT', &
                       this%input_mempath, found)
    if (found) write (this%iout, fmtiter) this%maxsfrit
    !
    ! -- MAXIMUM_DEPTH_CHANGE (mf6internal: dmaxchg)
    call mem_set_value(this%dmaxchg, 'DMAXCHG', &
                       this%input_mempath, found)
    if (found) then
      this%deps = dp999 * this%dmaxchg
      write (this%iout, fmtdmaxchg) this%dmaxchg
    end if
    !
    ! -- LENGTH_CONVERSION (mf6internal: lengthconv)
    call mem_set_value(this%lengthconv, 'LENGTHCONV', &
                       this%input_mempath, found)
    if (found) write (this%iout, fmtlengthconv) this%lengthconv
    !
    ! -- TIME_CONVERSION
    call mem_set_value(this%timeconv, 'TIME_CONVERSION', &
                       this%input_mempath, found)
    if (found) write (this%iout, fmttimeconv) this%timeconv
    !
    ! -- UNIT_CONVERSION (deprecated)
    call mem_set_value(this%unitconv, 'UNIT_CONVERSION', &
                       this%input_mempath, found)
    if (found) then
      write (warnmsg, '(a)') 'SETTING UNIT_CONVERSION DIRECTLY'
      call deprecation_warning('OPTIONS', 'UNIT_CONVERSION', '6.4.2', &
                               warnmsg, this%inunit)
    end if
    !
    ! -- DEV_STORAGE_WEIGHT (mf6internal: storage_weight)
    call mem_set_value(this%storage_weight, 'STORAGE_WEIGHT', &
                       this%input_mempath, found)
    if (found) then
      if (this%storage_weight < dhalf .or. this%storage_weight > done) then
        write (errmsg, '(a,g0,a)') &
          "STORAGE_WEIGHT SPECIFIED TO BE '", this%storage_weight, &
          "' BUT CANNOT BE LESS THAN 0.5 OR GREATER THAN 1.0"
        call store_error(errmsg)
        call store_error_filename(this%input_fname)
      else
        write (this%iout, fmtstoweight) this%storage_weight
      end if
    end if
    !
    ! -- DEV_NO_CHECK (mf6internal: icheck)
    call get_isize('ICHECK', this%input_mempath, isize)
    if (isize > 0) then
      this%icheck = 0
      write (this%iout, '(4x,a)') &
        'REACH GEOMETRY AND PARAMETER CHECKS WILL NOT BE MADE'
    end if
    !
    ! -- DEV_NO_FINAL_CHECK (mf6internal: iconvchk)
    call get_isize('ICONVCHK', this%input_mempath, isize)
    if (isize > 0) then
      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'
    end if
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' OPTIONS'

Here is the call graph for this function:

sfr_source_packagedata()

subroutine sfrmodule::sfr_source_packagedata

(

class(sfrtype), intent(inout)

this

)

Parameters

[in,out]

this

SfrType object

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

    ! -- modules
    use memorymanagermodule, only: mem_setptr
    use characterstringmodule, only: characterstringtype
    ! -- dummy
    class(SfrType), intent(inout) :: this !< SfrType object
    ! -- local
    integer(I4B) :: i, n, nconzero
    integer(I4B) :: ipos
    character(len=10) :: cnum
    character(len=LENBOUNDNAME) :: bndName
    character(len=LINELENGTH) :: cellid_str
    integer(I4B), allocatable :: nboundchk(:)
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    integer(I4B), dimension(:, :), pointer, contiguous :: cellid => null()
    real(DP), dimension(:), pointer, contiguous :: rlen => null()
    real(DP), dimension(:), pointer, contiguous :: rwid => null()
    real(DP), dimension(:), pointer, contiguous :: rgrd => null()
    real(DP), dimension(:), pointer, contiguous :: rtp => null()
    real(DP), dimension(:), pointer, contiguous :: rbth => null()
    integer(I4B), dimension(:), pointer, contiguous :: ncon => null()
    integer(I4B), dimension(:), pointer, contiguous :: ndv => null()
    type(CharacterStringType), dimension(:), pointer, contiguous :: &
      boundname => null()
    integer(I4B), dimension(:, :), pointer, contiguous :: tslink => null()
    type(CharacterStringType), dimension(:), pointer, contiguous :: &
      tslinkvars => null()
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' PACKAGEDATA'
    !
    ! -- get input context PACKAGEDATA column arrays
    call mem_setptr(ifno, 'PACKAGEDATA_IFNO', this%input_mempath)
    call mem_setptr(cellid, 'CELLID', this%input_mempath)
    call mem_setptr(rlen, 'RLEN', this%input_mempath)
    call mem_setptr(rwid, 'RWID', this%input_mempath)
    call mem_setptr(rgrd, 'RGRD', this%input_mempath)
    call mem_setptr(rtp, 'RTP', this%input_mempath)
    call mem_setptr(rbth, 'RBTH', this%input_mempath)
    call mem_setptr(this%input%rhk, 'RHK', this%input_mempath)
    call mem_setptr(this%input%man, 'MAN', this%input_mempath)
    call mem_setptr(ncon, 'NCON', this%input_mempath)
    call mem_setptr(this%input%ustrf, 'USTRF', this%input_mempath)
    call mem_setptr(tslink, 'TSLINK_PD', this%input_mempath)
    call mem_setptr(tslinkvars, 'TSLINKVARS_PD', this%input_mempath)
    do n = 1, size(tslinkvars)
      if (tslinkvars(n) == 'RHK') this%input%rhk_tslink => tslink(n, :)
      if (tslinkvars(n) == 'MAN') this%input%man_tslink => tslink(n, :)
      if (tslinkvars(n) == 'USTRF') this%input%ustrf_tslink => tslink(n, :)
    end do
    if (.not. (associated(this%input%rhk_tslink) .and. &
               associated(this%input%man_tslink) .and. &
               associated(this%input%ustrf_tslink))) then
      call store_error('Time series links for RHK, MAN and USTRF not found &
                       &in input.')
      call store_error_filename(this%input_fname)
    end if
    call mem_setptr(ndv, 'NDV', this%input_mempath)
    if (this%inamedbound /= 0) then
      call mem_setptr(boundname, 'BOUNDNAME', this%input_mempath)
    end if
    !
    ! -- allocate check array
    allocate (nboundchk(this%maxbound))
    nboundchk = 0
    nconzero = 0
    !
    ! -- populate reach data from input context arrays
    do i = 1, this%maxbound
      n = ifno(i)
      !
      n = this%validate_ifno(n, this%maxbound, nboundchk, &
                             'rno', 'PACKAGEDATA')
      if (n == 0) cycle
      !
      ! -- cellid → igwfnode: convert integer cellid array to string
      ! -- all-zero or negative cellid means unconnected reach (NONE keyword or 0 0 0)
      if (all(cellid(1:this%dis%ndim, i) <= 0)) then
        this%igwfnode(n) = 0
      else
        write (cellid_str, '(10(i0,1x))') cellid(1:this%dis%ndim, i)
        this%igwfnode(n) = this%dis%noder_from_cellid(cellid_str, &
                                                      this%inunit, this%iout, &
                                                      flag_string=.false., &
                                                      allow_zero=.true.)
      end if
      this%igwftopnode(n) = this%igwfnode(n)
      if (this%igwfnode(n) < 1) this%ianynone = this%ianynone + 1
      !
      ! -- scalar reach properties
      this%length(n) = rlen(i)
      this%width(n) = rwid(i)
      this%slope(n) = rgrd(i)
      this%strtop(n) = rtp(i)
      this%bthick(n) = rbth(i)
      !
      ! -- ncon and nconn
      this%nconnreach(n) = ncon(i)
      this%nconn = this%nconn + ncon(i)
      if (ncon(i) < 0) then
        write (errmsg, '(a,1x,i0,1x,a,i0,a)') &
          'NCON for reach', n, &
          'must be greater than or equal to 0 (', ncon(i), ').'
        call store_error(errmsg)
      else if (ncon(i) == 0) then
        nconzero = nconzero + 1
      end if
      !
      ! -- ndv (diversions)
      this%ndiv(n) = ndv(i)
      if (ndv(i) > 0) this%idiversions = 1
      !
      ! -- bound name
      write (cnum, '(i10.10)') n
      bndname = 'Reach'//cnum
      if (this%inamedbound /= 0) then
        bndname = boundname(i)
      end if
      this%sfrname(n) = bndname
      !
      ! -- hk, roughness (Manning), upstream fraction: copied from
      !    input context; sfr_ad re-syncs TS-linked values each timestep
      this%hk(n) = this%input%rhk(i)
      this%rough(n) = this%input%man(i)
      this%ustrf(n) = this%input%ustrf(i)
      !
      ! -- initialize sstage to reach top
      this%sstage(n) = this%strtop(n)
    end do
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' PACKAGEDATA'
    !
    ! -- check every reach is specified exactly once
    call this%report_ifno_coverage(nboundchk, this%maxbound, 'reach', &
                                   'PACKAGEDATA')
    deallocate (nboundchk)
    !
    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
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- initialize default 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
    !
    ! -- release input context packagedata memory; PACKAGEDATA_IFNO excluded
    !    (allocate_featureauxvar); RHK/MAN/USTRF excluded
    call memorystore_release('CELLID', this%input_mempath)
    call memorystore_release('RLEN', this%input_mempath)
    call memorystore_release('RWID', this%input_mempath)
    call memorystore_release('RGRD', this%input_mempath)
    call memorystore_release('RTP', this%input_mempath)
    call memorystore_release('RBTH', this%input_mempath)
    call memorystore_release('NCON', this%input_mempath)
    call memorystore_release('NDV', this%input_mempath)
    call memorystore_release('BOUNDNAME', this%input_mempath)

Here is the call graph for this function:

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 3621 of file gwf-sfr.f90.

    ! -- dummy
    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
    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

Variable Documentation

ftype

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

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

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

text

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

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

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