gwf-sfr.f90

Go to the documentation of this file.

!> @brief This module contains the SFR package methods
!!
!! This module contains the overridden methods for the streamflow routing (SFR)
!! package. Most of the methods in the base Boundary Package are overridden.
!!
!<
module sfrmodule
  !
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: linelength, lenboundname, &
                             maxadpit, &
                             dzero, dprec, dem30, dem6, dem5, dem4, dem2, &
                             donethird, dhalf, dp6, dtwothirds, dp7, &
                             dp9, dp99, dp999, &
                             done, d1p1, dfivethirds, dtwo, dpi, deight, &
                             dhundred, dep20, &
                             namedboundflag, lenboundname, lenftype, &
                             lenpackagename, lenpakloc, maxcharlen, &
                             lenbudtxt, &
                             dhnoflo, dhdry, dnodata, &
                             tableft, tabcenter, tabright, &
                             mnormal
  use smoothingmodule, only: squadraticsaturation, sqsaturation, &
                             squadraticsaturationderivative, &
                             sqsaturationderivative, &
                             scubicsaturation, schsmooth
  use characterstringmodule, only: characterstringtype
  use bndmodule, only: bndtype
  use bndextmodule, only: bndexttype
  use memorymanagerextmodule, only: memorystore_release
 
  use budgetobjectmodule, only: budgetobjecttype, budgetobject_cr
  use tablemodule, only: tabletype, table_cr
  use observemodule, only: observetype
  use inputoutputmodule, only: extract_idnum_or_bndname, upcase
  use basedismodule, only: disbasetype
  use simmodule, only: count_errors, store_error, store_error_unit, &
                       store_error_filename, store_warning, deprecation_warning
  use simvariablesmodule, only: errmsg, warnmsg
  use gwfsfrcrosssectionutilsmodule, only: get_saturated_topwidth, &
                                           get_wetted_topwidth, &
                                           get_wetted_perimeter, &
                                           get_cross_section_area, &
                                           get_mannings_section
  use dag_module, only: dag
  use matrixbasemodule
  !
  implicit none
  !
  character(len=LENFTYPE) :: ftype = 'SFR' !< package ftype string
  character(len=LENPACKAGENAME) :: text = '             SFR' !< package budget string
  !
  private
  public :: sfr_create
  public :: sfrtype
 
  !> @brief Input context array pointers for SFR.
  !!
  !! PACKAGEDATA pointers are set in sfr_source_packagedata; PERIOD
  !! pointers are set once in sfr_ar.
  !<
  type :: sfrinputtype
    ! -- packagedata (TS-linked; set in sfr_source_packagedata)
    real(dp), dimension(:), pointer, contiguous :: rhk => null()
    real(dp), dimension(:), pointer, contiguous :: man => null()
    real(dp), dimension(:), pointer, contiguous :: ustrf => null()
    ! -- period
    integer(I4B), pointer :: nbound => null()
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: status => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: setting => null()
    integer(I4B), dimension(:), pointer, contiguous :: idv => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: tab6_filename => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: auxname => null()
    real(dp), dimension(:), pointer, contiguous :: auxval => null()
  contains
    procedure :: init => sfr_input_init
    procedure :: destroy => sfr_input_destroy
  end type sfrinputtype
 
  type, extends(bndexttype) :: sfrtype
    ! -- scalars
    ! -- for budgets
    ! -- characters
    character(len=16), dimension(:), pointer, contiguous :: csfrbudget => null() !< advanced package budget names
    character(len=16), dimension(:), pointer, contiguous :: cauxcbc => null() !< aux names
    character(len=LENBOUNDNAME), dimension(:), pointer, &
      contiguous :: sfrname => null() !< internal SFR reach name
    ! -- integers
    integer(I4B), pointer :: istorage => null() !< flag for using kinematic wave approximation
    integer(I4B), pointer :: iprhed => null() !< flag for printing stages to listing file
    integer(I4B), pointer :: istageout => null() !< flag and unit number for binary stage output
    integer(I4B), pointer :: ibudgetout => null() !< flag and unit number for binary sfr budget output
    integer(I4B), pointer :: ibudcsv => null() !< unit number for csv budget output file
    integer(I4B), pointer :: ipakcsv => null() !< flag and unit number for package convergence information
    integer(I4B), pointer :: idiversions => null() !< flag indicating if there are any diversions
    integer(I4B), pointer :: nconn => null() !< number of reach connections
    integer(I4B), pointer :: maxsfrpicard => null() !< maximum number of Picard iteration calls to SFR solve
    integer(I4B), pointer :: maxsfrit => null() !< maximum number of iterations in SFR solve
    integer(I4B), pointer :: bditems => null() !< number of SFR budget items
    integer(I4B), pointer :: cbcauxitems => null() !< number of aux items in cell-by-cell budget file
    integer(I4B), pointer :: icheck => null() !< flag indicating if input should be checked (default is yes)
    integer(I4B), pointer :: iconvchk => null() !< flag indicating of final convergence run is executed
    integer(I4B), pointer :: gwfiss => null() !< groundwater model steady-state flag
    integer(I4B), pointer :: ianynone => null() !< number of reaches with 'none' connection
    ! -- double precision
    real(dp), pointer :: unitconv => null() !< unit conversion factor (SI to model units)
    real(dp), pointer :: lengthconv => null() !< length conversion factor (SI to model units)
    real(dp), pointer :: timeconv => null() !< time conversion factor (SI to model units)
    real(dp), pointer :: dmaxchg => null() !< maximum depth change allowed
    real(dp), pointer :: deps => null() !< perturbation value
    real(dp), pointer :: storage_weight => null() !< time weighting factor for kinematic wave approximation
    ! -- integer vectors
    integer(I4B), dimension(:), pointer, contiguous :: isfrorder => null() !< sfr reach order determined from DAG of upstream reaches
    integer(I4B), dimension(:), pointer, contiguous :: ia => null() !< CRS row pointer for SFR reaches
    integer(I4B), dimension(:), pointer, contiguous :: ja => null() !< CRS column pointers for SFR reach connections
    ! -- double precision output vectors
    real(dp), dimension(:), pointer, contiguous :: qoutflow => null() !< reach downstream flow
    real(dp), dimension(:), pointer, contiguous :: qextoutflow => null() !< reach discharge to external boundary
    real(dp), dimension(:), pointer, contiguous :: qauxcbc => null() !< aux value
    real(dp), dimension(:), pointer, contiguous :: dbuff => null() !< temporary vector
    !
    ! -- sfr budget object
    type(budgetobjecttype), pointer :: budobj => null() !< SFR budget object
    !
    ! -- sfr table objects
    type(tabletype), pointer :: stagetab => null() !< reach stage table written to the listing file
    type(tabletype), pointer :: couranttab => null() !< Courant number table written to the listing file
    type(tabletype), pointer :: pakcsvtab => null() !< SFR package convergence table
    !
    ! -- sfr reach data
    integer(I4B), dimension(:), pointer, contiguous :: iboundpak => null() !< ibound array for SFR reaches that defines active, inactive, and constant reaches
    integer(I4B), dimension(:), pointer, contiguous :: igwfnode => null() !< groundwater node connected to SFR reaches
    integer(I4B), dimension(:), pointer, contiguous :: igwftopnode => null() !< highest active groundwater node under SFR reaches
    real(dp), dimension(:), pointer, contiguous :: length => null() !< reach length
    real(dp), dimension(:), pointer, contiguous :: width => null() !< reach width
    real(dp), dimension(:), pointer, contiguous :: strtop => null() !< reach bed top elevation
    real(dp), dimension(:), pointer, contiguous :: bthick => null() !< reach bed thickness
    real(dp), dimension(:), pointer, contiguous :: hk => null() !< vertical hydraulic conductivity of reach bed sediments
    real(dp), dimension(:), pointer, contiguous :: slope => null() !< reach slope
    integer(I4B), dimension(:), pointer, contiguous :: nconnreach => null() !< number of connections for each reach
    real(dp), dimension(:), pointer, contiguous :: ustrf => null() !< upstream flow fraction for upstream connections
    logical(LGP), dimension(:), pointer, contiguous :: bedk_set => null() !< reach has an active PERIOD BEDK override
    logical(LGP), dimension(:), pointer, contiguous :: manning_set => null() !< reach has an active PERIOD MANNING override
    logical(LGP), dimension(:), pointer, contiguous :: ustrf_set => null() !< reach has an active PERIOD UPSTREAM_FRACTION override
    real(dp), dimension(:), pointer, contiguous :: ftotnd => null() !< total fraction of connected reaches that are not diversions
    integer(I4B), dimension(:), pointer, contiguous :: ndiv => null() !< number of diversions for each reach
    real(dp), dimension(:), pointer, contiguous :: usflow => null() !< upstream reach flow
    real(dp), dimension(:), pointer, contiguous :: dsflow => null() !< downstream reach flow
    real(dp), dimension(:), pointer, contiguous :: dsflowold => null() !< downstream reach flow for previous time step
    real(dp), dimension(:), pointer, contiguous :: usinflow => null() !< upstream reach flow for previous time step
    real(dp), dimension(:), pointer, contiguous :: usinflowold => null() !< upstream reach flow for previous time step
    real(dp), pointer :: ats_courant => null() !< target Courant number for ATS time step submission
    integer(I4B), dimension(:), pointer, contiguous :: itvd_upstream => null() !< upstream reach index for TVD limiter (0 = headwater or confluence)
    real(dp), dimension(:), pointer, contiguous :: crmin => null() !< simulation-wide minimum Courant number per reach
    real(dp), dimension(:), pointer, contiguous :: crmax => null() !< simulation-wide maximum Courant number per reach
    real(dp), dimension(:), pointer, contiguous :: crsum => null() !< simulation-wide sum of Courant numbers per reach
    integer(I4B), dimension(:), pointer, contiguous :: crcnt => null() !< simulation-wide count of non-zero Courant evaluations per reach
    real(dp), dimension(:), pointer, contiguous :: depth => null() !< reach depth
    real(dp), dimension(:), pointer, contiguous :: stage => null() !< reach stage
    real(dp), dimension(:), pointer, contiguous :: stageold => null() !< reach stage for last timestep
    real(dp), dimension(:), pointer, contiguous :: gwflow => null() !< flow from groundwater to reach
    real(dp), dimension(:), pointer, contiguous :: simevap => null() !< simulated reach evaporation
    real(dp), dimension(:), pointer, contiguous :: simrunoff => null() !< simulated reach runoff
    real(dp), dimension(:), pointer, contiguous :: stage0 => null() !< previous reach stage iterate
    real(dp), dimension(:), pointer, contiguous :: usflow0 => null() !< previous upstream reach flow iterate
    real(dp), dimension(:), pointer, contiguous :: storage => null() !< previous upstream reach flow iterate
    ! -- cross-section data
    integer(I4B), pointer :: ncrossptstot => null() !< total number of cross-section points
    integer(I4B), dimension(:), pointer, contiguous :: ncrosspts => null() !< number of cross-section points for each reach
    integer(I4B), dimension(:), pointer, contiguous :: iacross => null() !< pointers to cross-section data for each reach
    real(dp), dimension(:), pointer, contiguous :: station => null() !< cross-section station (x-position) data
    real(dp), dimension(:), pointer, contiguous :: xsheight => null() !< cross-section height data
    real(dp), dimension(:), pointer, contiguous :: xsrough => null() !< cross-section roughness data
    ! -- connection data
    integer(I4B), dimension(:), pointer, contiguous :: idir => null() !< reach connection direction
    integer(I4B), dimension(:), pointer, contiguous :: idiv => null() !< reach connection diversion number
    real(dp), dimension(:), pointer, contiguous :: qconn => null() !< reach connection flow
    ! -- boundary data
    real(dp), dimension(:), pointer, contiguous :: rough => null() !< reach Manning's roughness coefficient (SI units)
    real(dp), dimension(:), pointer, contiguous :: rain => null() !< reach rainfall
    real(dp), dimension(:), pointer, contiguous :: evap => null() !< reach potential evaporation
    real(dp), dimension(:), pointer, contiguous :: inflow => null() !< reach upstream inflow
    real(dp), dimension(:), pointer, contiguous :: runoff => null() !< reach maximum runoff
    real(dp), dimension(:), pointer, contiguous :: sstage => null() !< reach specified stage
    ! -- diversion data
    integer(I4B), dimension(:), pointer, contiguous :: iadiv => null() !< row pointer for reach diversions
    integer(I4B), dimension(:), pointer, contiguous :: divreach => null() !< diversion reach
    character(len=10), dimension(:), pointer, contiguous :: divcprior => null() !< diversion rule
    real(dp), dimension(:), pointer, contiguous :: divflow => null() !< specified diversion flow value
    real(dp), dimension(:), pointer, contiguous :: divq => null() !< simulated diversion flow
    !
    ! -- density variables
    integer(I4B), pointer :: idense !< flag indicating if density corrections are active
    real(dp), dimension(:, :), pointer, contiguous :: denseterms => null() !< density terms
    !
    ! -- viscosity variables
    real(dp), dimension(:, :), pointer, contiguous :: viscratios => null() !< viscosity ratios (1: sfr vsc ratio; 2: gwf vsc ratio)
    !
    ! -- input context data pointers (PACKAGEDATA + PERIOD arrays)
    type(sfrinputtype) :: input
    !
    ! -- type bound procedures
  contains
    procedure :: sfr_allocate_scalars
    procedure :: sfr_allocate_arrays
    procedure :: source_options => sfr_source_options
    procedure :: source_dimensions => sfr_source_dimensions
 
    ! procedure :: set_pointers => sfr_set_pointers
    procedure :: bnd_ar => sfr_ar
    procedure :: bnd_rp => sfr_rp
    procedure :: bnd_ad => sfr_ad
    procedure :: bnd_cf => sfr_cf
    procedure :: bnd_fc => sfr_fc
    procedure :: bnd_fn => sfr_fn
    procedure :: bnd_cc => sfr_cc
    procedure :: bnd_cq => sfr_cq
    procedure :: bnd_ot_package_flows => sfr_ot_package_flows
    procedure :: bnd_ot_dv => sfr_ot_dv
    procedure :: bnd_ot_bdsummary => sfr_ot_bdsummary
    procedure :: bnd_dt => sfr_dt
    procedure :: bnd_fp => sfr_fp
    procedure :: bnd_da => sfr_da
    procedure :: define_listlabel
    ! -- methods for observations
    procedure, public :: bnd_obs_supported => sfr_obs_supported
    procedure, public :: bnd_df_obs => sfr_df_obs
    procedure, public :: bnd_rp_obs => sfr_rp_obs
    procedure, public :: bnd_bd_obs => sfr_bd_obs
    ! -- private procedures
 
    procedure, private :: sfr_solve
    procedure, private :: sfr_calc_constant
    procedure, private :: sfr_calc_transient
    procedure, private :: sfr_calc_tvd
    procedure, private :: sfr_calc_celerity
    procedure, private :: sfr_calc_steady
    procedure, private :: sfr_precompute_tvd
    procedure, private :: sfr_update_flows
    procedure, private :: sfr_adjust_ro_ev
    procedure, private :: sfr_calc_qgwf
    procedure, private :: sfr_gwf_conn
    procedure, private :: sfr_calc_cond
    procedure, private :: sfr_calc_qman
    procedure, private :: sfr_calc_qd
    procedure, private :: sfr_calc_qsource
    procedure, private :: sfr_calc_div
    ! -- geometry
    procedure, private :: calc_area_wet
    procedure, private :: calc_perimeter_wet
    procedure, private :: calc_surface_area
    procedure, private :: calc_surface_area_wet
    procedure, private :: calc_top_width_wet
    ! -- input context sourcing
    procedure, private :: sfr_source_packagedata
    procedure, private :: sfr_source_crosssections
    procedure, private :: sfr_source_connectiondata
    procedure, private :: sfr_source_diversions
    procedure, private :: sfr_source_initialstages
    ! -- calculations
    procedure, private :: sfr_calc_reach_depth
    procedure, private :: sfr_calc_xs_depth
    ! -- error checking
    procedure, private :: sfr_check_conversion
    procedure, private :: sfr_check_storage_weight
    procedure, private :: sfr_check_reaches
    procedure, private :: sfr_check_connections
    procedure, private :: sfr_check_diversions
    procedure, private :: sfr_check_initialstages
    procedure, private :: sfr_check_ustrf
    procedure, private :: sfr_set_period_value
    ! -- budget
    procedure, private :: sfr_setup_budobj
    procedure, private :: sfr_fill_budobj
    ! -- table
    procedure, private :: sfr_setup_tableobj
    ! -- density
    procedure :: sfr_activate_density
    procedure, private :: sfr_calculate_density_exchange
    ! -- viscosity
    procedure :: sfr_activate_viscosity
  end type sfrtype
 
  interface
    module subroutine sfr_calc_steady(this, n, d1, hgwf, &
                                      qu, qi, qfrommvr, qr, qe, qro, &
                                      qgwf, qd)
      class(sfrtype) :: this !< SfrType object
      integer(I4B), intent(in) :: n !< reach number
      real(dp), intent(inout) :: d1 !< current reach depth estimate
      real(dp), intent(in) :: hgwf !< head in gw cell
      real(dp), intent(in) :: qu !< reach upstream flow
      real(dp), intent(in) :: qi !< reach specified inflow
      real(dp), intent(in) :: qfrommvr !< reach flow from mover
      real(dp), intent(in) :: qr !< reach rainfall
      real(dp), intent(in) :: qe !< reach evaporation
      real(dp), intent(in) :: qro !< reach runoff flow
      real(dp), intent(inout) :: qgwf !< reach-aquifer exchange
      real(dp), intent(inout) :: qd !< reach outflow
    end subroutine
  end interface
 
  interface
    module subroutine sfr_calc_transient(this, n, d1, hgwf, &
                                         qu, qi, qfrommvr, qr, qe, qro, &
                                         qgwf, qd)
      class(sfrtype) :: this !< SfrType object
      integer(I4B), intent(in) :: n !< reach number
      real(dp), intent(inout) :: d1 !< current reach depth estimate
      real(dp), intent(in) :: hgwf !< head in gw cell
      real(dp), intent(in) :: qu !< reach upstream flow
      real(dp), intent(in) :: qi !< reach specified inflow
      real(dp), intent(in) :: qfrommvr !< reach flow from mover
      real(dp), intent(in) :: qr !< reach rainfall
      real(dp), intent(in) :: qe !< reach evaporation
      real(dp), intent(in) :: qro !< reach runoff flow
      real(dp), intent(inout) :: qgwf !< reach-aquifer exchange
      real(dp), intent(inout) :: qd !< reach outflow
    end subroutine
  end interface
 
  interface
    module subroutine sfr_calc_tvd(this, n, d1, hgwf, &
                                   qu, qi, qfrommvr, qr, qe, qro, &
                                   qgwf, qd)
      class(sfrtype) :: this !< SfrType object
      integer(I4B), intent(in) :: n !< reach number
      real(dp), intent(inout) :: d1 !< current reach depth estimate
      real(dp), intent(in) :: hgwf !< head in gw cell
      real(dp), intent(in) :: qu !< reach upstream flow
      real(dp), intent(in) :: qi !< reach specified inflow
      real(dp), intent(in) :: qfrommvr !< reach flow from mover
      real(dp), intent(in) :: qr !< reach rainfall
      real(dp), intent(in) :: qe !< reach evaporation
      real(dp), intent(in) :: qro !< reach runoff flow
      real(dp), intent(inout) :: qgwf !< reach-aquifer exchange
      real(dp), intent(inout) :: qd !< reach outflow
    end subroutine
  end interface
 
  interface
    module subroutine sfr_calc_constant(this, n, d1, hgwf, qgwf, qd)
      class(sfrtype) :: this !< SfrType object
      integer(I4B), intent(in) :: n !< reach number
      real(dp), intent(inout) :: d1 !< current reach depth estimate
      real(dp), intent(in) :: hgwf !< head in gw cell
      real(dp), intent(inout) :: qgwf !< reach-aquifer exchange
      real(dp), intent(inout) :: qd !< reach outflow
    end subroutine
  end interface
 
contains
 
  !> @ brief Create a new package object
  !!
  !!  Create a new SFR Package object
  !<
  subroutine sfr_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
                        mempath)
    ! -- 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')
  end subroutine sfr_create
 
  !> @ brief Allocate scalars
  !!
  !! Allocate and initialize scalars for the SFR package. The base model
  !! allocate scalars method is also called.
  !<
  subroutine sfr_allocate_scalars(this)
    ! -- 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
  end subroutine sfr_allocate_scalars
 
  !> @ brief Allocate arrays
  !!
  !! Allocate and initialize array for the SFR package.
  !<
  subroutine sfr_allocate_arrays(this)
    ! -- 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)
  end subroutine sfr_allocate_arrays
 
  !> @brief Source OPTIONS from input context
  !<
  subroutine sfr_source_options(this)
    ! -- 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'
  end subroutine sfr_source_options
 
  !> @brief Source DIMENSIONS block from input context and load static data
  !<
  subroutine sfr_source_dimensions(this)
    ! -- 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()
    !
    ! -- 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()
  end subroutine sfr_source_dimensions
 
  !> @brief Source PACKAGEDATA block from input context
  !<
  subroutine sfr_source_packagedata(this)
    ! -- 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()
    !
    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(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)
  end subroutine sfr_source_packagedata
 
  !> @brief Source CROSSSECTIONS block from input context
  !<
  subroutine sfr_source_crosssections(this)
    ! -- 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)
  end subroutine sfr_source_crosssections
 
  !> @brief Source CONNECTIONDATA block from input context
  !<
  subroutine sfr_source_connectiondata(this)
    ! -- 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)
  end subroutine sfr_source_connectiondata
 
  !> @brief Source DIVERSIONS block from input context
  !<
  subroutine sfr_source_diversions(this)
    ! -- 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)
  end subroutine sfr_source_diversions
 
  !> @brief Source INITIALSTAGES block from input context
  !<
  subroutine sfr_source_initialstages(this)
    ! -- 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)
  end subroutine sfr_source_initialstages
 
  !> @ brief Allocate and read method for package
  !!
  !!  Method to read and prepare period data for the SFR package.
  !<
  subroutine sfr_ar(this)
    ! -- 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
  end subroutine sfr_ar
 
  !> @ brief Read and prepare period data for package
  !!
  !!  Method to read and prepare period data for the SFR package.
  !<
  subroutine sfr_rp(this)
    ! -- 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
  end subroutine sfr_rp
 
  !> @ brief Advance the package
  !!
  !!  Advance data in the SFR package. The method sets advances
  !!  time series, time array series, and observation data.
  !<
  subroutine sfr_ad(this)
    ! -- 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
    if (this%ts_active) then
      do i = 1, this%maxbound
        n = this%pkg_ifno(i)
        if (.not. this%bedk_set(n)) this%hk(n) = this%input%rhk(i)
        if (.not. this%manning_set(n)) this%rough(n) = this%input%man(i)
        if (.not. this%ustrf_set(n)) 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
  end subroutine sfr_ad
 
  !> @ brief Formulate the package hcof and rhs terms.
  !!
  !!  Formulate the hcof and rhs terms for the WEL package that will be
  !!  added to the coefficient matrix and right-hand side vector.
  !<
  subroutine sfr_cf(this)
    ! -- 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
  end subroutine sfr_cf
 
  !> @ brief Copy hcof and rhs terms into solution.
  !!
  !!  Add the hcof and rhs terms for the SFR package to the
  !!  coefficient matrix and right-hand side vector.
  !<
  subroutine sfr_fc(this, rhs, ia, idxglo, matrix_sln)
    ! -- 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
  end subroutine sfr_fc
 
  !> @ brief Add Newton-Raphson terms for package into solution.
  !!
  !!  Calculate and add the Newton-Raphson terms for the SFR package to the
  !!  coefficient matrix and right-hand side vector.
  !<
  subroutine sfr_fn(this, rhs, ia, idxglo, matrix_sln)
    ! -- 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
  end subroutine sfr_fn
 
  !> @ brief Convergence check for package.
  !!
  !!  Perform additional convergence checks on the flow between the SFR package
  !!  and the model it is attached to.
  !<
  subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
    ! -- 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
  end subroutine sfr_cc
 
  !> @ brief Calculate package flows.
  !!
  !!  Calculate the flow between connected SFR package control volumes.
  !<
  subroutine sfr_cq(this, x, flowja, iadv)
    ! -- 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()
  end subroutine sfr_cq
 
  !> @ brief Output package flow terms.
  !!
  !!  Output SFR package flow terms.
  !<
  subroutine sfr_ot_package_flows(this, icbcfl, ibudfl)
    ! -- 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
  end subroutine sfr_ot_package_flows
 
  !> @ brief Output package dependent-variable terms.
  !!
  !!  Output SFR boundary package dependent-variable terms.
  !<
  subroutine sfr_ot_dv(this, idvsave, idvprint)
    ! -- 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
  end subroutine sfr_ot_dv
 
  !> @ brief Output advanced package budget summary.
  !!
  !!  Output SFR package budget summary.
  !<
  subroutine sfr_ot_bdsummary(this, kstp, kper, iout, ibudfl)
    ! -- 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)
  end subroutine sfr_ot_bdsummary
 
  !> @brief Write the Courant-number summary table to the listing file
  !<
  subroutine sfr_fp(this)
    ! -- 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
  end subroutine sfr_fp
 
  !> @brief Submit the ATS time step for the most Courant-constraining reach
  !<
  subroutine sfr_dt(this)
    ! -- 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
  end subroutine sfr_dt
 
  !> @brief Kinematic-wave celerity from a flow perturbation (0 if dry)
  !<
  subroutine sfr_calc_celerity(this, n, q, celerity)
    ! -- 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
  end subroutine sfr_calc_celerity
 
  !> @ brief Deallocate package memory
  !!
  !!  Deallocate SFR package scalars and arrays.
  !<
  subroutine sfr_da(this)
    ! -- 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
    nullify (this%divflow)
    call mem_deallocate(this%divq)
    call mem_deallocate(this%ndiv)
    !
    ! -- deallocate cross-section data
    call mem_deallocate(this%ncrosspts)
    call mem_deallocate(this%iacross)
    call mem_deallocate(this%station)
    call mem_deallocate(this%xsheight)
    call mem_deallocate(this%xsrough)
    !
    ! -- deallocate budobj
    call this%budobj%budgetobject_da()
    deallocate (this%budobj)
    nullify (this%budobj)
    !
    ! -- deallocate stage table
    if (this%iprhed > 0) then
      call this%stagetab%table_da()
      deallocate (this%stagetab)
      nullify (this%stagetab)
    end if
    !
    ! -- deallocate 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()
  end subroutine sfr_da
 
  !> @brief Pre-compute the single upstream reach index for the TVD limiter
  !!
  !! Index 0 for reaches with zero or multiple upstream connections.
  !<
  subroutine sfr_precompute_tvd(this)
    ! -- 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
  end subroutine sfr_precompute_tvd
 
  !> @ brief Define the list label for the package
  !!
  !!  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.
  !<
  subroutine define_listlabel(this)
    ! -- 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
  end subroutine define_listlabel
 
  !
  ! -- Procedures related to observations (type-bound)
 
  !> @brief Determine if observations are supported.
  !!
  !! Function to determine if observations are supported by the SFR package.
  !! Observations are supported by the SFR package.
  !!
  !! @return  sfr_obs_supported       boolean indicating if observations are supported
  !<
  logical function sfr_obs_supported(this)
    ! -- dummy
    class(sfrtype) :: this !< SfrType object
    !
    ! -- set boolean
    sfr_obs_supported = .true.
  end function sfr_obs_supported
 
  !> @brief Define the observation types available in the package
  !!
  !! Method to define the observation types available in the SFR package.
  !<
  subroutine sfr_df_obs(this)
    ! -- 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
  end subroutine sfr_df_obs
 
  !> @brief Save observations for the package
  !!
  !! Method to save simulated values for the SFR package.
  !<
  subroutine sfr_bd_obs(this)
    ! -- 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
  end subroutine sfr_bd_obs
 
  !> @brief Read and prepare observations for a package
  !!
  !! Method to read and prepare observations for a SFR package.
  !<
  subroutine sfr_rp_obs(this)
    ! -- 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
  end subroutine sfr_rp_obs
 
  !
  ! -- Procedures related to observations (NOT type-bound)
 
  !> @brief Process observation IDs for a package
  !!
  !! Method to process observation ID strings for a SFR package.
  !<
  subroutine sfr_process_obsid(obsrv, dis, inunitobs, iout)
    ! -- 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
  end subroutine sfr_process_obsid
 
  !> @brief Solve reach continuity equation
  !!
  !! Method to solve the continuity equation for a SFR package reach.
  !<
  subroutine sfr_solve(this, n, h, hcof, rhs, update)
    ! -- 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
  end subroutine sfr_solve
 
  !> @brief Update flow terms
  !!
  !! Method to update downstream flow and groundwater leakage terms for
  !! a SFR package reach.
  !<
  subroutine sfr_update_flows(this, n, qd, qgwf)
    ! -- 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
  end subroutine sfr_update_flows
 
  !> @brief Adjust runoff and evaporation
  !!
  !! Method to adjust runoff and evaporation for a SFR package reach
  !! based on the total reach flow.
  !<
  subroutine sfr_adjust_ro_ev(this, qc, qu, qi, qr, qro, qe, qfrommvr)
    ! -- 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
  end subroutine sfr_adjust_ro_ev
 
  !> @brief Calculate downstream flow term
  !!
  !! Method to calculate downstream flow for a SFR package reach.
  !<
  subroutine sfr_calc_qd(this, n, depth, hgwf, qgwf, qd)
    ! -- 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
  end subroutine sfr_calc_qd
 
  !> @brief Calculate sum of sources
  !!
  !! Method to calculate the sum of sources for reach, excluding
  !! reach leakage, for a SFR package reach.
  !<
  subroutine sfr_calc_qsource(this, n, depth, qsrc)
    ! -- 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)
  end subroutine sfr_calc_qsource
 
  !> @brief Calculate streamflow
  !!
  !! Method to calculate the streamflow using Manning's equation for a
  !! SFR package reach.
  !<
  subroutine sfr_calc_qman(this, n, depth, qman)
    ! -- 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
  end subroutine sfr_calc_qman
 
  !> @brief Calculate reach-aquifer exchange
  !!
  !! 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.
  !<
  subroutine sfr_calc_qgwf(this, n, depth, hgwf, qgwf, gwfhcof, gwfrhs)
    ! -- 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
  end subroutine sfr_calc_qgwf
 
  !> @brief Determine if a reach is connected to a gwf cell
  !!
  !! 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.
  !<
  function sfr_gwf_conn(this, n)
    ! -- 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
  end function sfr_gwf_conn
 
  !> @brief Calculate reach-aquifer conductance
  !!
  !! Method to calculate the reach-aquifer conductance for a SFR package reach.
  !<
  subroutine sfr_calc_cond(this, n, depth, cond, hsfr, h_temp)
    ! -- 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
  end subroutine sfr_calc_cond
 
  !> @brief Calculate diversion flow
  !!
  !! 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.
  !<
  subroutine sfr_calc_div(this, n, i, qd, qdiv)
    ! -- 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
  end subroutine sfr_calc_div
 
  !> @brief Calculate the depth at the midpoint
  !!
  !! Method to calculate the depth at the midpoint of a reach.
  !<
  subroutine sfr_calc_reach_depth(this, n, q1, d1)
    ! -- 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
  end subroutine sfr_calc_reach_depth
 
  !> @brief Calculate the depth at the midpoint of a irregular cross-section
  !!
  !! Method to calculate the depth at the midpoint of a reach with a
  !! irregular cross-section using Newton-Raphson.
  !<
  subroutine sfr_calc_xs_depth(this, n, qrch, d)
    ! -- 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
  end subroutine sfr_calc_xs_depth
 
  !> @brief Check unit conversion data
  !!
  !! Method to check unit conversion data for a SFR package. This method
  !! also calculates unitconv that is used in the Manning's equation.
  !<
  subroutine sfr_check_conversion(this)
    ! -- 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
  end subroutine sfr_check_conversion
 
  !> @brief Check storage weight
  !!
  !! 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.
  !<
  subroutine sfr_check_storage_weight(this)
    ! -- 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
  end subroutine sfr_check_storage_weight
 
  !> @brief Check reach data
  !!
  !! 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.
  !<
  subroutine sfr_check_reaches(this)
    ! -- 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
  end subroutine sfr_check_reaches
 
  !> @brief Check connection data
  !!
  !! 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.
  !<
  subroutine sfr_check_connections(this)
    ! -- 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
  end subroutine sfr_check_connections
 
  !> @brief Check diversions data
  !!
  !! 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.
  !<
  subroutine sfr_check_diversions(this)
    ! -- 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
  end subroutine sfr_check_diversions
 
  !> @brief Check initial stage data
  !!
  !! 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.
  !<
  subroutine sfr_check_initialstages(this)
    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
  end subroutine sfr_check_initialstages
 
  !> @brief Check upstream fraction data
  !!
  !! 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.
  !<
  subroutine sfr_check_ustrf(this)
    ! -- 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
  end subroutine sfr_check_ustrf
 
  !> @brief Apply a PERIOD setting value to reach n.
  !!
  !! Shared by sfr_rp and sfr_ad; DIVERSION is handled separately by
  !! both callers because its target index requires its own validation.
  !<
  subroutine sfr_set_period_value(this, i, n, setting)
    ! -- 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
  end subroutine sfr_set_period_value
 
  !> @brief Setup budget object for package
  !!
  !! 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.
  !<
  subroutine sfr_setup_budobj(this)
    ! -- 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
  end subroutine sfr_setup_budobj
 
  !> @brief Copy flow terms into budget object for package
  !!
  !! 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.
  !<
  subroutine sfr_fill_budobj(this)
    ! -- 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()
  end subroutine sfr_fill_budobj
 
  !> @brief Setup stage table object for package
  !!
  !! 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.
  !<
  subroutine sfr_setup_tableobj(this)
    ! -- 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
  end subroutine sfr_setup_tableobj
 
  ! -- reach geometry functions
 
  !> @brief Calculate wetted area
  !!
  !! Function to calculate the wetted area for a SFR package reach.
  !<
  function calc_area_wet(this, n, depth)
    ! -- 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
  end function calc_area_wet
 
  !> @brief Calculate wetted perimeter
  !!
  !! Function to calculate the wetted perimeter for a SFR package reach.
  !<
  function calc_perimeter_wet(this, n, depth)
    ! -- 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
  end function calc_perimeter_wet
 
  !> @brief Calculate maximum surface area
  !!
  !! Function to calculate the maximum surface area for a SFR package reach.
  !<
  function calc_surface_area(this, n)
    ! -- 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)
  end function calc_surface_area
 
  !> @brief Calculate wetted surface area
  !!
  !! Function to calculate the wetted surface area for a SFR package reach.
  !<
  function calc_surface_area_wet(this, n, depth)
    ! -- 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)
  end function calc_surface_area_wet
 
  !> @brief Calculate wetted top width
  !!
  !! Function to calculate the wetted top width for a SFR package reach.
  !<
  function calc_top_width_wet(this, n, depth)
    ! -- 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
  end function calc_top_width_wet
 
  !> @brief Activate density terms
  !!
  !! Method to activate addition of density terms for a SFR package reach.
  !<
  subroutine sfr_activate_density(this)
    ! -- 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))
  end subroutine sfr_activate_density
 
  !> @brief Activate viscosity terms
  !!
  !! Method to activate addition of viscosity terms for exchange
  !! with groundwater along a SFR package reach.
  !<
  subroutine sfr_activate_viscosity(this)
    ! -- 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))
  end subroutine sfr_activate_viscosity
 
  !> @brief Calculate density terms
  !!
  !! 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
  !<
  subroutine sfr_calculate_density_exchange(this, n, stage, head, cond, &
                                            tops, flow, gwfhcof, gwfrhs)
    ! -- 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
  end subroutine sfr_calculate_density_exchange
 
  !> @brief Bind SfrInputType PERIOD pointers to input context arrays.
  !!
  !! Called once from sfr_ar; PACKAGEDATA pointers are set separately
  !! in sfr_source_packagedata.
  !<
  subroutine sfr_input_init(this, mempath, naux)
    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
  end subroutine sfr_input_init
 
  !> @brief Nullify all SfrInputType pointers.
  !<
  subroutine sfr_input_destroy(this)
    class(sfrinputtype), intent(inout) :: this
    nullify (this%rhk)
    nullify (this%man)
    nullify (this%ustrf)
    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)
  end subroutine sfr_input_destroy
 
end module sfrmodule