MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
sfrmodule Module Reference

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

Data Types

type  sfrtype
 

Functions/Subroutines

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

Variables

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

Detailed Description

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

Function/Subroutine Documentation

◆ calc_area_wet()

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

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

Returns
wetted area
Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth

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

5554  ! -- return variable
5555  real(DP) :: calc_area_wet !< wetted area
5556  ! -- dummy
5557  class(SfrType) :: this !< SfrType object
5558  integer(I4B), intent(in) :: n !< reach number
5559  real(DP), intent(in) :: depth !< reach depth
5560  ! -- local
5561  integer(I4B) :: npts
5562  integer(I4B) :: i0
5563  integer(I4B) :: i1
5564  !
5565  ! -- Calculate wetted area
5566  npts = this%ncrosspts(n)
5567  i0 = this%iacross(n)
5568  i1 = this%iacross(n + 1) - 1
5569  if (npts > 1) then
5570  calc_area_wet = get_cross_section_area(npts, this%station(i0:i1), &
5571  this%xsheight(i0:i1), depth)
5572  else
5573  calc_area_wet = this%station(i0) * depth
5574  end if
Here is the call graph for this function:

◆ calc_perimeter_wet()

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

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

Returns
wetted perimeter
Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth

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

5582  ! -- return variable
5583  real(DP) :: calc_perimeter_wet !< wetted perimeter
5584  ! -- dummy
5585  class(SfrType) :: this !< SfrType object
5586  integer(I4B), intent(in) :: n !< reach number
5587  real(DP), intent(in) :: depth !< reach depth
5588  ! -- local
5589  integer(I4B) :: npts
5590  integer(I4B) :: i0
5591  integer(I4B) :: i1
5592  !
5593  ! -- Calculate wetted perimeter
5594  npts = this%ncrosspts(n)
5595  i0 = this%iacross(n)
5596  i1 = this%iacross(n + 1) - 1
5597  if (npts > 1) then
5598  calc_perimeter_wet = get_wetted_perimeter(npts, this%station(i0:i1), &
5599  this%xsheight(i0:i1), depth)
5600  else
5601  calc_perimeter_wet = this%station(i0) ! no depth dependence in original implementation
5602  end if
Here is the call graph for this function:

◆ calc_surface_area()

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

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

Returns
surface area
Parameters
thisSfrType object
[in]nreach number

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

5610  ! -- return variable
5611  real(DP) :: calc_surface_area !< surface area
5612  ! -- dummy
5613  class(SfrType) :: this !< SfrType object
5614  integer(I4B), intent(in) :: n !< reach number
5615  ! -- local
5616  integer(I4B) :: npts
5617  integer(I4B) :: i0
5618  integer(I4B) :: i1
5619  real(DP) :: top_width
5620  !
5621  ! -- Calculate surface area
5622  npts = this%ncrosspts(n)
5623  i0 = this%iacross(n)
5624  i1 = this%iacross(n + 1) - 1
5625  if (npts > 1) then
5626  top_width = get_saturated_topwidth(npts, this%station(i0:i1))
5627  else
5628  top_width = this%station(i0)
5629  end if
5630  calc_surface_area = top_width * this%length(n)
Here is the call graph for this function:

◆ calc_surface_area_wet()

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

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

Returns
wetted surface area
Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth

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

5638  ! -- return variable
5639  real(DP) :: calc_surface_area_wet !< wetted surface area
5640  ! -- dummy
5641  class(SfrType) :: this !< SfrType object
5642  integer(I4B), intent(in) :: n !< reach number
5643  real(DP), intent(in) :: depth !< reach depth
5644  ! -- local
5645  real(DP) :: top_width
5646  !
5647  ! -- Calculate wetted surface area
5648  top_width = this%calc_top_width_wet(n, depth)
5649  calc_surface_area_wet = top_width * this%length(n)

◆ calc_top_width_wet()

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

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

Returns
wetted top width
Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth

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

5657  ! -- return variable
5658  real(DP) :: calc_top_width_wet !< wetted top width
5659  ! -- dummy
5660  class(SfrType) :: this !< SfrType object
5661  integer(I4B), intent(in) :: n !< reach number
5662  real(DP), intent(in) :: depth !< reach depth
5663  ! -- local
5664  integer(I4B) :: npts
5665  integer(I4B) :: i0
5666  integer(I4B) :: i1
5667  real(DP) :: sat
5668  !
5669  ! -- Calculate wetted top width
5670  npts = this%ncrosspts(n)
5671  i0 = this%iacross(n)
5672  i1 = this%iacross(n + 1) - 1
5673  sat = scubicsaturation(dem5, dzero, depth, dem5)
5674  if (npts > 1) then
5675  calc_top_width_wet = sat * get_wetted_topwidth(npts, &
5676  this%station(i0:i1), &
5677  this%xsheight(i0:i1), &
5678  depth)
5679  else
5680  calc_top_width_wet = sat * this%station(i0)
5681  end if
Here is the call graph for this function:

◆ define_listlabel()

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

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

Parameters
[in,out]thisSfrType object

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

2923  ! -- dummy
2924  class(SfrType), intent(inout) :: this !< SfrType object
2925  !
2926  ! -- create the header list label
2927  this%listlabel = trim(this%filtyp)//' NO.'
2928  if (this%dis%ndim == 3) then
2929  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
2930  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
2931  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
2932  elseif (this%dis%ndim == 2) then
2933  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
2934  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
2935  else
2936  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
2937  end if
2938  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'STRESS RATE'
2939  if (this%inamedbound == 1) then
2940  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
2941  end if

◆ sfr_activate_density()

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

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

Parameters
[in,out]thisSfrType object

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

5689  ! -- modules
5691  ! -- dummy
5692  class(SfrType), intent(inout) :: this !< SfrType object
5693  ! -- local
5694  integer(I4B) :: i
5695  integer(I4B) :: j
5696  !
5697  ! -- Set idense and reallocate denseterms to be of size MAXBOUND
5698  this%idense = 1
5699  call mem_reallocate(this%denseterms, 3, this%MAXBOUND, 'DENSETERMS', &
5700  this%memoryPath)
5701  do i = 1, this%maxbound
5702  do j = 1, 3
5703  this%denseterms(j, i) = dzero
5704  end do
5705  end do
5706  write (this%iout, '(/1x,a)') 'DENSITY TERMS HAVE BEEN ACTIVATED FOR SFR &
5707  &PACKAGE: '//trim(adjustl(this%packName))

◆ sfr_activate_viscosity()

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

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

Parameters
[in,out]thisSfrType object

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

5716  ! -- modules
5718  ! -- dummy
5719  class(SfrType), intent(inout) :: this !< SfrType object
5720  ! -- local
5721  integer(I4B) :: i
5722  integer(I4B) :: j
5723  !
5724  ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND
5725  this%ivsc = 1
5726  call mem_reallocate(this%viscratios, 2, this%MAXBOUND, 'VISCRATIOS', &
5727  this%memoryPath)
5728  do i = 1, this%maxbound
5729  do j = 1, 2
5730  this%viscratios(j, i) = done
5731  end do
5732  end do
5733  write (this%iout, '(/1x,a)') 'VISCOSITY HAS BEEN ACTIVATED FOR SFR &
5734  &PACKAGE: '//trim(adjustl(this%packName))

◆ sfr_ad()

subroutine sfrmodule::sfr_ad ( class(sfrtype this)

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

Parameters
thisSfrType object

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

2078  ! -- modules
2080  ! -- dummy
2081  class(SfrType) :: this !< SfrType object
2082  ! -- local
2083  integer(I4B) :: n
2084  integer(I4B) :: iaux
2085 
2086  ! -- update previous values
2087  if (this%istorage == 1) then
2088  do n = 1, this%maxbound
2089  this%stageold(n) = this%stage(n)
2090  this%usinflowold(n) = this%usinflow(n)
2091  this%dsflowold(n) = this%dsflow(n)
2092  end do
2093  end if
2094  !
2095  ! -- Most advanced package AD routines have to restore state if
2096  ! the solution failed and the time step is being retried with a smaller
2097  ! step size. This is not needed here because there is no old stage
2098  ! or storage effects in the stream.
2099  !
2100  ! -- Advance the time series manager
2101  call this%TsManager%ad()
2102  !
2103  ! -- check upstream fractions if time series are being used to
2104  ! define this variable
2105  if (var_timeseries(this%tsManager, this%packName, 'USTRF')) then
2106  call this%sfr_check_ustrf()
2107  end if
2108  !
2109  ! -- update auxiliary variables by copying from the derived-type time
2110  ! series variable into the bndpackage auxvar variable so that this
2111  ! information is properly written to the GWF budget file
2112  if (this%naux > 0) then
2113  do n = 1, this%maxbound
2114  do iaux = 1, this%naux
2115  if (this%noupdateauxvar(iaux) /= 0) cycle
2116  this%auxvar(iaux, n) = this%rauxvar(iaux, n)
2117  end do
2118  end do
2119  end if
2120  !
2121  ! -- reset upstream flow to zero and set specified stage
2122  do n = 1, this%maxbound
2123  this%usflow(n) = dzero
2124  if (this%iboundpak(n) < 0) then
2125  this%stage(n) = this%sstage(n)
2126  end if
2127  end do
2128  !
2129  ! -- pakmvrobj ad
2130  if (this%imover == 1) then
2131  call this%pakmvrobj%ad()
2132  end if
2133  !
2134  ! -- For each observation, push simulated value and corresponding
2135  ! simulation time from "current" to "preceding" and reset
2136  ! "current" value.
2137  call this%obs%obs_ad()
logical function, public var_timeseries(tsManager, pkgName, varName, auxOrBnd)
Determine if a timeseries link with varName is defined.
Here is the call graph for this function:

◆ sfr_adjust_ro_ev()

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

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

Parameters
thisSfrType object
[in,out]qctotal reach volumetric flow
[in]quupstream reach volumetric flow
[in]qireach volumetric inflow
[in]qrreach volumetric rainfall
[in,out]qroreach volumetric runoff
[in,out]qereach volumetric evaporation
[in]qfrommvrreach volumetric flow from mover

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

3748  ! -- dummy
3749  class(SfrType) :: this !< SfrType object
3750  real(DP), intent(inout) :: qc !< total reach volumetric flow
3751  real(DP), intent(in) :: qu !< upstream reach volumetric flow
3752  real(DP), intent(in) :: qi !< reach volumetric inflow
3753  real(DP), intent(in) :: qr !< reach volumetric rainfall
3754  real(DP), intent(inout) :: qro !< reach volumetric runoff
3755  real(DP), intent(inout) :: qe !< reach volumetric evaporation
3756  real(DP), intent(in) :: qfrommvr !< reach volumetric flow from mover
3757  ! -- local
3758  real(DP) :: qt
3759  !
3760  ! -- adjust runoff or evaporation if sum of sources is negative
3761  if (qc < dzero) then
3762  !
3763  ! -- calculate sources without evaporation
3764  qt = qu + qi + qr + qro + qfrommvr
3765  !
3766  ! -- runoff exceeds sources of water for reach
3767  if (qt < dzero) then
3768  if (qro < dzero) then
3769  qro = -(qu + qi + qr + qfrommvr)
3770  qe = dzero
3771  end if
3772  !
3773  ! -- evaporation exceeds sources of water for reach
3774  else
3775  if (qe > dzero) then
3776  qe = qu + qi + qr + qro + qfrommvr
3777  end if
3778  end if
3779  qc = qu + qi + qr - qe + qro + qfrommvr
3780  end if

◆ sfr_allocate_arrays()

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

Allocate and initialize array for the SFR package.

Parameters
[in,out]thisSfrType object

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

407  ! -- modules
409  ! -- dummy
410  class(SfrType), intent(inout) :: this !< SfrType object
411  ! -- local
412  integer(I4B) :: i
413  integer(I4B) :: j
414  !
415  ! -- allocate character array for budget text
416  allocate (this%csfrbudget(this%bditems))
417  call mem_allocate(this%sfrname, lenboundname, this%maxbound, &
418  'SFRNAME', this%memoryPath)
419  !
420  ! -- variables originally in SfrDataType
421  call mem_allocate(this%iboundpak, this%maxbound, 'IBOUNDPAK', &
422  this%memoryPath)
423  call mem_allocate(this%igwfnode, this%maxbound, 'IGWFNODE', this%memoryPath)
424  call mem_allocate(this%igwftopnode, this%maxbound, 'IGWFTOPNODE', &
425  this%memoryPath)
426  call mem_allocate(this%length, this%maxbound, 'LENGTH', this%memoryPath)
427  call mem_allocate(this%width, this%maxbound, 'WIDTH', this%memoryPath)
428  call mem_allocate(this%strtop, this%maxbound, 'STRTOP', this%memoryPath)
429  call mem_allocate(this%bthick, this%maxbound, 'BTHICK', this%memoryPath)
430  call mem_allocate(this%hk, this%maxbound, 'HK', this%memoryPath)
431  call mem_allocate(this%slope, this%maxbound, 'SLOPE', this%memoryPath)
432  call mem_allocate(this%nconnreach, this%maxbound, 'NCONNREACH', &
433  this%memoryPath)
434  call mem_allocate(this%ustrf, this%maxbound, 'USTRF', this%memoryPath)
435  call mem_allocate(this%ftotnd, this%maxbound, 'FTOTND', this%memoryPath)
436  call mem_allocate(this%ndiv, this%maxbound, 'NDIV', this%memoryPath)
437  call mem_allocate(this%usflow, this%maxbound, 'USFLOW', this%memoryPath)
438  call mem_allocate(this%dsflow, this%maxbound, 'DSFLOW', this%memoryPath)
439  call mem_allocate(this%depth, this%maxbound, 'DEPTH', this%memoryPath)
440  call mem_allocate(this%stage, this%maxbound, 'STAGE', this%memoryPath)
441  call mem_allocate(this%gwflow, this%maxbound, 'GWFLOW', this%memoryPath)
442  call mem_allocate(this%simevap, this%maxbound, 'SIMEVAP', this%memoryPath)
443  call mem_allocate(this%simrunoff, this%maxbound, 'SIMRUNOFF', &
444  this%memoryPath)
445  call mem_allocate(this%stage0, this%maxbound, 'STAGE0', this%memoryPath)
446  call mem_allocate(this%usflow0, this%maxbound, 'USFLOW0', this%memoryPath)
447  !
448  ! -- stage, usflow, inflow, and dsflow for previous timestep
449  if (this%istorage == 1) then
450  call mem_allocate(this%stageold, this%maxbound, 'STAGEOLD', &
451  this%memoryPath)
452  call mem_allocate(this%usinflow, this%maxbound, 'USINFLOW', &
453  this%memoryPath)
454  call mem_allocate(this%usinflowold, this%maxbound, 'USINFLOWOLD', &
455  this%memoryPath)
456  call mem_allocate(this%dsflowold, this%maxbound, 'DSFLOWOLD', &
457  this%memoryPath)
458  call mem_allocate(this%storage, this%maxbound, 'STORAGE', &
459  this%memoryPath)
460  end if
461  !
462  ! -- reach order and connection data
463  call mem_allocate(this%isfrorder, this%maxbound, 'ISFRORDER', &
464  this%memoryPath)
465  call mem_allocate(this%ia, this%maxbound + 1, 'IA', this%memoryPath)
466  call mem_allocate(this%ja, 0, 'JA', this%memoryPath)
467  call mem_allocate(this%idir, 0, 'IDIR', this%memoryPath)
468  call mem_allocate(this%idiv, 0, 'IDIV', this%memoryPath)
469  call mem_allocate(this%qconn, 0, 'QCONN', this%memoryPath)
470  !
471  ! -- boundary data
472  call mem_allocate(this%rough, this%maxbound, 'ROUGH', this%memoryPath)
473  call mem_allocate(this%rain, this%maxbound, 'RAIN', this%memoryPath)
474  call mem_allocate(this%evap, this%maxbound, 'EVAP', this%memoryPath)
475  call mem_allocate(this%inflow, this%maxbound, 'INFLOW', this%memoryPath)
476  call mem_allocate(this%runoff, this%maxbound, 'RUNOFF', this%memoryPath)
477  call mem_allocate(this%sstage, this%maxbound, 'SSTAGE', this%memoryPath)
478  !
479  ! -- aux variables
480  call mem_allocate(this%rauxvar, this%naux, this%maxbound, &
481  'RAUXVAR', this%memoryPath)
482  !
483  ! -- diversion variables
484  call mem_allocate(this%iadiv, this%maxbound + 1, 'IADIV', this%memoryPath)
485  call mem_allocate(this%divreach, 0, 'DIVREACH', this%memoryPath)
486  call mem_allocate(this%divflow, 0, 'DIVFLOW', this%memoryPath)
487  call mem_allocate(this%divq, 0, 'DIVQ', this%memoryPath)
488  !
489  ! -- cross-section data
490  call mem_allocate(this%ncrosspts, this%maxbound, 'NCROSSPTS', &
491  this%memoryPath)
492  call mem_allocate(this%iacross, this%maxbound + 1, 'IACROSS', &
493  this%memoryPath)
494  call mem_allocate(this%station, this%ncrossptstot, 'STATION', &
495  this%memoryPath)
496  call mem_allocate(this%xsheight, this%ncrossptstot, 'XSHEIGHT', &
497  this%memoryPath)
498  call mem_allocate(this%xsrough, this%ncrossptstot, 'XSROUGH', &
499  this%memoryPath)
500  !
501  ! -- initialize variables
502  this%iacross(1) = 0
503  do i = 1, this%maxbound
504  this%iboundpak(i) = 1
505  this%igwfnode(i) = 0
506  this%igwftopnode(i) = 0
507  this%length(i) = dzero
508  this%width(i) = dzero
509  this%strtop(i) = dzero
510  this%bthick(i) = dzero
511  this%hk(i) = dzero
512  this%slope(i) = dzero
513  this%nconnreach(i) = 0
514  this%ustrf(i) = dzero
515  this%ftotnd(i) = dzero
516  this%ndiv(i) = 0
517  this%usflow(i) = dzero
518  this%dsflow(i) = dzero
519  this%depth(i) = dzero
520  this%stage(i) = dzero
521  this%gwflow(i) = dzero
522  this%simevap(i) = dzero
523  this%simrunoff(i) = dzero
524  this%stage0(i) = dzero
525  this%usflow0(i) = dzero
526  !
527  ! -- stage
528  if (this%istorage == 1) then
529  this%stageold(i) = dzero
530  this%usinflow(i) = dzero
531  this%usinflowold(i) = dzero
532  this%dsflowold(i) = dzero
533  this%storage(i) = dzero
534  end if
535  !
536  ! -- boundary data
537  this%rough(i) = dzero
538  this%rain(i) = dzero
539  this%evap(i) = dzero
540  this%inflow(i) = dzero
541  this%runoff(i) = dzero
542  this%sstage(i) = dzero
543  !
544  ! -- aux variables
545  do j = 1, this%naux
546  this%rauxvar(j, i) = dzero
547  end do
548  !
549  ! -- cross-section data
550  this%ncrosspts(i) = 0
551  this%iacross(i + 1) = 0
552  end do
553  !
554  ! -- initialize additional cross-section data
555  do i = 1, this%ncrossptstot
556  this%station(i) = dzero
557  this%xsheight(i) = dzero
558  this%xsrough(i) = dzero
559  end do
560  !
561  !-- fill csfrbudget
562  this%csfrbudget(1) = ' RAINFALL'
563  this%csfrbudget(2) = ' EVAPORATION'
564  this%csfrbudget(3) = ' RUNOFF'
565  this%csfrbudget(4) = ' EXT-INFLOW'
566  this%csfrbudget(5) = ' GWF'
567  this%csfrbudget(6) = ' EXT-OUTFLOW'
568  this%csfrbudget(7) = ' FROM-MVR'
569  this%csfrbudget(8) = ' TO-MVR'
570  !
571  ! -- allocate and initialize budget output data
572  call mem_allocate(this%qoutflow, this%maxbound, 'QOUTFLOW', this%memoryPath)
573  call mem_allocate(this%qextoutflow, this%maxbound, 'QEXTOUTFLOW', &
574  this%memoryPath)
575  do i = 1, this%maxbound
576  this%qoutflow(i) = dzero
577  this%qextoutflow(i) = dzero
578  end do
579  !
580  ! -- allocate and initialize dbuff
581  if (this%istageout > 0) then
582  call mem_allocate(this%dbuff, this%maxbound, 'DBUFF', this%memoryPath)
583  do i = 1, this%maxbound
584  this%dbuff(i) = dzero
585  end do
586  else
587  call mem_allocate(this%dbuff, 0, 'DBUFF', this%memoryPath)
588  end if
589  !
590  ! -- allocate character array for budget text
591  allocate (this%cauxcbc(this%cbcauxitems))
592  !
593  ! -- allocate and initialize qauxcbc
594  call mem_allocate(this%qauxcbc, this%cbcauxitems, 'QAUXCBC', &
595  this%memoryPath)
596  do i = 1, this%cbcauxitems
597  this%qauxcbc(i) = dzero
598  end do
599  !
600  ! -- fill cauxcbc
601  this%cauxcbc(1) = 'FLOW-AREA '
602  !
603  ! -- allocate denseterms to size 0
604  call mem_allocate(this%denseterms, 3, 0, 'DENSETERMS', this%memoryPath)
605  !
606  ! -- allocate viscratios to size 0
607  call mem_allocate(this%viscratios, 2, 0, 'VISCRATIOS', this%memoryPath)

◆ sfr_allocate_scalars()

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

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

Parameters
[in,out]thisSfrType object

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

338  ! -- modules
341  ! -- dummy
342  class(SfrType), intent(inout) :: this !< SfrType object
343  !
344  ! -- call standard BndType allocate scalars
345  call this%BndType%allocate_scalars()
346  !
347  ! -- allocate the object and assign values to object variables
348  call mem_allocate(this%istorage, 'ISTORAGE', this%memoryPath)
349  call mem_allocate(this%iprhed, 'IPRHED', this%memoryPath)
350  call mem_allocate(this%istageout, 'ISTAGEOUT', this%memoryPath)
351  call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
352  call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
353  call mem_allocate(this%ipakcsv, 'IPAKCSV', this%memoryPath)
354  call mem_allocate(this%idiversions, 'IDIVERSIONS', this%memoryPath)
355  call mem_allocate(this%maxsfrpicard, 'MAXSFRPICARD', this%memoryPath)
356  call mem_allocate(this%maxsfrit, 'MAXSFRIT', this%memoryPath)
357  call mem_allocate(this%bditems, 'BDITEMS', this%memoryPath)
358  call mem_allocate(this%cbcauxitems, 'CBCAUXITEMS', this%memoryPath)
359  call mem_allocate(this%unitconv, 'UNITCONV', this%memoryPath)
360  call mem_allocate(this%lengthconv, 'LENGTHCONV', this%memoryPath)
361  call mem_allocate(this%timeconv, 'TIMECONV', this%memoryPath)
362  call mem_allocate(this%dmaxchg, 'DMAXCHG', this%memoryPath)
363  call mem_allocate(this%deps, 'DEPS', this%memoryPath)
364  call mem_allocate(this%storage_weight, 'STORAGE_WEIGHT', this%memoryPath)
365  call mem_allocate(this%nconn, 'NCONN', this%memoryPath)
366  call mem_allocate(this%icheck, 'ICHECK', this%memoryPath)
367  call mem_allocate(this%iconvchk, 'ICONVCHK', this%memoryPath)
368  call mem_allocate(this%idense, 'IDENSE', this%memoryPath)
369  call mem_allocate(this%ianynone, 'IANYNONE', this%memoryPath)
370  call mem_allocate(this%ncrossptstot, 'NCROSSPTSTOT', this%memoryPath)
371  !
372  ! -- set pointer to gwf iss
373  call mem_setptr(this%gwfiss, 'ISS', create_mem_path(this%name_model))
374  !
375  ! -- Set values
376  this%istorage = 0
377  this%iprhed = 0
378  this%istageout = 0
379  this%ibudgetout = 0
380  this%ibudcsv = 0
381  this%ipakcsv = 0
382  this%idiversions = 0
383  this%maxsfrpicard = 100
384  this%maxsfrit = maxadpit
385  this%bditems = 8
386  this%cbcauxitems = 1
387  this%unitconv = done
388  this%lengthconv = dnodata
389  this%timeconv = dnodata
390  this%dmaxchg = dem5
391  this%deps = dp999 * this%dmaxchg
392  this%storage_weight = dnodata
393  this%nconn = 0
394  this%icheck = 1
395  this%iconvchk = 1
396  this%idense = 0
397  this%ivsc = 0
398  this%ianynone = 0
399  this%ncrossptstot = 0
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
Here is the call graph for this function:

◆ sfr_ar()

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

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

Parameters
[in,out]thisSfrType object

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

867  ! -- dummy
868  class(SfrType), intent(inout) :: this !< SfrType object
869  ! -- local
870  integer(I4B) :: n
871  integer(I4B) :: ierr
872  !
873  ! -- allocate and read observations
874  call this%obs%obs_ar()
875  !
876  ! -- call standard BndType allocate scalars
877  call this%BndType%allocate_arrays()
878  !
879  ! -- set boundname for each connection
880  if (this%inamedbound /= 0) then
881  do n = 1, this%maxbound
882  this%boundname(n) = this%sfrname(n)
883  end do
884  end if
885  !
886  ! -- copy boundname into boundname_cst
887  call this%copy_boundname()
888  !
889  ! -- copy igwfnode into nodelist
890  do n = 1, this%maxbound
891  this%nodelist(n) = this%igwfnode(n)
892  end do
893  !
894  ! -- check the sfr unit conversion data
895  call this%sfr_check_conversion()
896  !
897  ! -- check the storage_weight
898  call this%sfr_check_storage_weight()
899  !
900  ! -- check the sfr reach data
901  call this%sfr_check_reaches()
902 
903  ! -- check the connection data
904  call this%sfr_check_connections()
905 
906  ! -- check the diversion data
907  if (this%idiversions /= 0) then
908  call this%sfr_check_diversions()
909  end if
910 
911  ! -- check the diversion data
912  if (this%istorage == 1) then
913  call this%sfr_check_initialstages()
914  end if
915  !
916  ! -- terminate if errors were detected in any of the static sfr data
917  ierr = count_errors()
918  if (ierr > 0) then
919  call this%parser%StoreErrorUnit()
920  end if
921  !
922  ! -- setup pakmvrobj
923  if (this%imover /= 0) then
924  allocate (this%pakmvrobj)
925  call this%pakmvrobj%ar(this%maxbound, this%maxbound, this%memoryPath)
926  end if
Here is the call graph for this function:

◆ sfr_bd_obs()

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

Method to save simulated values for the SFR package.

Parameters
thisSfrType object

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

3063  ! -- dummy
3064  class(SfrType) :: this !< SfrType object
3065  ! -- local
3066  integer(I4B) :: i
3067  integer(I4B) :: j
3068  integer(I4B) :: n
3069  real(DP) :: v
3070  character(len=100) :: msg
3071  type(ObserveType), pointer :: obsrv => null()
3072  !
3073  ! Write simulated values for all sfr observations
3074  if (this%obs%npakobs > 0) then
3075  call this%obs%obs_bd_clear()
3076  do i = 1, this%obs%npakobs
3077  obsrv => this%obs%pakobs(i)%obsrv
3078  do j = 1, obsrv%indxbnds_count
3079  n = obsrv%indxbnds(j)
3080  v = dzero
3081  select case (obsrv%ObsTypeId)
3082  case ('STAGE')
3083  v = this%stage(n)
3084  case ('TO-MVR')
3085  v = dnodata
3086  if (this%imover == 1) then
3087  v = this%pakmvrobj%get_qtomvr(n)
3088  if (v > dzero) then
3089  v = -v
3090  end if
3091  end if
3092  case ('FROM-MVR')
3093  v = dnodata
3094  if (this%imover == 1) then
3095  v = this%pakmvrobj%get_qfrommvr(n)
3096  end if
3097  case ('EXT-INFLOW')
3098  v = this%inflow(n)
3099  case ('INFLOW')
3100  v = this%usflow(n)
3101  case ('OUTFLOW')
3102  v = this%qoutflow(n)
3103  case ('EXT-OUTFLOW')
3104  v = this%qextoutflow(n)
3105  case ('RAINFALL')
3106  if (this%iboundpak(n) /= 0) then
3107  v = this%rain(n)
3108  else
3109  v = dzero
3110  end if
3111  case ('RUNOFF')
3112  v = this%simrunoff(n)
3113  case ('EVAPORATION')
3114  v = this%simevap(n)
3115  case ('SFR')
3116  v = this%gwflow(n)
3117  case ('UPSTREAM-FLOW')
3118  v = this%usflow(n)
3119  if (this%imover == 1) then
3120  v = v + this%pakmvrobj%get_qfrommvr(n)
3121  end if
3122  case ('DOWNSTREAM-FLOW')
3123  v = this%dsflow(n)
3124  if (v > dzero) then
3125  v = -v
3126  end if
3127  case ('DEPTH')
3128  v = this%depth(n)
3129  case ('WET-PERIMETER')
3130  v = this%calc_perimeter_wet(n, this%depth(n))
3131  case ('WET-AREA')
3132  v = this%calc_area_wet(n, this%depth(n))
3133  case ('WET-WIDTH')
3134  v = this%calc_top_width_wet(n, this%depth(n))
3135  case default
3136  msg = 'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
3137  call store_error(msg)
3138  end select
3139  call this%obs%SaveOneSimval(obsrv, v)
3140  end do
3141  end do
3142  !
3143  ! -- write summary of package error messages
3144  if (count_errors() > 0) then
3145  call this%parser%StoreErrorUnit()
3146  end if
3147  end if
Here is the call graph for this function:

◆ sfr_calc_cond()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth
[in,out]condreach-aquifer conductance
[in]hsfrstream stage
[in]h_temphead in gw cell

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

4026  ! -- dummy
4027  class(SfrType) :: this !< SfrType object
4028  integer(I4B), intent(in) :: n !< reach number
4029  real(DP), intent(in) :: depth !< reach depth
4030  real(DP), intent(inout) :: cond !< reach-aquifer conductance
4031  real(DP), intent(in), optional :: hsfr !< stream stage
4032  real(DP), intent(in), optional :: h_temp !< head in gw cell
4033  ! -- local
4034  integer(I4B) :: node
4035  real(DP) :: wp
4036  real(DP) :: vscratio
4037  !
4038  ! -- initialize conductance
4039  cond = dzero
4040  !
4041  ! -- initial viscosity ratio to 1
4042  vscratio = done
4043  !
4044  ! -- calculate conductance if GWF cell is active
4045  ! rch-gwf flow will not occur if reach connected to an constant head cell
4046  node = this%igwfnode(n)
4047  if (node > 0) then
4048  if (this%ibound(node) > 0) then
4049  !
4050  ! -- direction of gradient across streambed determines which vsc ratio
4051  if (this%ivsc == 1) then
4052  if (hsfr > h_temp) then
4053  ! strm stg > gw head
4054  vscratio = this%viscratios(1, n)
4055  else
4056  vscratio = this%viscratios(2, n)
4057  end if
4058  end if
4059  wp = this%calc_perimeter_wet(n, depth)
4060  cond = this%hk(n) * vscratio * this%length(n) * wp / this%bthick(n)
4061  end if
4062  end if

◆ sfr_calc_div()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]idiversion number in reach
[in,out]qdremaining downstream flow for reach
[in,out]qdivdiversion flow for diversion i

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

4073  ! -- dummy
4074  class(SfrType) :: this !< SfrType object
4075  integer(I4B), intent(in) :: n !< reach number
4076  integer(I4B), intent(in) :: i !< diversion number in reach
4077  real(DP), intent(inout) :: qd !< remaining downstream flow for reach
4078  real(DP), intent(inout) :: qdiv !< diversion flow for diversion i
4079  ! -- local
4080  character(len=10) :: cp
4081  integer(I4B) :: jpos
4082  integer(I4B) :: n2
4083  real(DP) :: v
4084  !
4085  ! -- set local variables
4086  jpos = this%iadiv(n) + i - 1
4087  n2 = this%divreach(jpos)
4088  cp = this%divcprior(jpos)
4089  v = this%divflow(jpos)
4090  !
4091  ! -- calculate diversion
4092  select case (cp)
4093  ! -- flood diversion
4094  case ('EXCESS')
4095  if (qd < v) then
4096  v = dzero
4097  else
4098  v = qd - v
4099  end if
4100  ! -- diversion percentage
4101  case ('FRACTION')
4102  v = qd * v
4103  ! -- STR priority algorithm
4104  case ('THRESHOLD')
4105  if (qd < v) then
4106  v = dzero
4107  end if
4108  ! -- specified diversion
4109  case ('UPTO')
4110  if (v > qd) then
4111  v = qd
4112  end if
4113  case default
4114  v = dzero
4115  end select
4116  !
4117  ! -- update upstream from for downstream reaches
4118  qd = qd - v
4119  qdiv = v

◆ sfr_calc_qd()

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

Method to calculate downstream flow for a SFR package reach.

Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth
[in]hgwfgroundwater head in connected GWF cell
[in,out]qgwfgroundwater leakage for reach
[in,out]qdresidual

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

3788  ! -- dummy
3789  class(SfrType) :: this !< SfrType object
3790  integer(I4B), intent(in) :: n !< reach number
3791  real(DP), intent(in) :: depth !< reach depth
3792  real(DP), intent(in) :: hgwf !< groundwater head in connected GWF cell
3793  real(DP), intent(inout) :: qgwf !< groundwater leakage for reach
3794  real(DP), intent(inout) :: qd !< residual
3795  ! -- local
3796  real(DP) :: qsrc
3797  !
3798  ! -- initialize residual
3799  qd = dzero
3800  !
3801  ! -- calculate total water sources excluding groundwater leakage
3802  call this%sfr_calc_qsource(n, depth, qsrc)
3803  !
3804  ! -- estimate groundwater leakage
3805  call this%sfr_calc_qgwf(n, depth, hgwf, qgwf)
3806  if (-qgwf > qsrc) qgwf = -qsrc
3807  !
3808  ! -- calculate down stream flow
3809  qd = qsrc + qgwf
3810  !
3811  ! -- limit downstream flow to a positive value
3812  if (qd < dem30) qd = dzero

◆ sfr_calc_qgwf()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth
[in]hgwfhead in GWF cell connected to reach
[in,out]qgwfreach-aquifer exchange
[in,out]gwfhcofdiagonal coefficient term for reach
[in,out]gwfrhsright-hand side term for reach

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

3939  ! -- dummy
3940  class(SfrType) :: this !< SfrType object
3941  integer(I4B), intent(in) :: n !< reach number
3942  real(DP), intent(in) :: depth !< reach depth
3943  real(DP), intent(in) :: hgwf !< head in GWF cell connected to reach
3944  real(DP), intent(inout) :: qgwf !< reach-aquifer exchange
3945  real(DP), intent(inout), optional :: gwfhcof !< diagonal coefficient term for reach
3946  real(DP), intent(inout), optional :: gwfrhs !< right-hand side term for reach
3947  ! -- local
3948  integer(I4B) :: node
3949  real(DP) :: tp
3950  real(DP) :: bt
3951  real(DP) :: hsfr
3952  real(DP) :: h_temp
3953  real(DP) :: cond
3954  real(DP) :: sat
3955  real(DP) :: derv
3956  real(DP) :: gwfhcof0
3957  real(DP) :: gwfrhs0
3958  !
3959  ! -- initialize qgwf
3960  qgwf = dzero
3961  !
3962  ! -- skip sfr-aquifer exchange in external cells
3963  node = this%igwfnode(n)
3964  if (node < 1) return
3965  !
3966  ! -- skip sfr-aquifer exchange in inactive cells
3967  if (this%ibound(node) == 0) return
3968  !
3969  ! -- calculate saturation
3970  call schsmooth(depth, sat, derv)
3971  !
3972  ! -- terms for calculating direction of gradient across streambed
3973  tp = this%strtop(n)
3974  bt = tp - this%bthick(n)
3975  hsfr = tp + depth
3976  h_temp = hgwf
3977  if (h_temp < bt) then
3978  h_temp = bt
3979  end if
3980  !
3981  ! -- calculate conductance
3982  call this%sfr_calc_cond(n, depth, cond, hsfr, h_temp)
3983  !
3984  ! -- calculate groundwater leakage
3985  qgwf = sat * cond * (h_temp - hsfr)
3986  gwfrhs0 = -sat * cond * hsfr
3987  gwfhcof0 = -sat * cond
3988  !
3989  ! Add density contributions, if active
3990  if (this%idense /= 0) then
3991  call this%sfr_calculate_density_exchange(n, hsfr, hgwf, cond, tp, &
3992  qgwf, gwfhcof0, gwfrhs0)
3993  end if
3994  !
3995  ! -- Set gwfhcof and gwfrhs if present
3996  if (present(gwfhcof)) gwfhcof = gwfhcof0
3997  if (present(gwfrhs)) gwfrhs = gwfrhs0
Here is the call graph for this function:

◆ sfr_calc_qman()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth
[in,out]qmanstreamflow

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

3869  ! -- dummy
3870  class(SfrType) :: this !< SfrType object
3871  integer(I4B), intent(in) :: n !< reach number
3872  real(DP), intent(in) :: depth !< reach depth
3873  real(DP), intent(inout) :: qman !< streamflow
3874  ! -- local
3875  integer(I4B) :: npts
3876  integer(I4B) :: i0
3877  integer(I4B) :: i1
3878  real(DP) :: sat
3879  real(DP) :: derv
3880  real(DP) :: s
3881  real(DP) :: r
3882  real(DP) :: aw
3883  real(DP) :: wp
3884  real(DP) :: rh
3885  !
3886  ! -- initialize variables
3887  qman = dzero
3888  !
3889  ! -- calculate Manning's discharge for non-zero depths
3890  if (depth > dzero) then
3891  npts = this%ncrosspts(n)
3892  !
3893  ! -- set constant terms for Manning's equation
3894  call schsmooth(depth, sat, derv)
3895  s = this%slope(n)
3896  !
3897  ! -- calculate the mannings coefficient that is a
3898  ! function of depth
3899  if (npts > 1) then
3900  !
3901  ! -- get the location of the cross-section data for the reach
3902  i0 = this%iacross(n)
3903  i1 = this%iacross(n + 1) - 1
3904  !
3905  ! -- get the Manning's sum of the Manning's discharge
3906  ! for each section
3907  qman = get_mannings_section(npts, &
3908  this%station(i0:i1), &
3909  this%xsheight(i0:i1), &
3910  this%xsrough(i0:i1), &
3911  this%rough(n), &
3912  this%unitconv, &
3913  s, &
3914  depth)
3915  else
3916  r = this%rough(n)
3917  aw = this%calc_area_wet(n, depth)
3918  wp = this%calc_perimeter_wet(n, depth)
3919  if (wp > dzero) then
3920  rh = aw / wp
3921  else
3922  rh = dzero
3923  end if
3924  qman = this%unitconv * aw * (rh**dtwothirds) * sqrt(s) / r
3925  end if
3926  !
3927  ! -- calculate stream flow
3928  qman = sat * qman
3929  end if
Here is the call graph for this function:

◆ sfr_calc_qsource()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]depthreach depth
[in,out]qsrcsum of sources for reach

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

3821  ! -- dummy
3822  class(SfrType) :: this !< SfrType object
3823  integer(I4B), intent(in) :: n !< reach number
3824  real(DP), intent(in) :: depth !< reach depth
3825  real(DP), intent(inout) :: qsrc !< sum of sources for reach
3826  ! -- local
3827  real(DP) :: qu
3828  real(DP) :: qi
3829  real(DP) :: qr
3830  real(DP) :: qe
3831  real(DP) :: qro
3832  real(DP) :: qfrommvr
3833  real(DP) :: a
3834  real(DP) :: ae
3835  !
3836  ! -- initialize residual
3837  qsrc = dzero
3838  !
3839  ! -- calculate flow terms
3840  qu = this%usflow(n)
3841  qi = this%inflow(n)
3842  qro = this%runoff(n)
3843  !
3844  ! -- calculate rainfall and evap
3845  a = this%calc_surface_area(n)
3846  ae = this%calc_surface_area_wet(n, depth)
3847  qr = this%rain(n) * a
3848  qe = this%evap(n) * ae
3849  !
3850  ! -- calculate mover term
3851  qfrommvr = dzero
3852  if (this%imover == 1) then
3853  qfrommvr = this%pakmvrobj%get_qfrommvr(n)
3854  end if
3855  !
3856  ! -- calculate down stream flow
3857  qsrc = qu + qi + qr - qe + qro + qfrommvr
3858  !
3859  ! -- adjust runoff or evaporation if sum of sources is negative
3860  call this%sfr_adjust_ro_ev(qsrc, qu, qi, qr, qro, qe, qfrommvr)

◆ sfr_calc_reach_depth()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]q1streamflow
[in,out]d1stream depth at midpoint of reach

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

4127  ! -- dummy
4128  class(SfrType) :: this !< SfrType object
4129  integer(I4B), intent(in) :: n !< reach number
4130  real(DP), intent(in) :: q1 !< streamflow
4131  real(DP), intent(inout) :: d1 !< stream depth at midpoint of reach
4132  ! -- local
4133  real(DP) :: w
4134  real(DP) :: s
4135  real(DP) :: r
4136  real(DP) :: qconst
4137  !
4138  ! -- initialize slope and roughness
4139  s = this%slope(n)
4140  r = this%rough(n)
4141  !
4142  ! -- calculate stream depth at the midpoint
4143  if (q1 > dzero) then
4144  if (this%ncrosspts(n) > 1) then
4145  call this%sfr_calc_xs_depth(n, q1, d1)
4146  else
4147  w = this%station(this%iacross(n))
4148  qconst = this%unitconv * w * sqrt(s) / r
4149  d1 = (q1 / qconst)**dp6
4150  end if
4151  else
4152  d1 = dzero
4153  end if

◆ sfr_calc_xs_depth()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]qrchstreamflow
[in,out]dstream depth at midpoint of reach

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

4162  ! -- dummy
4163  class(SfrType) :: this !< SfrType object
4164  integer(I4B), intent(in) :: n !< reach number
4165  real(DP), intent(in) :: qrch !< streamflow
4166  real(DP), intent(inout) :: d !< stream depth at midpoint of reach
4167  ! -- local
4168  integer(I4B) :: iter
4169  real(DP) :: perturbation
4170  real(DP) :: q0
4171  real(DP) :: q1
4172  real(DP) :: dq
4173  real(DP) :: derv
4174  real(DP) :: dd
4175  real(DP) :: residual
4176  !
4177  ! -- initialize variables
4178  perturbation = this%deps * dtwo
4179  d = dzero
4180  q0 = dzero
4181  residual = q0 - qrch
4182  !
4183  ! -- Newton-Raphson iteration
4184  nriter: do iter = 1, this%maxsfrit
4185  call this%sfr_calc_qman(n, d + perturbation, q1)
4186  dq = (q1 - q0)
4187  if (dq /= dzero) then
4188  derv = perturbation / (q1 - q0)
4189  else
4190  derv = dzero
4191  end if
4192  dd = derv * residual
4193  d = d - dd
4194  call this%sfr_calc_qman(n, d, q0)
4195  residual = q0 - qrch
4196  !
4197  ! -- check for convergence
4198  if (abs(dd) < this%dmaxchg) then
4199  exit nriter
4200  end if
4201  end do nriter

◆ sfr_calculate_density_exchange()

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

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

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

Parameters
[in,out]thisSfrType object
[in]nreach number
[in]stagereach stage
[in]headhead in connected GWF cell
[in]condreach conductance
[in]botsbottom elevation of reach
[in,out]flowcalculated flow, updated here with density terms
[in,out]gwfhcofGWF diagonal coefficient, updated here with density terms
[in,out]gwfrhsGWF right-hand-side value, updated here with density terms

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

5750  ! -- dummy
5751  class(SfrType), intent(inout) :: this !< SfrType object
5752  integer(I4B), intent(in) :: n !< reach number
5753  real(DP), intent(in) :: stage !< reach stage
5754  real(DP), intent(in) :: head !< head in connected GWF cell
5755  real(DP), intent(in) :: cond !< reach conductance
5756  real(DP), intent(in) :: bots !< bottom elevation of reach
5757  real(DP), intent(inout) :: flow !< calculated flow, updated here with density terms
5758  real(DP), intent(inout) :: gwfhcof !< GWF diagonal coefficient, updated here with density terms
5759  real(DP), intent(inout) :: gwfrhs !< GWF right-hand-side value, updated here with density terms
5760  ! -- local
5761  real(DP) :: ss
5762  real(DP) :: hh
5763  real(DP) :: havg
5764  real(DP) :: rdensesfr
5765  real(DP) :: rdensegwf
5766  real(DP) :: rdenseavg
5767  real(DP) :: elevsfr
5768  real(DP) :: elevgwf
5769  real(DP) :: elevavg
5770  real(DP) :: d1
5771  real(DP) :: d2
5772  logical(LGP) :: stage_below_bot
5773  logical(LGP) :: head_below_bot
5774  !
5775  ! -- Set sfr density to sfr density or gwf density
5776  if (stage >= bots) then
5777  ss = stage
5778  stage_below_bot = .false.
5779  rdensesfr = this%denseterms(1, n) ! sfr rel density
5780  else
5781  ss = bots
5782  stage_below_bot = .true.
5783  rdensesfr = this%denseterms(2, n) ! gwf rel density
5784  end if
5785  !
5786  ! -- set hh to head or bots
5787  if (head >= bots) then
5788  hh = head
5789  head_below_bot = .false.
5790  rdensegwf = this%denseterms(2, n) ! gwf rel density
5791  else
5792  hh = bots
5793  head_below_bot = .true.
5794  rdensegwf = this%denseterms(1, n) ! sfr rel density
5795  end if
5796  !
5797  ! -- todo: hack because denseterms not updated in a cf calculation
5798  if (rdensegwf == dzero) return
5799  !
5800  ! -- Update flow
5801  if (stage_below_bot .and. head_below_bot) then
5802  !
5803  ! -- flow is zero, so no terms are updated
5804  !
5805  else
5806  !
5807  ! -- calculate average relative density
5808  rdenseavg = dhalf * (rdensesfr + rdensegwf)
5809  !
5810  ! -- Add contribution of first density term:
5811  ! cond * (denseavg/denseref - 1) * (hgwf - hsfr)
5812  d1 = cond * (rdenseavg - done)
5813  gwfhcof = gwfhcof - d1
5814  gwfrhs = gwfrhs - d1 * ss
5815  d1 = d1 * (hh - ss)
5816  flow = flow + d1
5817  !
5818  ! -- Add second density term if stage and head not below bottom
5819  if (.not. stage_below_bot .and. .not. head_below_bot) then
5820  !
5821  ! -- Add contribution of second density term:
5822  ! cond * (havg - elevavg) * (densegwf - densesfr) / denseref
5823  elevgwf = this%denseterms(3, n)
5824  elevsfr = bots
5825  elevavg = dhalf * (elevsfr + elevgwf)
5826  havg = dhalf * (hh + ss)
5827  d2 = cond * (havg - elevavg) * (rdensegwf - rdensesfr)
5828  gwfrhs = gwfrhs + d2
5829  flow = flow + d2
5830  end if
5831  end if

◆ sfr_cc()

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

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

Parameters
[in,out]thisSfrType object
[in]innertottotal number of inner iterations
[in]kiterPicard iteration number
[in]iendflag indicating if this is the last Picard iteration
[in]icnvgmodflag inficating if the model has met specific convergence criteria
[in,out]cpakstring for user node
[in,out]ipaklocation of the maximum dependent variable change
[in,out]dpakmaximum dependent variable change

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

2324  ! -- modules
2325  use tdismodule, only: totim, kstp, kper, delt
2326  ! -- dummy
2327  class(SfrType), intent(inout) :: this !< SfrType object
2328  integer(I4B), intent(in) :: innertot !< total number of inner iterations
2329  integer(I4B), intent(in) :: kiter !< Picard iteration number
2330  integer(I4B), intent(in) :: iend !< flag indicating if this is the last Picard iteration
2331  integer(I4B), intent(in) :: icnvgmod !< flag inficating if the model has met specific convergence criteria
2332  character(len=LENPAKLOC), intent(inout) :: cpak !< string for user node
2333  integer(I4B), intent(inout) :: ipak !< location of the maximum dependent variable change
2334  real(DP), intent(inout) :: dpak !< maximum dependent variable change
2335  ! -- local
2336  character(len=LENPAKLOC) :: cloc
2337  character(len=LINELENGTH) :: tag
2338  integer(I4B) :: icheck
2339  integer(I4B) :: ipakfail
2340  integer(I4B) :: locdhmax
2341  integer(I4B) :: locrmax
2342  integer(I4B) :: locdqfrommvrmax
2343  integer(I4B) :: ntabrows
2344  integer(I4B) :: ntabcols
2345  integer(I4B) :: n
2346  real(DP) :: q
2347  real(DP) :: q0
2348  real(DP) :: qtolfact
2349  real(DP) :: dh
2350  real(DP) :: r
2351  real(DP) :: dhmax
2352  real(DP) :: rmax
2353  real(DP) :: dqfrommvr
2354  real(DP) :: dqfrommvrmax
2355  !
2356  ! -- initialize local variables
2357  icheck = this%iconvchk
2358  ipakfail = 0
2359  locdhmax = 0
2360  locrmax = 0
2361  r = dzero
2362  dhmax = dzero
2363  rmax = dzero
2364  locdqfrommvrmax = 0
2365  dqfrommvrmax = dzero
2366  !
2367  ! -- if not saving package convergence data on check convergence if
2368  ! the model is considered converged
2369  if (this%ipakcsv == 0) then
2370  if (icnvgmod == 0) then
2371  icheck = 0
2372  end if
2373  !
2374  ! -- saving package convergence data
2375  else
2376  !
2377  ! -- header for package csv
2378  if (.not. associated(this%pakcsvtab)) then
2379  !
2380  ! -- determine the number of columns and rows
2381  ntabrows = 1
2382  ntabcols = 9
2383  if (this%imover == 1) then
2384  ntabcols = ntabcols + 2
2385  end if
2386  !
2387  ! -- setup table
2388  call table_cr(this%pakcsvtab, this%packName, '')
2389  call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
2390  lineseparator=.false., separator=',', &
2391  finalize=.false.)
2392  !
2393  ! -- add columns to package csv
2394  tag = 'total_inner_iterations'
2395  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
2396  tag = 'totim'
2397  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
2398  tag = 'kper'
2399  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
2400  tag = 'kstp'
2401  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
2402  tag = 'nouter'
2403  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
2404  tag = 'dvmax'
2405  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
2406  tag = 'dvmax_loc'
2407  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
2408  tag = 'dinflowmax'
2409  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
2410  tag = 'dinflowmax_loc'
2411  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
2412  if (this%imover == 1) then
2413  tag = 'dqfrommvrmax'
2414  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
2415  tag = 'dqfrommvrmax_loc'
2416  call this%pakcsvtab%initialize_column(tag, 16, alignment=tableft)
2417  end if
2418  end if
2419  end if
2420  !
2421  ! -- perform package convergence check
2422  if (icheck /= 0) then
2423  final_check: do n = 1, this%maxbound
2424  if (this%iboundpak(n) == 0) cycle
2425  !
2426  ! -- set the Q to length factor
2427  qtolfact = delt / this%calc_surface_area(n)
2428  !
2429  ! -- calculate stage change
2430  dh = this%stage0(n) - this%stage(n)
2431  !
2432  ! -- evaluate flow difference if the time step is transient
2433  if (this%gwfiss == 0) then
2434  r = this%usflow0(n) - this%usflow(n)
2435  !
2436  ! -- normalize flow difference and convert to a depth
2437  r = r * qtolfact
2438  end if
2439  !
2440  ! -- q from mvr
2441  dqfrommvr = dzero
2442  if (this%imover == 1) then
2443  q = this%pakmvrobj%get_qfrommvr(n)
2444  q0 = this%pakmvrobj%get_qfrommvr0(n)
2445  dqfrommvr = qtolfact * (q0 - q)
2446  end if
2447  !
2448  ! -- evaluate magnitude of differences
2449  if (n == 1) then
2450  locdhmax = n
2451  dhmax = dh
2452  locrmax = n
2453  rmax = r
2454  dqfrommvrmax = dqfrommvr
2455  locdqfrommvrmax = n
2456  else
2457  if (abs(dh) > abs(dhmax)) then
2458  locdhmax = n
2459  dhmax = dh
2460  end if
2461  if (abs(r) > abs(rmax)) then
2462  locrmax = n
2463  rmax = r
2464  end if
2465  if (abs(dqfrommvr) > abs(dqfrommvrmax)) then
2466  dqfrommvrmax = dqfrommvr
2467  locdqfrommvrmax = n
2468  end if
2469  end if
2470  end do final_check
2471  !
2472  ! -- set dpak and cpak
2473  if (abs(dhmax) > abs(dpak)) then
2474  ipak = locdhmax
2475  dpak = dhmax
2476  write (cloc, "(a,'-',a)") trim(this%packName), 'stage'
2477  cpak = trim(cloc)
2478  end if
2479  if (abs(rmax) > abs(dpak)) then
2480  ipak = locrmax
2481  dpak = rmax
2482  write (cloc, "(a,'-',a)") trim(this%packName), 'inflow'
2483  cpak = trim(cloc)
2484  end if
2485  if (this%imover == 1) then
2486  if (abs(dqfrommvrmax) > abs(dpak)) then
2487  ipak = locdqfrommvrmax
2488  dpak = dqfrommvrmax
2489  write (cloc, "(a,'-',a)") trim(this%packName), 'qfrommvr'
2490  cpak = trim(cloc)
2491  end if
2492  end if
2493  !
2494  ! -- write convergence data to package csv
2495  if (this%ipakcsv /= 0) then
2496  !
2497  ! -- write the data
2498  call this%pakcsvtab%add_term(innertot)
2499  call this%pakcsvtab%add_term(totim)
2500  call this%pakcsvtab%add_term(kper)
2501  call this%pakcsvtab%add_term(kstp)
2502  call this%pakcsvtab%add_term(kiter)
2503  call this%pakcsvtab%add_term(dhmax)
2504  call this%pakcsvtab%add_term(locdhmax)
2505  call this%pakcsvtab%add_term(rmax)
2506  call this%pakcsvtab%add_term(locrmax)
2507  if (this%imover == 1) then
2508  call this%pakcsvtab%add_term(dqfrommvrmax)
2509  call this%pakcsvtab%add_term(locdqfrommvrmax)
2510  end if
2511  !
2512  ! -- finalize the package csv
2513  if (iend == 1) then
2514  call this%pakcsvtab%finalize_table()
2515  end if
2516  end if
2517  end if
real(dp), pointer, public totim
time relative to start of simulation
Definition: tdis.f90:32
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
Here is the call graph for this function:

◆ sfr_cf()

subroutine sfrmodule::sfr_cf ( class(sfrtype this)

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

Parameters
thisSfrType object

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

2146  ! -- dummy
2147  class(SfrType) :: this !< SfrType object
2148  ! -- local
2149  integer(I4B) :: n
2150  integer(I4B) :: igwfnode
2151  !
2152  ! -- return if no sfr reaches
2153  if (this%nbound == 0) return
2154  !
2155  ! -- find highest active cell
2156  do n = 1, this%nbound
2157  igwfnode = this%igwftopnode(n)
2158  if (igwfnode > 0) then
2159  if (this%ibound(igwfnode) == 0) then
2160  call this%dis%highest_active(igwfnode, this%ibound)
2161  end if
2162  end if
2163  this%igwfnode(n) = igwfnode
2164  this%nodelist(n) = igwfnode
2165  end do

◆ sfr_check_connections()

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

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

Parameters
thisSfrType object

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

4389  ! -- dummy
4390  class(SfrType) :: this !< SfrType object
4391  ! -- local
4392  logical(LGP) :: lreorder
4393  character(len=5) :: crch
4394  character(len=5) :: crch2
4395  character(len=LINELENGTH) :: text
4396  character(len=LINELENGTH) :: title
4397  integer(I4B) :: n
4398  integer(I4B) :: nn
4399  integer(I4B) :: nc
4400  integer(I4B) :: i
4401  integer(I4B) :: ii
4402  integer(I4B) :: j
4403  integer(I4B) :: ifound
4404  integer(I4B) :: ierr
4405  integer(I4B) :: maxconn
4406  integer(I4B) :: ntabcol
4407  !
4408  ! -- determine if the reaches have been reordered
4409  lreorder = .false.
4410  do j = 1, this%MAXBOUND
4411  n = this%isfrorder(j)
4412  if (n /= j) then
4413  lreorder = .true.
4414  exit
4415  end if
4416  end do
4417  !
4418  ! -- write message that the solution order h
4419  if (lreorder) then
4420  write (this%iout, '(/,1x,a)') &
4421  trim(adjustl(this%text))//' PACKAGE ('// &
4422  trim(adjustl(this%packName))//') REACH SOLUTION HAS BEEN '// &
4423  'REORDERED USING A DAG'
4424  !
4425  ! -- print table
4426  if (this%iprpak /= 0) then
4427  !
4428  ! -- reset the input table object
4429  ntabcol = 2
4430  title = trim(adjustl(this%text))//' PACKAGE ('// &
4431  trim(adjustl(this%packName))//') REACH SOLUTION ORDER'
4432  call table_cr(this%inputtab, this%packName, title)
4433  call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4434  text = 'ORDER'
4435  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4436  text = 'REACH'
4437  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4438  !
4439  ! -- upstream connection data
4440  do j = 1, this%maxbound
4441  n = this%isfrorder(j)
4442  call this%inputtab%add_term(j)
4443  call this%inputtab%add_term(n)
4444  end do
4445  end if
4446  end if
4447  !
4448  ! -- create input table for reach connections data
4449  if (this%iprpak /= 0) then
4450  !
4451  ! -- calculate the maximum number of connections
4452  maxconn = 0
4453  do n = 1, this%maxbound
4454  maxconn = max(maxconn, this%nconnreach(n))
4455  end do
4456  ntabcol = 1 + maxconn
4457  !
4458  ! -- reset the input table object
4459  title = trim(adjustl(this%text))//' PACKAGE ('// &
4460  trim(adjustl(this%packName))//') STATIC REACH CONNECTION DATA'
4461  call table_cr(this%inputtab, this%packName, title)
4462  call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4463  text = 'REACH'
4464  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4465  do n = 1, maxconn
4466  write (text, '(a,1x,i6)') 'CONN', n
4467  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4468  end do
4469  end if
4470  !
4471  ! -- check the reach connections for simple errors
4472  ! -- connection check
4473  do n = 1, this%MAXBOUND
4474  write (crch, '(i5)') n
4475  eachconn: do i = this%ia(n) + 1, this%ia(n + 1) - 1
4476  nn = this%ja(i)
4477  write (crch2, '(i5)') nn
4478  ifound = 0
4479  connreach: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4480  nc = this%ja(ii)
4481  if (nc == n) then
4482  ifound = 1
4483  exit connreach
4484  end if
4485  end do connreach
4486  if (ifound /= 1) then
4487  errmsg = 'Reach '//crch//' is connected to '// &
4488  'reach '//crch2//' but reach '//crch2// &
4489  ' is not connected to reach '//crch//'.'
4490  call store_error(errmsg)
4491  end if
4492  end do eachconn
4493  !
4494  ! -- write connection data to the table
4495  if (this%iprpak /= 0) then
4496  call this%inputtab%add_term(n)
4497  do i = this%ia(n) + 1, this%ia(n + 1) - 1
4498  call this%inputtab%add_term(this%ja(i))
4499  end do
4500  nn = maxconn - this%nconnreach(n)
4501  do i = 1, nn
4502  call this%inputtab%add_term(' ')
4503  end do
4504  end if
4505  end do
4506  !
4507  ! -- check for incorrect connections between upstream connections
4508  !
4509  ! -- check upstream connections for each reach
4510  ierr = 0
4511  do n = 1, this%maxbound
4512  write (crch, '(i5)') n
4513  eachconnv: do i = this%ia(n) + 1, this%ia(n + 1) - 1
4514  !
4515  ! -- skip downstream connections
4516  if (this%idir(i) < 0) cycle eachconnv
4517  nn = this%ja(i)
4518  write (crch2, '(i5)') nn
4519  connreachv: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4520  ! -- skip downstream connections
4521  if (this%idir(ii) < 0) cycle connreachv
4522  nc = this%ja(ii)
4523  !
4524  ! -- if nc == n then that means reach n is an upstream connection for
4525  ! reach nn and reach nn is an upstream connection for reach n
4526  if (nc == n) then
4527  ierr = ierr + 1
4528  errmsg = 'Reach '//crch//' is connected to '// &
4529  'reach '//crch2//' but streamflow from reach '// &
4530  crch//' to reach '//crch2//' is not permitted.'
4531  call store_error(errmsg)
4532  exit connreachv
4533  end if
4534  end do connreachv
4535  end do eachconnv
4536  end do
4537  !
4538  ! -- terminate if connectivity errors
4539  if (count_errors() > 0) then
4540  call this%parser%StoreErrorUnit()
4541  end if
4542  !
4543  ! -- check that downstream reaches for a reach are
4544  ! the upstream reaches for the reach
4545  do n = 1, this%maxbound
4546  write (crch, '(i5)') n
4547  eachconnds: do i = this%ia(n) + 1, this%ia(n + 1) - 1
4548  nn = this%ja(i)
4549  if (this%idir(i) > 0) cycle eachconnds
4550  write (crch2, '(i5)') nn
4551  ifound = 0
4552  connreachds: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4553  nc = this%ja(ii)
4554  if (nc == n) then
4555  if (this%idir(i) /= this%idir(ii)) then
4556  ifound = 1
4557  end if
4558  exit connreachds
4559  end if
4560  end do connreachds
4561  if (ifound /= 1) then
4562  errmsg = 'Reach '//crch//' downstream connected reach '// &
4563  'is reach '//crch2//' but reach '//crch//' is not'// &
4564  ' the upstream connected reach for reach '//crch2//'.'
4565  call store_error(errmsg)
4566  end if
4567  end do eachconnds
4568  end do
4569  !
4570  ! -- create input table for upstream and downstream connections
4571  if (this%iprpak /= 0) then
4572  !
4573  ! -- calculate the maximum number of upstream connections
4574  maxconn = 0
4575  do n = 1, this%maxbound
4576  ii = 0
4577  do i = this%ia(n) + 1, this%ia(n + 1) - 1
4578  if (this%idir(i) > 0) then
4579  ii = ii + 1
4580  end if
4581  end do
4582  maxconn = max(maxconn, ii)
4583  end do
4584  ntabcol = 1 + maxconn
4585  !
4586  ! -- reset the input table object
4587  title = trim(adjustl(this%text))//' PACKAGE ('// &
4588  trim(adjustl(this%packName))//') STATIC UPSTREAM REACH '// &
4589  'CONNECTION DATA'
4590  call table_cr(this%inputtab, this%packName, title)
4591  call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4592  text = 'REACH'
4593  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4594  do n = 1, maxconn
4595  write (text, '(a,1x,i6)') 'UPSTREAM CONN', n
4596  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4597  end do
4598  !
4599  ! -- upstream connection data
4600  do n = 1, this%maxbound
4601  call this%inputtab%add_term(n)
4602  ii = 0
4603  do i = this%ia(n) + 1, this%ia(n + 1) - 1
4604  if (this%idir(i) > 0) then
4605  call this%inputtab%add_term(this%ja(i))
4606  ii = ii + 1
4607  end if
4608  end do
4609  nn = maxconn - ii
4610  do i = 1, nn
4611  call this%inputtab%add_term(' ')
4612  end do
4613  end do
4614  !
4615  ! -- calculate the maximum number of downstream connections
4616  maxconn = 0
4617  do n = 1, this%maxbound
4618  ii = 0
4619  do i = this%ia(n) + 1, this%ia(n + 1) - 1
4620  if (this%idir(i) < 0) then
4621  ii = ii + 1
4622  end if
4623  end do
4624  maxconn = max(maxconn, ii)
4625  end do
4626  ntabcol = 1 + maxconn
4627  !
4628  ! -- reset the input table object
4629  title = trim(adjustl(this%text))//' PACKAGE ('// &
4630  trim(adjustl(this%packName))//') STATIC DOWNSTREAM '// &
4631  'REACH CONNECTION DATA'
4632  call table_cr(this%inputtab, this%packName, title)
4633  call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4634  text = 'REACH'
4635  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4636  do n = 1, maxconn
4637  write (text, '(a,1x,i6)') 'DOWNSTREAM CONN', n
4638  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4639  end do
4640  !
4641  ! -- downstream connection data
4642  do n = 1, this%maxbound
4643  call this%inputtab%add_term(n)
4644  ii = 0
4645  do i = this%ia(n) + 1, this%ia(n + 1) - 1
4646  if (this%idir(i) < 0) then
4647  call this%inputtab%add_term(this%ja(i))
4648  ii = ii + 1
4649  end if
4650  end do
4651  nn = maxconn - ii
4652  do i = 1, nn
4653  call this%inputtab%add_term(' ')
4654  end do
4655  end do
4656  end if
Here is the call graph for this function:

◆ sfr_check_conversion()

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

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

Parameters
thisSfrType object

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

4210  ! -- dummy
4211  class(SfrType) :: this !< SfrType object
4212  ! -- local
4213  ! -- formats
4214  character(len=*), parameter :: fmtunitconv_error = &
4215  &"('SFR (',a,') UNIT_CONVERSION SPECIFIED VALUE (',g0,') AND', &
4216  &1x,'LENGTH_CONVERSION OR TIME_CONVERSION SPECIFIED.')"
4217  character(len=*), parameter :: fmtunitconv = &
4218  &"(1x,'SFR PACKAGE (',a,') CONVERSION DATA',&
4219  &/4x,'UNIT CONVERSION VALUE (',g0,').',/)"
4220  !
4221  ! -- check the reach data for simple errors
4222  if (this%lengthconv /= dnodata .or. this%timeconv /= dnodata) then
4223  if (this%unitconv /= done) then
4224  write (errmsg, fmtunitconv_error) &
4225  trim(adjustl(this%packName)), this%unitconv
4226  call store_error(errmsg)
4227  else
4228  if (this%lengthconv /= dnodata) then
4229  this%unitconv = this%unitconv * this%lengthconv**donethird
4230  end if
4231  if (this%timeconv /= dnodata) then
4232  this%unitconv = this%unitconv * this%timeconv
4233  end if
4234  write (this%iout, fmtunitconv) &
4235  trim(adjustl(this%packName)), this%unitconv
4236  end if
4237  end if
Here is the call graph for this function:

◆ sfr_check_diversions()

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

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

Parameters
thisSfrType object

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

4666  ! -- dummy
4667  class(SfrType) :: this !< SfrType object
4668  ! -- local
4669  character(len=LINELENGTH) :: title
4670  character(len=LINELENGTH) :: text
4671  character(len=5) :: crch
4672  character(len=5) :: cdiv
4673  character(len=5) :: crch2
4674  character(len=10) :: cprior
4675  integer(I4B) :: maxdiv
4676  integer(I4B) :: n
4677  integer(I4B) :: nn
4678  integer(I4B) :: nc
4679  integer(I4B) :: ii
4680  integer(I4B) :: idiv
4681  integer(I4B) :: ifound
4682  integer(I4B) :: jpos
4683  ! -- format
4684 10 format('Diversion ', i0, ' of reach ', i0, &
4685  ' is invalid or has not been defined.')
4686  !
4687  ! -- write header
4688  if (this%iprpak /= 0) then
4689  !
4690  ! -- determine the maximum number of diversions
4691  maxdiv = 0
4692  do n = 1, this%maxbound
4693  maxdiv = maxdiv + this%ndiv(n)
4694  end do
4695  !
4696  ! -- reset the input table object
4697  title = trim(adjustl(this%text))//' PACKAGE ('// &
4698  trim(adjustl(this%packName))//') REACH DIVERSION DATA'
4699  call table_cr(this%inputtab, this%packName, title)
4700  call this%inputtab%table_df(maxdiv, 4, this%iout)
4701  text = 'REACH'
4702  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4703  text = 'DIVERSION'
4704  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4705  text = 'REACH 2'
4706  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4707  text = 'CPRIOR'
4708  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4709  end if
4710  !
4711  ! -- check that diversion data are correct
4712  do n = 1, this%maxbound
4713  if (this%ndiv(n) < 1) cycle
4714  write (crch, '(i5)') n
4715 
4716  do idiv = 1, this%ndiv(n)
4717  !
4718  ! -- determine diversion index
4719  jpos = this%iadiv(n) + idiv - 1
4720  !
4721  ! -- write idiv to cdiv
4722  write (cdiv, '(i5)') idiv
4723  !
4724  !
4725  nn = this%divreach(jpos)
4726  write (crch2, '(i5)') nn
4727  !
4728  ! -- make sure diversion reach is connected to current reach
4729  ifound = 0
4730  if (nn < 1 .or. nn > this%maxbound) then
4731  write (errmsg, 10) idiv, n
4732  call store_error(errmsg)
4733  cycle
4734  end if
4735  connreach: do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4736  nc = this%ja(ii)
4737  if (nc == n) then
4738  if (this%idir(ii) > 0) then
4739  ifound = 1
4740  end if
4741  exit connreach
4742  end if
4743  end do connreach
4744  if (ifound /= 1) then
4745  errmsg = 'Reach '//crch//' is not a upstream reach for '// &
4746  'reach '//crch2//' as a result diversion '//cdiv// &
4747  ' from reach '//crch//' to reach '//crch2// &
4748  ' is not possible. Check reach connectivity.'
4749  call store_error(errmsg)
4750  end if
4751  ! -- iprior
4752  cprior = this%divcprior(jpos)
4753  !
4754  ! -- add terms to the table
4755  if (this%iprpak /= 0) then
4756  call this%inputtab%add_term(n)
4757  call this%inputtab%add_term(idiv)
4758  call this%inputtab%add_term(nn)
4759  call this%inputtab%add_term(cprior)
4760  end if
4761  end do
4762  end do
Here is the call graph for this function:

◆ sfr_check_initialstages()

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

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

Parameters
thisSfrType object

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

4773  class(SfrType) :: this !< SfrType object
4774 
4775  character(len=LINELENGTH) :: title
4776  character(len=LINELENGTH) :: text
4777  character(len=5) :: crch
4778  integer(I4B) :: n
4779  real(DP) :: qman
4780 
4781  ! skip check if storage is not activated
4782  if (this%istorage == 0) return
4783 
4784  ! write header
4785  if (this%iprpak /= 0) then
4786  !
4787  ! -- reset the input table object
4788  title = trim(adjustl(this%text))//' PACKAGE ('// &
4789  trim(adjustl(this%packName))//') REACH INITIAL STAGE DATA'
4790  call table_cr(this%inputtab, this%packName, title)
4791  call this%inputtab%table_df(this%maxbound, 4, this%iout)
4792  text = 'REACH'
4793  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4794  text = 'INITIAL STAGE'
4795  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4796  text = 'INITIAL DEPTH'
4797  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4798  text = 'INITIAL FLOW'
4799  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4800  end if
4801  !
4802  ! -- check that data are correct
4803  do n = 1, this%maxbound
4804  write (crch, '(i5)') n
4805 
4806  ! calculate the initial flows
4807  call this%sfr_calc_qman(n, this%depth(n), qman)
4808  this%usinflow(n) = qman
4809  this%dsflow(n) = qman
4810 
4811  ! add terms to the table
4812  if (this%iprpak /= 0) then
4813  call this%inputtab%add_term(n)
4814  call this%inputtab%add_term(this%stage(n))
4815  call this%inputtab%add_term(this%depth(n))
4816  call this%inputtab%add_term(qman)
4817  end if
4818  end do
Here is the call graph for this function:

◆ sfr_check_reaches()

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

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

Parameters
thisSfrType object

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

4271  ! -- dummy
4272  class(SfrType) :: this !< SfrType object
4273  ! -- local
4274  character(len=5) :: crch
4275  character(len=10) :: cval
4276  character(len=30) :: nodestr
4277  character(len=LINELENGTH) :: title
4278  character(len=LINELENGTH) :: text
4279  integer(I4B) :: n
4280  integer(I4B) :: nn
4281  real(DP) :: btgwf
4282  real(DP) :: bt
4283  !
4284  ! -- setup inputtab tableobj
4285  if (this%iprpak /= 0) then
4286  title = trim(adjustl(this%text))//' PACKAGE ('// &
4287  trim(adjustl(this%packName))//') STATIC REACH DATA'
4288  call table_cr(this%inputtab, this%packName, title)
4289  call this%inputtab%table_df(this%maxbound, 10, this%iout)
4290  text = 'NUMBER'
4291  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4292  text = 'CELLID'
4293  call this%inputtab%initialize_column(text, 20, alignment=tableft)
4294  text = 'LENGTH'
4295  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4296  text = 'WIDTH'
4297  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4298  text = 'SLOPE'
4299  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4300  text = 'TOP'
4301  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4302  text = 'THICKNESS'
4303  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4304  text = 'HK'
4305  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4306  text = 'ROUGHNESS'
4307  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4308  text = 'UPSTREAM FRACTION'
4309  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4310  end if
4311  !
4312  ! -- check the reach data for simple errors
4313  do n = 1, this%maxbound
4314  write (crch, '(i5)') n
4315  nn = this%igwfnode(n)
4316  if (nn > 0) then
4317  btgwf = this%dis%bot(nn)
4318  call this%dis%noder_to_string(nn, nodestr)
4319  else
4320  nodestr = 'none'
4321  end if
4322  ! -- check reach length
4323  if (this%length(n) <= dzero) then
4324  errmsg = 'Reach '//crch//' length must be greater than 0.0.'
4325  call store_error(errmsg)
4326  end if
4327  ! -- check reach width
4328  if (this%width(n) <= dzero) then
4329  errmsg = 'Reach '//crch//' width must be greater than 0.0.'
4330  call store_error(errmsg)
4331  end if
4332  ! -- check reach slope
4333  if (this%slope(n) <= dzero) then
4334  errmsg = 'Reach '//crch//' slope must be greater than 0.0.'
4335  call store_error(errmsg)
4336  end if
4337  ! -- check bed thickness and bed hk for reaches connected to GWF
4338  if (nn > 0) then
4339  bt = this%strtop(n) - this%bthick(n)
4340  if (bt <= btgwf .and. this%icheck /= 0) then
4341  write (cval, '(f10.4)') bt
4342  errmsg = 'Reach '//crch//' bed bottom (rtp-rbth ='// &
4343  cval//') must be greater than the bottom of cell ('// &
4344  nodestr
4345  write (cval, '(f10.4)') btgwf
4346  errmsg = trim(adjustl(errmsg))//'='//cval//').'
4347  call store_error(errmsg)
4348  end if
4349  if (this%hk(n) < dzero) then
4350  errmsg = 'Reach '//crch//' hk must be greater than or equal to 0.0.'
4351  call store_error(errmsg)
4352  end if
4353  end if
4354  ! -- check reach roughness
4355  if (this%rough(n) <= dzero) then
4356  errmsg = 'Reach '//crch//" Manning's roughness "// &
4357  'coefficient must be greater than 0.0.'
4358  call store_error(errmsg)
4359  end if
4360  ! -- check reach upstream fraction
4361  if (this%ustrf(n) < dzero) then
4362  errmsg = 'Reach '//crch//' upstream fraction must be greater '// &
4363  'than or equal to 0.0.'
4364  call store_error(errmsg)
4365  end if
4366  ! -- write summary of reach information
4367  if (this%iprpak /= 0) then
4368  call this%inputtab%add_term(n)
4369  call this%inputtab%add_term(nodestr)
4370  call this%inputtab%add_term(this%length(n))
4371  call this%inputtab%add_term(this%width(n))
4372  call this%inputtab%add_term(this%slope(n))
4373  call this%inputtab%add_term(this%strtop(n))
4374  call this%inputtab%add_term(this%bthick(n))
4375  call this%inputtab%add_term(this%hk(n))
4376  call this%inputtab%add_term(this%rough(n))
4377  call this%inputtab%add_term(this%ustrf(n))
4378  end if
4379  end do
Here is the call graph for this function:

◆ sfr_check_storage_weight()

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

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

Parameters
thisSfrType object

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

4247  ! -- dummy
4248  class(SfrType) :: this !< SfrType object
4249  ! -- formats
4250  character(len=*), parameter :: fmtweight = &
4251  &"(1x,'SFR PACKAGE (',a,') SETTING DEFAULT',&
4252  &/4x,'STORAGE_WEIGHT VALUE (',g0,').',/)"
4253  !
4254  ! -- set storage weight if it has not been defined yet
4255  if (this%istorage == 1) then
4256  if (this%storage_weight == dnodata) then
4257  this%storage_weight = done
4258  write (this%iout, fmtweight) &
4259  trim(adjustl(this%packName)), this%storage_weight
4260  end if
4261  end if

◆ sfr_check_ustrf()

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

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

Parameters
thisSfrType object

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

4828  ! -- dummy
4829  class(SfrType) :: this !< SfrType object
4830  ! -- local
4831  character(len=LINELENGTH) :: title
4832  character(len=LINELENGTH) :: text
4833  logical(LGP) :: lcycle
4834  logical(LGP) :: ladd
4835  character(len=5) :: crch
4836  character(len=5) :: crch2
4837  character(len=10) :: cval
4838  integer(I4B) :: maxcols
4839  integer(I4B) :: npairs
4840  integer(I4B) :: ipair
4841  integer(I4B) :: i
4842  integer(I4B) :: n
4843  integer(I4B) :: n2
4844  integer(I4B) :: idiv
4845  integer(I4B) :: i0
4846  integer(I4B) :: i1
4847  integer(I4B) :: jpos
4848  integer(I4B) :: ids
4849  real(DP) :: f
4850  real(DP) :: rval
4851  !
4852  ! -- write table header
4853  if (this%iprpak /= 0) then
4854  !
4855  ! -- determine the maximum number of columns
4856  npairs = 0
4857  do n = 1, this%maxbound
4858  ipair = 0
4859  ec: do i = this%ia(n) + 1, this%ia(n + 1) - 1
4860  !
4861  ! -- skip upstream connections
4862  if (this%idir(i) > 0) cycle ec
4863  n2 = this%ja(i)
4864  !
4865  ! -- skip inactive downstream reaches
4866  if (this%iboundpak(n2) == 0) cycle ec
4867  !
4868  ! -- increment ipair and see if it exceeds npairs
4869  ipair = ipair + 1
4870  npairs = max(npairs, ipair)
4871  end do ec
4872  end do
4873  maxcols = 1 + npairs * 2
4874  !
4875  ! -- reset the input table object
4876  title = trim(adjustl(this%text))//' PACKAGE ('// &
4877  trim(adjustl(this%packName))//') CONNECTED REACH UPSTREAM '// &
4878  'FRACTION DATA'
4879  call table_cr(this%inputtab, this%packName, title)
4880  call this%inputtab%table_df(this%maxbound, maxcols, this%iout)
4881  text = 'REACH'
4882  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4883  do i = 1, npairs
4884  write (cval, '(i10)') i
4885  text = 'DOWNSTREAM REACH '//trim(adjustl(cval))
4886  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
4887  text = 'FRACTION '//trim(adjustl(cval))
4888  call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
4889  end do
4890  end if
4891  !
4892  ! -- fill diversion number for each connection
4893  do n = 1, this%maxbound
4894  do idiv = 1, this%ndiv(n)
4895  i0 = this%iadiv(n)
4896  i1 = this%iadiv(n + 1) - 1
4897  do jpos = i0, i1
4898  do i = this%ia(n) + 1, this%ia(n + 1) - 1
4899  n2 = this%ja(i)
4900  if (this%divreach(jpos) == n2) then
4901  this%idiv(i) = jpos - i0 + 1
4902  exit
4903  end if
4904  end do
4905  end do
4906  end do
4907  end do
4908  !
4909  ! -- check that the upstream fraction for reach connected by
4910  ! a diversion is zero
4911  do n = 1, this%maxbound
4912  !
4913  ! -- determine the number of downstream reaches
4914  ids = 0
4915  do i = this%ia(n) + 1, this%ia(n + 1) - 1
4916  if (this%idir(i) < 0) then
4917  ids = ids + 1
4918  end if
4919  end do
4920  !
4921  ! -- evaluate the diversions
4922  do idiv = 1, this%ndiv(n)
4923  jpos = this%iadiv(n) + idiv - 1
4924  n2 = this%divreach(jpos)
4925  f = this%ustrf(n2)
4926  if (f /= dzero) then
4927  write (errmsg, '(a,2(1x,i0,1x,a),1x,a,g0,a,2(1x,a))') &
4928  'Reach', n, 'is connected to reach', n2, 'by a diversion', &
4929  'but the upstream fraction is not equal to zero (', f, '). Check', &
4930  trim(this%packName), 'package diversion and package data.'
4931  if (ids > 1) then
4932  call store_error(errmsg)
4933  else
4934  write (warnmsg, '(a,3(1x,a))') &
4935  trim(warnmsg), &
4936  'A warning instead of an error is issued because', &
4937  'the reach is only connected to the diversion reach in the ', &
4938  'downstream direction.'
4939  call store_warning(warnmsg)
4940  end if
4941  end if
4942  end do
4943  end do
4944  !
4945  ! -- calculate the total fraction of connected reaches that are
4946  ! not diversions and check that the sum of upstream fractions
4947  ! is equal to 1 for each reach
4948  do n = 1, this%maxbound
4949  ids = 0
4950  rval = dzero
4951  f = dzero
4952  write (crch, '(i5)') n
4953  if (this%iprpak /= 0) then
4954  call this%inputtab%add_term(n)
4955  end if
4956  ipair = 0
4957  eachconn: do i = this%ia(n) + 1, this%ia(n + 1) - 1
4958  lcycle = .false.
4959  !
4960  ! -- initialize downstream connection q
4961  this%qconn(i) = dzero
4962  !
4963  ! -- skip upstream connections
4964  if (this%idir(i) > 0) then
4965  lcycle = .true.
4966  end if
4967  n2 = this%ja(i)
4968  !
4969  ! -- skip inactive downstream reaches
4970  if (this%iboundpak(n2) == 0) then
4971  lcycle = .true.
4972  end if
4973  if (lcycle) then
4974  cycle eachconn
4975  end if
4976  ipair = ipair + 1
4977  write (crch2, '(i5)') n2
4978  ids = ids + 1
4979  ladd = .true.
4980  f = f + this%ustrf(n2)
4981  write (cval, '(f10.4)') this%ustrf(n2)
4982  !
4983  ! -- write upstream fractions
4984  if (this%iprpak /= 0) then
4985  call this%inputtab%add_term(n2)
4986  call this%inputtab%add_term(this%ustrf(n2))
4987  end if
4988  eachdiv: do idiv = 1, this%ndiv(n)
4989  jpos = this%iadiv(n) + idiv - 1
4990  if (this%divreach(jpos) == n2) then
4991  ladd = .false.
4992  exit eachdiv
4993  end if
4994  end do eachdiv
4995  if (ladd) then
4996  rval = rval + this%ustrf(n2)
4997  end if
4998  end do eachconn
4999  this%ftotnd(n) = rval
5000  !
5001  ! -- write remaining table columns
5002  if (this%iprpak /= 0) then
5003  ipair = ipair + 1
5004  do i = ipair, npairs
5005  call this%inputtab%add_term(' ')
5006  call this%inputtab%add_term(' ')
5007  end do
5008  end if
5009  !
5010  ! -- evaluate if an error condition has occurred
5011  ! the sum of fractions is not equal to 1
5012  if (ids /= 0) then
5013  if (abs(f - done) > dem6) then
5014  write (errmsg, '(a,1x,i0,1x,a,g0,a,3(1x,a))') &
5015  'Upstream fractions for reach ', n, 'is not equal to one (', f, &
5016  '). Check', trim(this%packName), 'package reach connectivity and', &
5017  'package data.'
5018  call store_error(errmsg)
5019  end if
5020  end if
5021  end do
Here is the call graph for this function:

◆ sfr_cq()

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

Calculate the flow between connected SFR package control volumes.

Parameters
[in,out]thisSfrType object
[in]xcurrent dependent-variable value
[in,out]flowjaflow between two connected control volumes
[in]iadvflag that indicates if this is an advance package

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

2525  ! -- modules
2526  use budgetmodule, only: budgettype
2527  ! -- dummy
2528  class(SfrType), intent(inout) :: this !< SfrType object
2529  real(DP), dimension(:), intent(in) :: x !< current dependent-variable value
2530  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
2531  integer(I4B), optional, intent(in) :: iadv !< flag that indicates if this is an advance package
2532  ! -- local
2533  integer(I4B) :: i
2534  real(DP) :: qext
2535  ! -- for budget
2536  integer(I4B) :: n
2537  integer(I4B) :: n2
2538  real(DP) :: qoutflow
2539  real(DP) :: qfrommvr
2540  real(DP) :: qtomvr
2541  !
2542  ! -- call base functionality in bnd_cq. This will calculate sfr-gwf flows
2543  ! and put them into this%simvals
2544  call this%BndType%bnd_cq(x, flowja, iadv=1)
2545  !
2546  ! -- Calculate qextoutflow and qoutflow for subsequent budgets
2547  do n = 1, this%maxbound
2548  !
2549  ! -- mover
2550  qfrommvr = dzero
2551  qtomvr = dzero
2552  if (this%imover == 1) then
2553  qfrommvr = this%pakmvrobj%get_qfrommvr(n)
2554  qtomvr = this%pakmvrobj%get_qtomvr(n)
2555  if (qtomvr > dzero) then
2556  qtomvr = -qtomvr
2557  end if
2558  end if
2559  !
2560  ! -- external downstream stream flow
2561  qext = this%dsflow(n)
2562  qoutflow = dzero
2563  if (qext > dzero) then
2564  qext = -qext
2565  end if
2566  do i = this%ia(n) + 1, this%ia(n + 1) - 1
2567  if (this%idir(i) > 0) cycle
2568  n2 = this%ja(i)
2569  if (this%iboundpak(n2) == 0) cycle
2570  qext = dzero
2571  exit
2572  end do
2573  !
2574  ! -- adjust external downstream stream flow using qtomvr
2575  if (qext < dzero) then
2576  if (qtomvr < dzero) then
2577  qext = qext - qtomvr
2578  end if
2579  else
2580  qoutflow = this%dsflow(n)
2581  if (qoutflow > dzero) then
2582  qoutflow = -qoutflow
2583  end if
2584  end if
2585  !
2586  ! -- set qextoutflow and qoutflow for cell by cell budget
2587  ! output and observations
2588  this%qextoutflow(n) = qext
2589  this%qoutflow(n) = qoutflow
2590  !
2591  end do
2592  !
2593  ! -- fill the budget object
2594  call this%sfr_fill_budobj()
This module contains the BudgetModule.
Definition: Budget.f90:20
Derived type for the Budget object.
Definition: Budget.f90:39

◆ sfr_create()

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

Create a new SFR Package object

Parameters
packobjpointer to default package type
[in]idpackage id
[in]ibcnumboundary condition number
[in]inunitunit number of SFR package input file
[in]ioutunit number of model listing file
[in]namemodelmodel name
[in]paknamepackage name

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

295  ! -- modules
297  ! -- dummy
298  class(BndType), pointer :: packobj !< pointer to default package type
299  integer(I4B), intent(in) :: id !< package id
300  integer(I4B), intent(in) :: ibcnum !< boundary condition number
301  integer(I4B), intent(in) :: inunit !< unit number of SFR package input file
302  integer(I4B), intent(in) :: iout !< unit number of model listing file
303  character(len=*), intent(in) :: namemodel !< model name
304  character(len=*), intent(in) :: pakname !< package name
305  ! -- local
306  type(SfrType), pointer :: sfrobj
307  !
308  ! -- allocate the object and assign values to object variables
309  allocate (sfrobj)
310  packobj => sfrobj
311  !
312  ! -- create name and memory path
313  call packobj%set_names(ibcnum, namemodel, pakname, ftype)
314  packobj%text = text
315  !
316  ! -- allocate scalars
317  call sfrobj%sfr_allocate_scalars()
318  !
319  ! -- initialize package
320  call packobj%pack_initialize()
321 
322  packobj%inunit = inunit
323  packobj%iout = iout
324  packobj%id = id
325  packobj%ibcnum = ibcnum
326  packobj%ncolbnd = 4
327  packobj%iscloc = 0 ! not supported
328  packobj%isadvpak = 1
329  packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
Here is the call graph for this function:
Here is the caller graph for this function:

◆ sfr_da()

subroutine sfrmodule::sfr_da ( class(sfrtype this)

Deallocate SFR package scalars and arrays.

Parameters
thisSfrType object

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

2784  ! -- modules
2786  ! -- dummy
2787  class(SfrType) :: this !< SfrType object
2788  !
2789  ! -- deallocate arrays
2790  call mem_deallocate(this%qoutflow)
2791  call mem_deallocate(this%qextoutflow)
2792  deallocate (this%csfrbudget)
2793  call mem_deallocate(this%sfrname, 'SFRNAME', this%memoryPath)
2794  call mem_deallocate(this%dbuff)
2795  deallocate (this%cauxcbc)
2796  call mem_deallocate(this%qauxcbc)
2797  call mem_deallocate(this%iboundpak)
2798  call mem_deallocate(this%igwfnode)
2799  call mem_deallocate(this%igwftopnode)
2800  call mem_deallocate(this%length)
2801  call mem_deallocate(this%width)
2802  call mem_deallocate(this%strtop)
2803  call mem_deallocate(this%bthick)
2804  call mem_deallocate(this%hk)
2805  call mem_deallocate(this%slope)
2806  call mem_deallocate(this%nconnreach)
2807  call mem_deallocate(this%ustrf)
2808  call mem_deallocate(this%ftotnd)
2809  call mem_deallocate(this%usflow)
2810  call mem_deallocate(this%dsflow)
2811  call mem_deallocate(this%depth)
2812  call mem_deallocate(this%stage)
2813  call mem_deallocate(this%gwflow)
2814  call mem_deallocate(this%simevap)
2815  call mem_deallocate(this%simrunoff)
2816  call mem_deallocate(this%stage0)
2817  call mem_deallocate(this%usflow0)
2818  call mem_deallocate(this%denseterms)
2819  call mem_deallocate(this%viscratios)
2820  !
2821  ! -- stage, usflow, and dsflow for previous timestep
2822  if (this%istorage == 1) then
2823  call mem_deallocate(this%stageold)
2824  call mem_deallocate(this%dsflowold)
2825  call mem_deallocate(this%storage)
2826  call mem_deallocate(this%usinflow)
2827  call mem_deallocate(this%usinflowold)
2828  end if
2829  !
2830  ! -- deallocate reach order and connection data
2831  call mem_deallocate(this%isfrorder)
2832  call mem_deallocate(this%ia)
2833  call mem_deallocate(this%ja)
2834  call mem_deallocate(this%idir)
2835  call mem_deallocate(this%idiv)
2836  call mem_deallocate(this%qconn)
2837  !
2838  ! -- deallocate boundary data
2839  call mem_deallocate(this%rough)
2840  call mem_deallocate(this%rain)
2841  call mem_deallocate(this%evap)
2842  call mem_deallocate(this%inflow)
2843  call mem_deallocate(this%runoff)
2844  call mem_deallocate(this%sstage)
2845  !
2846  ! -- deallocate aux variables
2847  call mem_deallocate(this%rauxvar)
2848  !
2849  ! -- deallocate diversion variables
2850  call mem_deallocate(this%iadiv)
2851  call mem_deallocate(this%divreach)
2852  if (associated(this%divcprior)) then
2853  deallocate (this%divcprior)
2854  end if
2855  call mem_deallocate(this%divflow)
2856  call mem_deallocate(this%divq)
2857  call mem_deallocate(this%ndiv)
2858  !
2859  ! -- deallocate cross-section data
2860  call mem_deallocate(this%ncrosspts)
2861  call mem_deallocate(this%iacross)
2862  call mem_deallocate(this%station)
2863  call mem_deallocate(this%xsheight)
2864  call mem_deallocate(this%xsrough)
2865  !
2866  ! -- deallocate budobj
2867  call this%budobj%budgetobject_da()
2868  deallocate (this%budobj)
2869  nullify (this%budobj)
2870  !
2871  ! -- deallocate stage table
2872  if (this%iprhed > 0) then
2873  call this%stagetab%table_da()
2874  deallocate (this%stagetab)
2875  nullify (this%stagetab)
2876  end if
2877  !
2878  ! -- deallocate package csv table
2879  if (this%ipakcsv > 0) then
2880  if (associated(this%pakcsvtab)) then
2881  call this%pakcsvtab%table_da()
2882  deallocate (this%pakcsvtab)
2883  nullify (this%pakcsvtab)
2884  end if
2885  end if
2886  !
2887  ! -- deallocate scalars
2888  call mem_deallocate(this%istorage)
2889  call mem_deallocate(this%storage_weight)
2890  call mem_deallocate(this%iprhed)
2891  call mem_deallocate(this%istageout)
2892  call mem_deallocate(this%ibudgetout)
2893  call mem_deallocate(this%ibudcsv)
2894  call mem_deallocate(this%ipakcsv)
2895  call mem_deallocate(this%idiversions)
2896  call mem_deallocate(this%maxsfrpicard)
2897  call mem_deallocate(this%maxsfrit)
2898  call mem_deallocate(this%bditems)
2899  call mem_deallocate(this%cbcauxitems)
2900  call mem_deallocate(this%unitconv)
2901  call mem_deallocate(this%lengthconv)
2902  call mem_deallocate(this%timeconv)
2903  call mem_deallocate(this%dmaxchg)
2904  call mem_deallocate(this%deps)
2905  call mem_deallocate(this%nconn)
2906  call mem_deallocate(this%icheck)
2907  call mem_deallocate(this%iconvchk)
2908  call mem_deallocate(this%idense)
2909  call mem_deallocate(this%ianynone)
2910  call mem_deallocate(this%ncrossptstot)
2911  nullify (this%gwfiss)
2912  !
2913  ! -- call base BndType deallocate
2914  call this%BndType%bnd_da()

◆ sfr_df_obs()

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

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

Parameters
thisSfrType object

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

2967  ! -- dummy
2968  class(SfrType) :: this !< SfrType object
2969  ! -- local
2970  integer(I4B) :: indx
2971  !
2972  ! -- Store obs type and assign procedure pointer
2973  ! for stage observation type.
2974  call this%obs%StoreObsType('stage', .false., indx)
2975  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
2976  !
2977  ! -- Store obs type and assign procedure pointer
2978  ! for inflow observation type.
2979  call this%obs%StoreObsType('inflow', .true., indx)
2980  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
2981  !
2982  ! -- Store obs type and assign procedure pointer
2983  ! for inflow observation type.
2984  call this%obs%StoreObsType('ext-inflow', .true., indx)
2985  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
2986  !
2987  ! -- Store obs type and assign procedure pointer
2988  ! for rainfall observation type.
2989  call this%obs%StoreObsType('rainfall', .true., indx)
2990  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
2991  !
2992  ! -- Store obs type and assign procedure pointer
2993  ! for runoff observation type.
2994  call this%obs%StoreObsType('runoff', .true., indx)
2995  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
2996  !
2997  ! -- Store obs type and assign procedure pointer
2998  ! for evaporation observation type.
2999  call this%obs%StoreObsType('evaporation', .true., indx)
3000  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3001  !
3002  ! -- Store obs type and assign procedure pointer
3003  ! for outflow observation type.
3004  call this%obs%StoreObsType('outflow', .true., indx)
3005  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3006  !
3007  ! -- Store obs type and assign procedure pointer
3008  ! for ext-outflow observation type.
3009  call this%obs%StoreObsType('ext-outflow', .true., indx)
3010  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3011  !
3012  ! -- Store obs type and assign procedure pointer
3013  ! for to-mvr observation type.
3014  call this%obs%StoreObsType('to-mvr', .true., indx)
3015  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3016  !
3017  ! -- Store obs type and assign procedure pointer
3018  ! for sfr-frommvr observation type.
3019  call this%obs%StoreObsType('from-mvr', .true., indx)
3020  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3021  !
3022  ! -- Store obs type and assign procedure pointer
3023  ! for sfr observation type.
3024  call this%obs%StoreObsType('sfr', .true., indx)
3025  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3026  !
3027  ! -- Store obs type and assign procedure pointer
3028  ! for upstream flow observation type.
3029  call this%obs%StoreObsType('upstream-flow', .true., indx)
3030  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3031  !
3032  ! -- Store obs type and assign procedure pointer
3033  ! for downstream flow observation type.
3034  call this%obs%StoreObsType('downstream-flow', .true., indx)
3035  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3036  !
3037  ! -- Store obs type and assign procedure pointer
3038  ! for depth observation type.
3039  call this%obs%StoreObsType('depth', .false., indx)
3040  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3041  !
3042  ! -- Store obs type and assign procedure pointer
3043  ! for wetted-perimeter observation type.
3044  call this%obs%StoreObsType('wet-perimeter', .false., indx)
3045  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3046  !
3047  ! -- Store obs type and assign procedure pointer
3048  ! for wetted-area observation type.
3049  call this%obs%StoreObsType('wet-area', .false., indx)
3050  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid
3051  !
3052  ! -- Store obs type and assign procedure pointer
3053  ! for wetted-width observation type.
3054  call this%obs%StoreObsType('wet-width', .false., indx)
3055  this%obs%obsData(indx)%ProcessIdPtr => sfr_process_obsid

◆ sfr_fc()

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

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

Parameters
thisSfrType object
[in,out]rhsright-hand side vector for model
[in]iasolution CRS row pointers
[in]idxglomapping vector for model (local) to solution (global)
matrix_slnsolution coefficient matrix

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

2174  ! -- dummy
2175  class(SfrType) :: this !< SfrType object
2176  real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model
2177  integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers
2178  integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global)
2179  class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
2180  ! -- local
2181  integer(I4B) :: i
2182  integer(I4B) :: j
2183  integer(I4B) :: n
2184  integer(I4B) :: ipos
2185  integer(I4B) :: node
2186  real(DP) :: s0
2187  real(DP) :: ds
2188  real(DP) :: dsmax
2189  real(DP) :: hgwf
2190  real(DP) :: v
2191  real(DP) :: hhcof
2192  real(DP) :: rrhs
2193  !
2194  ! -- picard iterations for sfr to achieve good solution regardless
2195  ! of reach order
2196  sfrpicard: do i = 1, this%maxsfrpicard
2197  !
2198  ! -- initialize maximum stage change for iteration to zero
2199  dsmax = dzero
2200  !
2201  ! -- pakmvrobj fc - reset qformvr to zero
2202  if (this%imover == 1) then
2203  call this%pakmvrobj%fc()
2204  end if
2205  !
2206  ! -- solve for each sfr reach
2207  reachsolve: do j = 1, this%nbound
2208  n = this%isfrorder(j)
2209  node = this%igwfnode(n)
2210  if (node > 0) then
2211  hgwf = this%xnew(node)
2212  else
2213  hgwf = dep20
2214  end if
2215  !
2216  ! -- save previous stage and upstream flow
2217  if (i == 1) then
2218  this%stage0(n) = this%stage(n)
2219  this%usflow0(n) = this%usflow(n)
2220  end if
2221  !
2222  ! -- set initial stage to calculate stage change
2223  s0 = this%stage(n)
2224  !
2225  ! -- solve for flow in swr
2226  if (this%iboundpak(n) /= 0) then
2227  call this%sfr_solve(n, hgwf, hhcof, rrhs)
2228  else
2229  this%depth(n) = dzero
2230  this%stage(n) = this%strtop(n)
2231  v = dzero
2232  call this%sfr_update_flows(n, v, v)
2233  hhcof = dzero
2234  rrhs = dzero
2235  end if
2236  !
2237  ! -- set package hcof and rhs
2238  this%hcof(n) = hhcof
2239  this%rhs(n) = rrhs
2240  !
2241  ! -- calculate stage change
2242  ds = s0 - this%stage(n)
2243  !
2244  ! -- evaluate if stage change exceeds dsmax
2245  if (abs(ds) > abs(dsmax)) then
2246  dsmax = ds
2247  end if
2248 
2249  end do reachsolve
2250  !
2251  ! -- evaluate if the sfr picard iterations should be terminated
2252  if (abs(dsmax) <= this%dmaxchg) then
2253  exit sfrpicard
2254  end if
2255 
2256  end do sfrpicard
2257  !
2258  ! -- Copy package rhs and hcof into solution rhs and amat
2259  do n = 1, this%nbound
2260  node = this%nodelist(n)
2261  if (node < 1) cycle
2262  rhs(node) = rhs(node) + this%rhs(n)
2263  ipos = ia(node)
2264  call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(n))
2265  end do

◆ sfr_fill_budobj()

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

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

Parameters
thisSfrType object

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

5246  ! -- dummy
5247  class(SfrType) :: this !< SfrType object
5248  ! -- local
5249  integer(I4B) :: naux
5250  integer(I4B) :: i
5251  integer(I4B) :: n
5252  integer(I4B) :: n1
5253  integer(I4B) :: n2
5254  integer(I4B) :: ii
5255  integer(I4B) :: idx
5256  integer(I4B) :: idiv
5257  integer(I4B) :: jpos
5258  real(DP) :: q
5259  real(DP) :: qt
5260  real(DP) :: d
5261  real(DP) :: ca
5262  real(DP) :: a
5263  real(DP) :: wp
5264  real(DP) :: l
5265  !
5266  ! -- initialize counter
5267  idx = 0
5268  !
5269  ! -- FLOW JA FACE
5270  idx = idx + 1
5271  call this%budobj%budterm(idx)%reset(this%nconn)
5272  do n = 1, this%maxbound
5273  n1 = n
5274  q = dzero
5275  ca = dzero
5276  do i = this%ia(n) + 1, this%ia(n + 1) - 1
5277  n2 = this%ja(i)
5278  if (this%iboundpak(n) /= 0) then
5279  ! flow to downstream reaches
5280  if (this%idir(i) < 0) then
5281  qt = this%dsflow(n)
5282  q = -this%qconn(i)
5283  ! flow from upstream reaches
5284  else
5285  qt = this%usflow(n)
5286  do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
5287  if (this%idir(ii) > 0) cycle
5288  if (this%ja(ii) /= n) cycle
5289  q = this%qconn(ii)
5290  exit
5291  end do
5292  end if
5293  ! calculate flow area
5294  call this%sfr_calc_reach_depth(n, qt, d)
5295  ca = this%calc_area_wet(n, d)
5296  else
5297  q = dzero
5298  ca = dzero
5299  end if
5300  this%qauxcbc(1) = ca
5301  call this%budobj%budterm(idx)%update_term(n1, n2, q, this%qauxcbc)
5302  end do
5303  end do
5304  !
5305  ! -- GWF (LEAKAGE)
5306  idx = idx + 1
5307  call this%budobj%budterm(idx)%reset(this%maxbound - this%ianynone)
5308  do n = 1, this%maxbound
5309  n2 = this%igwfnode(n)
5310  if (n2 > 0) then
5311  if (this%iboundpak(n) /= 0) then
5312  ! -- calc_perimeter_wet() does not enforce depth dependence
5313  if (this%depth(n) > dzero) then
5314  wp = this%calc_perimeter_wet(n, this%depth(n))
5315  else
5316  wp = dzero
5317  end if
5318  l = this%length(n)
5319  a = wp * l
5320  this%qauxcbc(1) = a
5321  q = -this%gwflow(n)
5322  else
5323  this%qauxcbc(1) = dzero
5324  q = dzero
5325  end if
5326  call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
5327  end if
5328  end do
5329  !
5330  ! -- RAIN
5331  idx = idx + 1
5332  call this%budobj%budterm(idx)%reset(this%maxbound)
5333  do n = 1, this%maxbound
5334  if (this%iboundpak(n) /= 0) then
5335  a = this%calc_surface_area(n)
5336  q = this%rain(n) * a
5337  else
5338  q = dzero
5339  end if
5340  call this%budobj%budterm(idx)%update_term(n, n, q)
5341  end do
5342  !
5343  ! -- EVAPORATION
5344  idx = idx + 1
5345  call this%budobj%budterm(idx)%reset(this%maxbound)
5346  do n = 1, this%maxbound
5347  if (this%iboundpak(n) /= 0) then
5348  q = -this%simevap(n)
5349  else
5350  q = dzero
5351  end if
5352  call this%budobj%budterm(idx)%update_term(n, n, q)
5353  end do
5354  !
5355  ! -- RUNOFF
5356  idx = idx + 1
5357  call this%budobj%budterm(idx)%reset(this%maxbound)
5358  do n = 1, this%maxbound
5359  if (this%iboundpak(n) /= 0) then
5360  q = this%simrunoff(n)
5361  else
5362  q = dzero
5363  end if
5364  call this%budobj%budterm(idx)%update_term(n, n, q)
5365  end do
5366  !
5367  ! -- INFLOW
5368  idx = idx + 1
5369  call this%budobj%budterm(idx)%reset(this%maxbound)
5370  do n = 1, this%maxbound
5371  if (this%iboundpak(n) /= 0) then
5372  q = this%inflow(n)
5373  else
5374  q = dzero
5375  end if
5376  call this%budobj%budterm(idx)%update_term(n, n, q)
5377  end do
5378  !
5379  ! -- EXTERNAL OUTFLOW
5380  idx = idx + 1
5381  call this%budobj%budterm(idx)%reset(this%maxbound)
5382  do n = 1, this%maxbound
5383  q = dzero
5384  if (this%iboundpak(n) /= 0) then
5385  do i = this%ia(n) + 1, this%ia(n + 1) - 1
5386  if (this%idir(i) > 0) cycle
5387  idiv = this%idiv(i)
5388  if (idiv > 0) then
5389  jpos = this%iadiv(n) + idiv - 1
5390  q = q + this%divq(jpos)
5391  else
5392  q = q + this%qconn(i)
5393  end if
5394  end do
5395  q = q - this%dsflow(n)
5396  if (this%imover == 1) then
5397  q = q + this%pakmvrobj%get_qtomvr(n)
5398  end if
5399  else
5400  if (this%imover == 1) then
5401  q = this%pakmvrobj%get_qfrommvr(n)
5402  end if
5403  end if
5404  call this%budobj%budterm(idx)%update_term(n, n, q)
5405  end do
5406  !
5407  ! -- STORAGE
5408  idx = idx + 1
5409  call this%budobj%budterm(idx)%reset(this%maxbound)
5410  do n = 1, this%maxbound
5411  q = dzero
5412  if (this%iboundpak(n) /= 0) then
5413  d = this%depth(n)
5414  a = this%calc_surface_area_wet(n, d)
5415  this%qauxcbc(1) = a * d
5416  if (this%gwfiss == 0 .and. this%istorage == 1) then
5417  q = this%storage(n)
5418  end if
5419  else
5420  q = dzero
5421  this%qauxcbc(1) = dzero
5422  end if
5423  call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc)
5424  end do
5425  !
5426  ! -- MOVER
5427  if (this%imover == 1) then
5428  !
5429  ! -- FROM MOVER
5430  idx = idx + 1
5431  call this%budobj%budterm(idx)%reset(this%maxbound)
5432  do n = 1, this%maxbound
5433  q = dzero
5434  if (this%iboundpak(n) /= 0) then
5435  q = this%pakmvrobj%get_qfrommvr(n)
5436  end if
5437  call this%budobj%budterm(idx)%update_term(n, n, q)
5438  end do
5439  !
5440  ! -- TO MOVER
5441  idx = idx + 1
5442  call this%budobj%budterm(idx)%reset(this%maxbound)
5443  do n = 1, this%maxbound
5444  if (this%iboundpak(n) /= 0) then
5445  q = this%pakmvrobj%get_qtomvr(n)
5446  if (q > dzero) then
5447  q = -q
5448  end if
5449  else
5450  q = dzero
5451  end if
5452  call this%budobj%budterm(idx)%update_term(n, n, q)
5453  end do
5454  end if
5455  !
5456  ! -- AUXILIARY VARIABLES
5457  naux = this%naux
5458  if (naux > 0) then
5459  idx = idx + 1
5460  call this%budobj%budterm(idx)%reset(this%maxbound)
5461  do n = 1, this%maxbound
5462  q = dzero
5463  call this%budobj%budterm(idx)%update_term(n, n, q, this%auxvar(:, n))
5464  end do
5465  end if
5466  !
5467  ! --Terms are filled, now accumulate them for this time step
5468  call this%budobj%accumulate_terms()

◆ sfr_fn()

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

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

Parameters
thisSfrType object
[in,out]rhsright-hand side vector for model
[in]iasolution CRS row pointers
[in]idxglomapping vector for model (local) to solution (global)
matrix_slnsolution coefficient matrix

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

2274  ! -- dummy
2275  class(SfrType) :: this !< SfrType object
2276  real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model
2277  integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers
2278  integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global)
2279  class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
2280  ! -- local
2281  integer(I4B) :: i
2282  integer(I4B) :: j
2283  integer(I4B) :: n
2284  integer(I4B) :: ipos
2285  real(DP) :: rterm
2286  real(DP) :: drterm
2287  real(DP) :: rhs1
2288  real(DP) :: hcof1
2289  real(DP) :: q1
2290  real(DP) :: q2
2291  real(DP) :: hgwf
2292  !
2293  ! -- Copy package rhs and hcof into solution rhs and amat
2294  do j = 1, this%nbound
2295  i = this%isfrorder(j)
2296  ! -- skip inactive reaches
2297  if (this%iboundpak(i) < 1) cycle
2298  ! -- skip if reach is not connected to gwf
2299  n = this%nodelist(i)
2300  if (n < 1) cycle
2301  ipos = ia(n)
2302  rterm = this%hcof(i) * this%xnew(n)
2303  ! -- calculate perturbed head
2304  hgwf = this%xnew(n) + dem4
2305  call this%sfr_solve(i, hgwf, hcof1, rhs1, update=.false.)
2306  q1 = rhs1 - hcof1 * hgwf
2307  ! -- calculate unperturbed head
2308  q2 = this%rhs(i) - this%hcof(i) * this%xnew(n)
2309  ! -- calculate derivative
2310  drterm = (q2 - q1) / dem4
2311  ! -- add terms to convert conductance formulation into
2312  ! newton-raphson formulation
2313  call matrix_sln%add_value_pos(idxglo(ipos), drterm - this%hcof(i))
2314  rhs(n) = rhs(n) - rterm + drterm * this%xnew(n)
2315  end do

◆ sfr_gwf_conn()

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

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

Returns
flag indicating if reach is connected to a gwf cell
Parameters
thisSfrType object
[in]nreach number

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

4006  ! -- return variable
4007  integer(I4B) :: sfr_gwf_conn !< flag indicating if reach is connected to a gwf cell
4008  ! -- dummy
4009  class(SfrType) :: this !< SfrType object
4010  integer(I4B), intent(in) :: n !< reach number
4011  ! -- local
4012  integer(I4B) :: node
4013 
4014  sfr_gwf_conn = 0
4015  node = this%igwfnode(n)
4016  if (node > 0 .and. this%hk(n) > dzero) then
4017  sfr_gwf_conn = 1
4018  end if

◆ sfr_obs_supported()

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

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

Returns
sfr_obs_supported boolean indicating if observations are supported
Parameters
thisSfrType object

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

2955  ! -- dummy
2956  class(SfrType) :: this !< SfrType object
2957  !
2958  ! -- set boolean
2959  sfr_obs_supported = .true.

◆ sfr_options()

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

Read additional options for SFR package.

Parameters
[in,out]thisSfrType object
[in,out]optionoption keyword string
[in,out]foundboolean indicating if option found

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

703  ! -- modules
704  use openspecmodule, only: access, form
706  ! -- dummy
707  class(SfrType), intent(inout) :: this !< SfrType object
708  character(len=*), intent(inout) :: option !< option keyword string
709  logical(LGP), intent(inout) :: found !< boolean indicating if option found
710  ! -- local
711  real(DP) :: r
712  character(len=MAXCHARLEN) :: fname
713  character(len=MAXCHARLEN) :: keyword
714  ! -- formats
715  character(len=*), parameter :: fmttimeconv = &
716  &"(4x, 'TIME CONVERSION VALUE (',g0,') SPECIFIED.')"
717  character(len=*), parameter :: fmtlengthconv = &
718  &"(4x, 'LENGTH CONVERSION VALUE (',g0,') SPECIFIED.')"
719  character(len=*), parameter :: fmtpicard = &
720  &"(4x, 'MAXIMUM SFR PICARD ITERATION VALUE (',i0,') SPECIFIED.')"
721  character(len=*), parameter :: fmtiter = &
722  &"(4x, 'MAXIMUM SFR ITERATION VALUE (',i0,') SPECIFIED.')"
723  character(len=*), parameter :: fmtdmaxchg = &
724  &"(4x, 'MAXIMUM DEPTH CHANGE VALUE (',g0,') SPECIFIED.')"
725  character(len=*), parameter :: fmtsfrbin = &
726  "(4x, 'SFR ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, /4x, &
727  &'OPENED ON UNIT: ', I0)"
728  character(len=*), parameter :: fmtstoweight = &
729  &"(4x, 'KINEMATIC STORAGE WEIGHT (',g0,') SPECIFIED.')"
730  !
731  ! -- Check for SFR options
732  found = .true.
733  select case (option)
734  case ('STORAGE')
735  this%istorage = 1
736  write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
737  ' REACH STORAGE IS ACTIVE.'
738  case ('PRINT_STAGE')
739  this%iprhed = 1
740  write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
741  ' STAGES WILL BE PRINTED TO LISTING FILE.'
742  case ('STAGE')
743  call this%parser%GetStringCaps(keyword)
744  if (keyword == 'FILEOUT') then
745  call this%parser%GetString(fname)
746  this%istageout = getunit()
747  call openfile(this%istageout, this%iout, fname, 'DATA(BINARY)', &
748  form, access, 'REPLACE', mnormal)
749  write (this%iout, fmtsfrbin) &
750  'STAGE', trim(adjustl(fname)), this%istageout
751  else
752  call store_error('Optional stage keyword must &
753  &be followed by fileout.')
754  end if
755  case ('BUDGET')
756  call this%parser%GetStringCaps(keyword)
757  if (keyword == 'FILEOUT') then
758  call this%parser%GetString(fname)
759  call assign_iounit(this%ibudgetout, this%inunit, "BUDGET fileout")
760  call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
761  form, access, 'REPLACE', mnormal)
762  write (this%iout, fmtsfrbin) &
763  'BUDGET', trim(adjustl(fname)), this%ibudgetout
764  else
765  call store_error('Optional budget keyword must be '// &
766  'followed by fileout.')
767  end if
768  case ('BUDGETCSV')
769  call this%parser%GetStringCaps(keyword)
770  if (keyword == 'FILEOUT') then
771  call this%parser%GetString(fname)
772  call assign_iounit(this%ibudcsv, this%inunit, "BUDGETCSV fileout")
773  call openfile(this%ibudcsv, this%iout, fname, 'CSV', &
774  filstat_opt='REPLACE')
775  write (this%iout, fmtsfrbin) &
776  'BUDGET CSV', trim(adjustl(fname)), this%ibudcsv
777  else
778  call store_error('OPTIONAL BUDGETCSV KEYWORD MUST BE FOLLOWED BY &
779  &FILEOUT')
780  end if
781  case ('PACKAGE_CONVERGENCE')
782  call this%parser%GetStringCaps(keyword)
783  if (keyword == 'FILEOUT') then
784  call this%parser%GetString(fname)
785  this%ipakcsv = getunit()
786  call openfile(this%ipakcsv, this%iout, fname, 'CSV', &
787  filstat_opt='REPLACE', mode_opt=mnormal)
788  write (this%iout, fmtsfrbin) &
789  'PACKAGE_CONVERGENCE', trim(adjustl(fname)), this%ipakcsv
790  else
791  call store_error('Optional package_convergence keyword must be '// &
792  'followed by fileout.')
793  end if
794  case ('UNIT_CONVERSION')
795  this%unitconv = this%parser%GetDouble()
796  !
797  ! -- create warning message
798  write (warnmsg, '(a)') &
799  'SETTING UNIT_CONVERSION DIRECTLY'
800  !
801  ! -- create deprecation warning
802  call deprecation_warning('OPTIONS', 'UNIT_CONVERSION', '6.4.2', &
803  warnmsg, this%parser%GetUnit())
804  case ('LENGTH_CONVERSION')
805  this%lengthconv = this%parser%GetDouble()
806  write (this%iout, fmtlengthconv) this%lengthconv
807  case ('TIME_CONVERSION')
808  this%timeconv = this%parser%GetDouble()
809  write (this%iout, fmttimeconv) this%timeconv
810  case ('MAXIMUM_PICARD_ITERATIONS')
811  this%maxsfrpicard = this%parser%GetInteger()
812  write (this%iout, fmtpicard) this%maxsfrpicard
813  case ('MAXIMUM_ITERATIONS')
814  this%maxsfrit = this%parser%GetInteger()
815  write (this%iout, fmtiter) this%maxsfrit
816  case ('MAXIMUM_DEPTH_CHANGE')
817  r = this%parser%GetDouble()
818  this%dmaxchg = r
819  this%deps = dp999 * r
820  write (this%iout, fmtdmaxchg) this%dmaxchg
821  case ('MOVER')
822  this%imover = 1
823  write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
824  !
825  ! -- right now these are options that are only available in the
826  ! development version and are not included in the documentation.
827  ! These options are only available when IDEVELOPMODE in
828  ! constants module is set to 1
829  case ('DEV_NO_CHECK')
830  call this%parser%DevOpt()
831  this%icheck = 0
832  write (this%iout, '(4x,A)') 'SFR CHECKS OF REACH GEOMETRY '// &
833  'RELATIVE TO MODEL GRID AND '// &
834  'REASONABLE PARAMETERS WILL NOT '// &
835  'BE PERFORMED.'
836  case ('DEV_NO_FINAL_CHECK')
837  call this%parser%DevOpt()
838  this%iconvchk = 0
839  write (this%iout, '(4x,a)') &
840  'A FINAL CONVERGENCE CHECK OF THE CHANGE IN STREAM FLOW ROUTING &
841  &STAGES AND FLOWS WILL NOT BE MADE'
842  case ('DEV_STORAGE_WEIGHT')
843  r = this%parser%GetDouble()
844  if (r < dhalf .or. r > done) then
845  write (errmsg, '(a,g0,a)') &
846  "STORAGE_WEIGHT SPECIFIED TO BE '", r, &
847  "' BUT CANNOT BE LESS THAN 0.5 OR GREATER THAN 1.0"
848  call store_error(errmsg)
849  else
850  this%storage_weight = r
851  write (this%iout, fmtstoweight) this%storage_weight
852  end if
853  !
854  ! -- no valid options found
855  case default
856  !
857  ! -- No options found
858  found = .false.
859  end select
subroutine, public assign_iounit(iounit, errunit, description)
@ brief assign io unit number
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
Here is the call graph for this function:

◆ sfr_ot_bdsummary()

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

Output SFR package budget summary.

Parameters
thisSfrType object
[in]kstptime step number
[in]kperperiod number
[in]ioutflag and unit number for the model listing file
[in]ibudflflag indicating budget should be written

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

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

◆ sfr_ot_dv()

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

Output SFR boundary package dependent-variable terms.

Parameters
thisSfrType object
[in]idvsaveflag and unit number for dependent-variable output
[in]idvprintflag indicating if dependent-variable should be written to the model listing file

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

2655  ! -- modules
2656  use tdismodule, only: kstp, kper, pertim, totim
2657  use inputoutputmodule, only: ulasav
2658  ! -- dummy
2659  class(SfrType) :: this !< SfrType object
2660  integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output
2661  integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file
2662  ! -- local
2663  character(len=20) :: cellid
2664  integer(I4B) :: ibinun
2665  integer(I4B) :: n
2666  integer(I4B) :: node
2667  real(DP) :: d
2668  real(DP) :: v
2669  real(DP) :: hgwf
2670  real(DP) :: sbot
2671  real(DP) :: depth
2672  real(DP) :: stage
2673  real(DP) :: w
2674  real(DP) :: cond
2675  real(DP) :: grad
2676  !
2677  ! -- set unit number for binary dependent variable output
2678  ibinun = 0
2679  if (this%istageout /= 0) then
2680  ibinun = this%istageout
2681  end if
2682  if (idvsave == 0) ibinun = 0
2683  !
2684  ! -- write sfr binary output
2685  if (ibinun > 0) then
2686  do n = 1, this%maxbound
2687  d = this%depth(n)
2688  v = this%stage(n)
2689  if (this%iboundpak(n) == 0) then
2690  v = dhnoflo
2691  else if (d == dzero) then
2692  v = dhdry
2693  end if
2694  this%dbuff(n) = v
2695  end do
2696  call ulasav(this%dbuff, ' STAGE', kstp, kper, pertim, totim, &
2697  this%maxbound, 1, 1, ibinun)
2698  end if
2699  !
2700  ! -- print sfr stage and depth table
2701  if (idvprint /= 0 .and. this%iprhed /= 0) then
2702  !
2703  ! -- set table kstp and kper
2704  call this%stagetab%set_kstpkper(kstp, kper)
2705  !
2706  ! -- fill stage data
2707  do n = 1, this%maxbound
2708  node = this%igwfnode(n)
2709  if (node > 0) then
2710  call this%dis%noder_to_string(node, cellid)
2711  hgwf = this%xnew(node)
2712  else
2713  cellid = 'NONE'
2714  end if
2715  if (this%inamedbound == 1) then
2716  call this%stagetab%add_term(this%boundname(n))
2717  end if
2718  call this%stagetab%add_term(n)
2719  call this%stagetab%add_term(cellid)
2720  if (this%iboundpak(n) /= 0) then
2721  depth = this%depth(n)
2722  stage = this%stage(n)
2723  w = this%calc_top_width_wet(n, depth)
2724  call this%sfr_calc_cond(n, depth, cond, stage, hgwf)
2725  else
2726  depth = dhnoflo
2727  stage = dhnoflo
2728  w = dhnoflo
2729  cond = dhnoflo
2730  end if
2731  if (depth == dzero) then
2732  call this%stagetab%add_term(dhdry)
2733  else
2734  call this%stagetab%add_term(stage)
2735  end if
2736  call this%stagetab%add_term(depth)
2737  call this%stagetab%add_term(w)
2738  if (node > 0) then
2739  if (this%iboundpak(n) /= 0) then
2740  sbot = this%strtop(n) - this%bthick(n)
2741  if (hgwf < sbot) then
2742  grad = stage - sbot
2743  else
2744  grad = stage - hgwf
2745  end if
2746  grad = grad / this%bthick(n)
2747  else
2748  grad = dhnoflo
2749  end if
2750  call this%stagetab%add_term(hgwf)
2751  call this%stagetab%add_term(cond)
2752  call this%stagetab%add_term(grad)
2753  else
2754  call this%stagetab%add_term('--')
2755  call this%stagetab%add_term('--')
2756  call this%stagetab%add_term('--')
2757  end if
2758  end do
2759  end if
subroutine, public ulasav(buf, text, kstp, kper, pertim, totim, ncol, nrow, ilay, ichn)
Save 1 layer array on disk.
real(dp), pointer, public pertim
time relative to start of stress period
Definition: tdis.f90:30
Here is the call graph for this function:

◆ sfr_ot_package_flows()

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

Output SFR package flow terms.

Parameters
thisSfrType object
[in]icbcflflag and unit number for cell-by-cell output
[in]ibudflflag indication if cell-by-cell data should be saved

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

2602  ! -- modules
2603  use tdismodule, only: kstp, kper, delt, pertim, totim
2604  ! -- dummy
2605  class(SfrType) :: this !< SfrType object
2606  integer(I4B), intent(in) :: icbcfl !< flag and unit number for cell-by-cell output
2607  integer(I4B), intent(in) :: ibudfl !< flag indication if cell-by-cell data should be saved
2608  ! -- local
2609  integer(I4B) :: ibinun
2610  character(len=20), dimension(:), allocatable :: cellidstr
2611  integer(I4B) :: n
2612  integer(I4B) :: node
2613  !
2614  ! -- write the flows from the budobj
2615  ibinun = 0
2616  if (this%ibudgetout /= 0) then
2617  ibinun = this%ibudgetout
2618  end if
2619  if (icbcfl == 0) ibinun = 0
2620  if (ibinun > 0) then
2621  call this%budobj%save_flows(this%dis, ibinun, kstp, kper, delt, &
2622  pertim, totim, this%iout)
2623  end if
2624  !
2625  ! -- Print sfr flows table
2626  if (ibudfl /= 0 .and. this%iprflow /= 0) then
2627  !
2628  ! -- If there are any 'none' gwf connections then need to calculate
2629  ! a vector of cellids and pass that in to the budget flow table because
2630  ! the table assumes that there are maxbound gwf entries, which is not
2631  ! the case if any 'none's are specified.
2632  if (this%ianynone > 0) then
2633  allocate (cellidstr(this%maxbound))
2634  do n = 1, this%maxbound
2635  node = this%igwfnode(n)
2636  if (node > 0) then
2637  call this%dis%noder_to_string(node, cellidstr(n))
2638  else
2639  cellidstr(n) = 'NONE'
2640  end if
2641  end do
2642  call this%budobj%write_flowtable(this%dis, kstp, kper, cellidstr)
2643  deallocate (cellidstr)
2644  else
2645  call this%budobj%write_flowtable(this%dis, kstp, kper)
2646  end if
2647  end if

◆ sfr_process_obsid()

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

Method to process observation ID strings for a SFR package.

Parameters
[in,out]obsrvObservation object
[in]disDiscretization object
[in]inunitobsfile unit number for the package observation file
[in]ioutmodel listing file unit number

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

3265  ! -- dummy
3266  type(ObserveType), intent(inout) :: obsrv !< Observation object
3267  class(DisBaseType), intent(in) :: dis !< Discretization object
3268  integer(I4B), intent(in) :: inunitobs !< file unit number for the package observation file
3269  integer(I4B), intent(in) :: iout !< model listing file unit number
3270  ! -- local
3271  integer(I4B) :: nn1
3272  integer(I4B) :: icol
3273  integer(I4B) :: istart
3274  integer(I4B) :: istop
3275  character(len=LINELENGTH) :: string
3276  character(len=LENBOUNDNAME) :: bndname
3277  !
3278  ! -- initialize local variables
3279  string = obsrv%IDstring
3280  !
3281  ! -- Extract reach number from string and store it.
3282  ! If 1st item is not an integer(I4B), it should be a
3283  ! boundary name--deal with it.
3284  icol = 1
3285  !
3286  ! -- get reach number or boundary name
3287  call extract_idnum_or_bndname(string, icol, istart, istop, nn1, bndname)
3288  if (nn1 == namedboundflag) then
3289  obsrv%FeatureName = bndname
3290  end if
3291  !
3292  ! -- store reach number (NodeNumber)
3293  obsrv%NodeNumber = nn1
Here is the call graph for this function:

◆ sfr_read_connectiondata()

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

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

Parameters
[in,out]thisSfrType object

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

1322  ! -- modules
1324  use sparsemodule, only: sparsematrix
1325  ! -- dummy
1326  class(SfrType), intent(inout) :: this !< SfrType object
1327  ! -- local
1328  character(len=LINELENGTH) :: line
1329  logical(LGP) :: isfound
1330  logical(LGP) :: endOfBlock
1331  integer(I4B) :: n
1332  integer(I4B) :: i
1333  integer(I4B) :: j
1334  integer(I4B) :: jj
1335  integer(I4B) :: jcol
1336  integer(I4B) :: jcol2
1337  integer(I4B) :: nja
1338  integer(I4B) :: ival
1339  integer(I4B) :: idir
1340  integer(I4B) :: ierr
1341  integer(I4B) :: nconnmax
1342  integer(I4B) :: nup
1343  integer(I4B) :: ipos
1344  integer(I4B) :: istat
1345  integer(I4B), dimension(:), pointer, contiguous :: rowmaxnnz => null()
1346  integer, allocatable, dimension(:) :: nboundchk
1347  integer, allocatable, dimension(:, :) :: iconndata
1348  type(sparsematrix), pointer :: sparse => null()
1349  integer(I4B), dimension(:), allocatable :: iup
1350  integer(I4B), dimension(:), allocatable :: order
1351  type(dag) :: sfr_dag
1352  !
1353  ! -- allocate and initialize local variables for reach connections
1354  allocate (nboundchk(this%maxbound))
1355  do n = 1, this%maxbound
1356  nboundchk(n) = 0
1357  end do
1358  !
1359  ! -- calculate the number of non-zero entries (size of ja maxtrix)
1360  nja = 0
1361  nconnmax = 0
1362  allocate (rowmaxnnz(this%maxbound))
1363  do n = 1, this%maxbound
1364  ival = this%nconnreach(n)
1365  if (ival < 0) ival = 0
1366  rowmaxnnz(n) = ival + 1
1367  nja = nja + ival + 1
1368  if (ival > nconnmax) then
1369  nconnmax = ival
1370  end if
1371  end do
1372  !
1373  ! -- reallocate connection data for package
1374  call mem_reallocate(this%ja, nja, 'JA', this%memoryPath)
1375  call mem_reallocate(this%idir, nja, 'IDIR', this%memoryPath)
1376  call mem_reallocate(this%idiv, nja, 'IDIV', this%memoryPath)
1377  call mem_reallocate(this%qconn, nja, 'QCONN', this%memoryPath)
1378  !
1379  ! -- initialize connection data
1380  do n = 1, nja
1381  this%idir(n) = 0
1382  this%idiv(n) = 0
1383  this%qconn(n) = dzero
1384  end do
1385  !
1386  ! -- allocate space for iconndata
1387  allocate (iconndata(nconnmax, this%maxbound))
1388  !
1389  ! -- initialize iconndata
1390  do n = 1, this%maxbound
1391  do j = 1, nconnmax
1392  iconndata(j, n) = 0
1393  end do
1394  end do
1395  !
1396  ! -- allocate space for connectivity
1397  allocate (sparse)
1398  !
1399  ! -- set up sparse
1400  call sparse%init(this%maxbound, this%maxbound, rowmaxnnz)
1401  !
1402  ! -- read connection data
1403  call this%parser%GetBlock('CONNECTIONDATA', isfound, ierr, &
1404  supportopenclose=.true.)
1405  !
1406  ! -- parse reach connectivity block if detected
1407  if (isfound) then
1408  write (this%iout, '(/1x,a)') &
1409  'PROCESSING '//trim(adjustl(this%text))//' CONNECTIONDATA'
1410  do
1411  call this%parser%GetNextLine(endofblock)
1412  if (endofblock) exit
1413  !
1414  ! -- get reach number
1415  n = this%parser%GetInteger()
1416  !
1417  ! -- check for error
1418  if (n < 1 .or. n > this%maxbound) then
1419  write (errmsg, '(a,1x,a,1x,i0)') &
1420  'SFR reach in connectiondata block is less than one or greater', &
1421  'than NREACHES:', n
1422  call store_error(errmsg)
1423  cycle
1424  end if
1425  !
1426  ! -- increment nboundchk
1427  nboundchk(n) = nboundchk(n) + 1
1428  !
1429  ! -- add diagonal connection for reach
1430  call sparse%addconnection(n, n, 1)
1431  !
1432  ! -- fill off diagonals
1433  do i = 1, this%nconnreach(n)
1434  !
1435  ! -- get connected reach
1436  ival = this%parser%GetInteger()
1437  !
1438  ! -- save connection data to temporary iconndata
1439  iconndata(i, n) = ival
1440  !
1441  ! -- determine idir
1442  if (ival < 0) then
1443  idir = -1
1444  ival = abs(ival)
1445  elseif (ival == 0) then
1446  call store_error('Missing or zero connection reach in line:')
1447  call store_error(line)
1448  else
1449  idir = 1
1450  end if
1451  if (ival > this%maxbound) then
1452  call store_error('Reach number exceeds NREACHES in line:')
1453  call store_error(line)
1454  end if
1455  !
1456  ! -- add connection to sparse
1457  call sparse%addconnection(n, ival, 1)
1458  end do
1459  end do
1460 
1461  write (this%iout, '(1x,a)') &
1462  'END OF '//trim(adjustl(this%text))//' CONNECTIONDATA'
1463 
1464  do n = 1, this%maxbound
1465  !
1466  ! -- check for missing or duplicate sfr connections
1467  if (nboundchk(n) == 0) then
1468  write (errmsg, '(a,1x,i0)') &
1469  'No connection data specified for reach', n
1470  call store_error(errmsg)
1471  else if (nboundchk(n) > 1) then
1472  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1473  'Connection data for reach', n, &
1474  'specified', nboundchk(n), 'times.'
1475  call store_error(errmsg)
1476  end if
1477  end do
1478  else
1479  call store_error('Required connectiondata block not found.')
1480  end if
1481  !
1482  ! -- terminate if errors encountered in connectiondata block
1483  if (count_errors() > 0) then
1484  call this%parser%StoreErrorUnit()
1485  end if
1486  !
1487  ! -- create ia and ja from sparse
1488  call sparse%filliaja(this%ia, this%ja, ierr, sort=.true.)
1489  !
1490  ! -- test for error condition
1491  if (ierr /= 0) then
1492  write (errmsg, '(a,3(1x,a))') &
1493  'Could not fill', trim(this%packName), &
1494  'package IA and JA connection data.', &
1495  'Check connectivity data in connectiondata block.'
1496  call store_error(errmsg)
1497  end if
1498  !
1499  ! -- fill flat connection storage
1500  do n = 1, this%maxbound
1501  do j = this%ia(n) + 1, this%ia(n + 1) - 1
1502  jcol = this%ja(j)
1503  do jj = 1, this%nconnreach(n)
1504  jcol2 = iconndata(jj, n)
1505  if (abs(jcol2) == jcol) then
1506  idir = 1
1507  if (jcol2 < 0) then
1508  idir = -1
1509  end if
1510  this%idir(j) = idir
1511  exit
1512  end if
1513  end do
1514  end do
1515  end do
1516  !
1517  ! -- deallocate temporary local storage for reach connections
1518  deallocate (rowmaxnnz)
1519  deallocate (nboundchk)
1520  deallocate (iconndata)
1521  !
1522  ! -- destroy sparse
1523  call sparse%destroy()
1524  deallocate (sparse)
1525  !
1526  ! -- calculate reach order using DAG
1527  !
1528  ! -- initialize the DAG
1529  call sfr_dag%set_vertices(this%maxbound)
1530  !
1531  ! -- fill DAG
1532  fill_dag: do n = 1, this%maxbound
1533  !
1534  ! -- determine the number of upstream reaches
1535  nup = 0
1536  do j = this%ia(n) + 1, this%ia(n + 1) - 1
1537  if (this%idir(j) > 0) then
1538  nup = nup + 1
1539  end if
1540  end do
1541  !
1542  ! -- cycle if nu upstream reacches
1543  if (nup == 0) cycle fill_dag
1544  !
1545  ! -- allocate local storage
1546  allocate (iup(nup))
1547  !
1548  ! -- fill local storage
1549  ipos = 1
1550  do j = this%ia(n) + 1, this%ia(n + 1) - 1
1551  if (this%idir(j) > 0) then
1552  iup(ipos) = this%ja(j)
1553  ipos = ipos + 1
1554  end if
1555  end do
1556  !
1557  ! -- add upstream connections to DAG
1558  call sfr_dag%set_edges(n, iup)
1559  !
1560  ! -- clean up local storage
1561  deallocate (iup)
1562  end do fill_dag
1563  !
1564  ! -- perform toposort on DAG
1565  call sfr_dag%toposort(order, istat)
1566  !
1567  ! -- write warning if circular dependency
1568  if (istat == -1) then
1569  write (warnmsg, '(a)') &
1570  trim(adjustl(this%text))//' PACKAGE ('// &
1571  trim(adjustl(this%packName))//') cannot calculate a '// &
1572  'Directed Asyclic Graph for reach connectivity because '// &
1573  'of circular dependency. Using the reach number for '// &
1574  'solution ordering.'
1575  call store_warning(warnmsg)
1576  end if
1577  !
1578  ! -- fill isfrorder
1579  do n = 1, this%maxbound
1580  if (istat == 0) then
1581  this%isfrorder(n) = order(n)
1582  else
1583  this%isfrorder(n) = n
1584  end if
1585  end do
1586  !
1587  ! -- clean up DAG and remaining local storage
1588  call sfr_dag%destroy()
1589  if (istat == 0) then
1590  deallocate (order)
1591  end if
Here is the call graph for this function:

◆ sfr_read_crossection()

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

Method to read crosssection data for the SFR package.

Parameters
[in,out]thisSfrType object

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

1184  ! -- modules
1187  ! -- dummy
1188  class(SfrType), intent(inout) :: this !< SfrType object
1189  ! -- local
1190  character(len=LINELENGTH) :: keyword
1191  character(len=LINELENGTH) :: line
1192  logical(LGP) :: isfound
1193  logical(LGP) :: endOfBlock
1194  integer(I4B) :: n
1195  integer(I4B) :: ierr
1196  integer(I4B) :: ncrossptstot
1197  integer, allocatable, dimension(:) :: nboundchk
1198  type(SfrCrossSection), pointer :: cross_data => null()
1199  !
1200  ! -- read cross-section data
1201  call this%parser%GetBlock('CROSSSECTIONS', isfound, ierr, &
1202  supportopenclose=.true., &
1203  blockrequired=.false.)
1204  !
1205  ! -- parse reach connectivity block if detected
1206  if (isfound) then
1207  write (this%iout, '(/1x,a)') &
1208  'PROCESSING '//trim(adjustl(this%text))//' CROSSSECTIONS'
1209  !
1210  ! -- allocate and initialize local variables for reach cross-sections
1211  allocate (nboundchk(this%maxbound))
1212  do n = 1, this%maxbound
1213  nboundchk(n) = 0
1214  end do
1215  !
1216  ! -- create and initialize cross-section data
1217  call cross_section_cr(cross_data, this%iout, this%iprpak, this%maxbound)
1218  call cross_data%initialize(this%ncrossptstot, this%ncrosspts, &
1219  this%iacross, &
1220  this%station, this%xsheight, &
1221  this%xsrough)
1222  !
1223  ! -- read all of the entries in the block
1224  readtable: do
1225  call this%parser%GetNextLine(endofblock)
1226  if (endofblock) exit
1227  !
1228  ! -- get reach number
1229  n = this%parser%GetInteger()
1230  !
1231  ! -- check for reach number error
1232  if (n < 1 .or. n > this%maxbound) then
1233  write (errmsg, '(a,1x,a,1x,i0)') &
1234  'SFR reach in crosssections block is less than one or greater', &
1235  'than NREACHES:', n
1236  call store_error(errmsg)
1237  cycle readtable
1238  end if
1239  !
1240  ! -- increment nboundchk
1241  nboundchk(n) = nboundchk(n) + 1
1242  !
1243  ! -- read FILE keyword
1244  call this%parser%GetStringCaps(keyword)
1245  select case (keyword)
1246  case ('TAB6')
1247  call this%parser%GetStringCaps(keyword)
1248  if (trim(adjustl(keyword)) /= 'FILEIN') then
1249  errmsg = 'TAB6 keyword must be followed by "FILEIN" '// &
1250  'then by filename.'
1251  call store_error(errmsg)
1252  cycle readtable
1253  end if
1254  call this%parser%GetString(line)
1255  call cross_data%read_table(n, this%width(n), &
1256  trim(adjustl(line)))
1257  case default
1258  write (errmsg, '(a,1x,i4,1x,a)') &
1259  'CROSS-SECTION TABLE ENTRY for REACH ', n, &
1260  'MUST INCLUDE TAB6 KEYWORD'
1261  call store_error(errmsg)
1262  cycle readtable
1263  end select
1264  end do readtable
1265 
1266  write (this%iout, '(1x,a)') &
1267  'END OF '//trim(adjustl(this%text))//' CROSSSECTIONS'
1268 
1269  !
1270  ! -- check for duplicate sfr crosssections
1271  do n = 1, this%maxbound
1272  if (nboundchk(n) > 1) then
1273  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1274  'Cross-section data for reach', n, &
1275  'specified', nboundchk(n), 'times.'
1276  call store_error(errmsg)
1277  end if
1278  end do
1279  !
1280  ! -- terminate if errors encountered in cross-sections block
1281  if (count_errors() > 0) then
1282  call this%parser%StoreErrorUnit()
1283  end if
1284  !
1285  ! -- determine the current size of cross-section data
1286  ncrossptstot = cross_data%get_ncrossptstot()
1287  !
1288  ! -- reallocate sfr package cross-section data
1289  if (ncrossptstot /= this%ncrossptstot) then
1290  this%ncrossptstot = ncrossptstot
1291  call mem_reallocate(this%station, this%ncrossptstot, 'STATION', &
1292  this%memoryPath)
1293  call mem_reallocate(this%xsheight, this%ncrossptstot, 'XSHEIGHT', &
1294  this%memoryPath)
1295  call mem_reallocate(this%xsrough, this%ncrossptstot, 'XSROUGH', &
1296  this%memoryPath)
1297  end if
1298  !
1299  ! -- write cross-section data to the model listing file
1300  call cross_data%output(this%width, this%rough)
1301  !
1302  ! -- pack cross-section data
1303  call cross_data%pack(this%ncrossptstot, this%ncrosspts, &
1304  this%iacross, &
1305  this%station, &
1306  this%xsheight, &
1307  this%xsrough)
1308  !
1309  ! -- deallocate temporary local storage for reach cross-sections
1310  deallocate (nboundchk)
1311  call cross_data%destroy()
1312  deallocate (cross_data)
1313  nullify (cross_data)
1314  end if
subroutine, public cross_section_cr(this, iout, iprpak, nreaches)
Create a cross-section object.
Here is the call graph for this function:

◆ sfr_read_dimensions()

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

Read dimensions for the SFR package.

Parameters
[in,out]thisSfrType object

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

615  ! -- dummy
616  class(SfrType), intent(inout) :: this !< SfrType object
617  ! -- local
618  character(len=LINELENGTH) :: keyword
619  integer(I4B) :: ierr
620  logical(LGP) :: isfound
621  logical(LGP) :: endOfBlock
622  !
623  ! -- initialize dimensions to 0
624  this%maxbound = 0
625  !
626  ! -- get dimensions block
627  call this%parser%GetBlock('DIMENSIONS', isfound, ierr, &
628  supportopenclose=.true.)
629  !
630  ! -- parse dimensions block if detected
631  if (isfound) then
632  write (this%iout, '(/1x,a)') &
633  'PROCESSING '//trim(adjustl(this%text))//' DIMENSIONS'
634  do
635  call this%parser%GetNextLine(endofblock)
636  if (endofblock) exit
637  call this%parser%GetStringCaps(keyword)
638  select case (keyword)
639  case ('NREACHES')
640  this%maxbound = this%parser%GetInteger()
641  write (this%iout, '(4x,a,i0)') 'NREACHES = ', this%maxbound
642  case default
643  write (errmsg, '(2a)') &
644  'Unknown '//trim(this%text)//' dimension: ', trim(keyword)
645  call store_error(errmsg)
646  end select
647  end do
648  write (this%iout, '(1x,a)') &
649  'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
650  else
651  call store_error('Required dimensions block not found.')
652  end if
653  !
654  ! -- verify dimensions were set
655  if (this%maxbound < 1) then
656  write (errmsg, '(a)') &
657  'NREACHES was not specified or was specified incorrectly.'
658  call store_error(errmsg)
659  end if
660  !
661  ! -- write summary of error messages for block
662  if (count_errors() > 0) then
663  call this%parser%StoreErrorUnit()
664  end if
665  !
666  ! -- Call define_listlabel to construct the list label that is written
667  ! when PRINT_INPUT option is used.
668  call this%define_listlabel()
669  !
670  ! -- Define default cross-section data size
671  this%ncrossptstot = this%maxbound
672  !
673  ! -- Allocate arrays in package superclass
674  call this%sfr_allocate_arrays()
675  !
676  ! -- read package data
677  call this%sfr_read_packagedata()
678  !
679  ! -- read cross-section data
680  call this%sfr_read_crossection()
681  !
682  ! -- read connection data
683  call this%sfr_read_connectiondata()
684  !
685  ! -- read diversion data
686  call this%sfr_read_diversions()
687  !
688  ! -- read initial stage data
689  call this%sfr_read_initial_stages()
690  !
691  ! -- setup the budget object
692  call this%sfr_setup_budobj()
693  !
694  ! -- setup the stage table object
695  call this%sfr_setup_tableobj()
Here is the call graph for this function:

◆ sfr_read_diversions()

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

Method to read diversions for the SFR package.

Parameters
[in,out]thisSfrType object

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

1599  ! -- modules
1601  ! -- dummy
1602  class(SfrType), intent(inout) :: this !< SfrType object
1603  ! -- local
1604  character(len=10) :: cnum
1605  character(len=10) :: cval
1606  integer(I4B) :: j
1607  integer(I4B) :: n
1608  integer(I4B) :: ierr
1609  integer(I4B) :: ival
1610  integer(I4B) :: i0
1611  integer(I4B) :: ipos
1612  integer(I4B) :: jpos
1613  integer(I4B) :: ndiv
1614  integer(I4B) :: ndiversions
1615  integer(I4B) :: idivreach
1616  logical(LGP) :: isfound
1617  logical(LGP) :: endOfBlock
1618  integer(I4B) :: idiv
1619  integer, allocatable, dimension(:) :: iachk
1620  integer, allocatable, dimension(:) :: nboundchk
1621  !
1622  ! -- determine the total number of diversions and fill iadiv
1623  ndiversions = 0
1624  i0 = 1
1625  this%iadiv(1) = i0
1626  do n = 1, this%maxbound
1627  ndiversions = ndiversions + this%ndiv(n)
1628  i0 = i0 + this%ndiv(n)
1629  this%iadiv(n + 1) = i0
1630  end do
1631  !
1632  ! -- reallocate memory for diversions
1633  if (ndiversions > 0) then
1634  call mem_reallocate(this%divreach, ndiversions, 'DIVREACH', &
1635  this%memoryPath)
1636  allocate (this%divcprior(ndiversions))
1637  call mem_reallocate(this%divflow, ndiversions, 'DIVFLOW', this%memoryPath)
1638  call mem_reallocate(this%divq, ndiversions, 'DIVQ', this%memoryPath)
1639  end if
1640  !
1641  ! -- initialize diversion flow
1642  do n = 1, ndiversions
1643  this%divflow(n) = dzero
1644  this%divq(n) = dzero
1645  end do
1646  !
1647  ! -- read diversions
1648  call this%parser%GetBlock('DIVERSIONS', isfound, ierr, &
1649  supportopenclose=.true., &
1650  blockrequired=.false.)
1651  !
1652  ! -- parse reach connectivity block if detected
1653  if (isfound) then
1654  if (this%idiversions /= 0) then
1655  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
1656  ' DIVERSIONS'
1657  !
1658  ! -- allocate and initialize local variables for diversions
1659  ndiv = 0
1660  do n = 1, this%maxbound
1661  ndiv = ndiv + this%ndiv(n)
1662  end do
1663  allocate (iachk(this%maxbound + 1))
1664  allocate (nboundchk(ndiv))
1665  iachk(1) = 1
1666  do n = 1, this%maxbound
1667  iachk(n + 1) = iachk(n) + this%ndiv(n)
1668  end do
1669  do n = 1, ndiv
1670  nboundchk(n) = 0
1671  end do
1672  !
1673  ! -- read diversion data
1674  do
1675  call this%parser%GetNextLine(endofblock)
1676  if (endofblock) exit
1677  !
1678  ! -- get reach number
1679  n = this%parser%GetInteger()
1680  if (n < 1 .or. n > this%maxbound) then
1681  write (cnum, '(i0)') n
1682  errmsg = 'Reach number should be between 1 and '// &
1683  trim(cnum)//'.'
1684  call store_error(errmsg)
1685  cycle
1686  end if
1687  !
1688  ! -- make sure reach has at least one diversion
1689  if (this%ndiv(n) < 1) then
1690  write (cnum, '(i0)') n
1691  errmsg = 'Diversions cannot be specified '// &
1692  'for reach '//trim(cnum)
1693  call store_error(errmsg)
1694  cycle
1695  end if
1696  !
1697  ! -- read diversion number
1698  ival = this%parser%GetInteger()
1699  if (ival < 1 .or. ival > this%ndiv(n)) then
1700  write (cnum, '(i0)') n
1701  errmsg = 'Reach '//trim(cnum)
1702  write (cnum, '(i0)') this%ndiv(n)
1703  errmsg = trim(errmsg)//' diversion number should be between '// &
1704  '1 and '//trim(cnum)//'.'
1705  call store_error(errmsg)
1706  cycle
1707  end if
1708 
1709  ! -- increment nboundchk
1710  ipos = iachk(n) + ival - 1
1711  nboundchk(ipos) = nboundchk(ipos) + 1
1712 
1713  idiv = ival
1714  !
1715  ! -- get target reach for diversion
1716  ival = this%parser%GetInteger()
1717  if (ival < 1 .or. ival > this%maxbound) then
1718  write (cnum, '(i0)') ival
1719  errmsg = 'Diversion target reach number should be '// &
1720  'between 1 and '//trim(cnum)//'.'
1721  call store_error(errmsg)
1722  cycle
1723  end if
1724  idivreach = ival
1725  jpos = this%iadiv(n) + idiv - 1
1726  this%divreach(jpos) = idivreach
1727  !
1728  ! -- get cprior
1729  call this%parser%GetStringCaps(cval)
1730  ival = -1
1731  select case (cval)
1732  case ('UPTO')
1733  ival = 0
1734  case ('THRESHOLD')
1735  ival = -1
1736  case ('FRACTION')
1737  ival = -2
1738  case ('EXCESS')
1739  ival = -3
1740  case default
1741  errmsg = 'Invalid cprior type '//trim(cval)//'.'
1742  call store_error(errmsg)
1743  end select
1744  !
1745  ! -- set cprior for diversion
1746  this%divcprior(jpos) = cval
1747  end do
1748 
1749  write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))// &
1750  ' DIVERSIONS'
1751 
1752  do n = 1, this%maxbound
1753  do j = 1, this%ndiv(n)
1754  ipos = iachk(n) + j - 1
1755  !
1756  ! -- check for missing or duplicate reach diversions
1757  if (nboundchk(ipos) == 0) then
1758  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
1759  'No data specified for reach', n, 'diversion', j
1760  call store_error(errmsg)
1761  else if (nboundchk(ipos) > 1) then
1762  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
1763  'Data for reach', n, 'diversion', j, &
1764  'specified', nboundchk(ipos), 'times'
1765  call store_error(errmsg)
1766  end if
1767  end do
1768  end do
1769  !
1770  ! -- deallocate local variables
1771  deallocate (iachk)
1772  deallocate (nboundchk)
1773  else
1774  !
1775  ! -- error condition
1776  write (errmsg, '(a,1x,a)') &
1777  'A diversions block should not be', &
1778  'specified if diversions are not specified.'
1779  call store_error(errmsg)
1780  end if
1781  else
1782  if (this%idiversions /= 0) then
1783  call store_error('REQUIRED DIVERSIONS BLOCK NOT FOUND.')
1784  end if
1785  end if
1786  !
1787  ! -- write summary of diversion error messages
1788  if (count_errors() > 0) then
1789  call this%parser%StoreErrorUnit()
1790  end if
Here is the call graph for this function:

◆ sfr_read_initial_stages()

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

Method to read initialstages data for each reach for the SFR package.

Parameters
[in,out]thisSfrType object

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

1798  ! -- modules
1800  ! -- dummy
1801  class(SfrType), intent(inout) :: this !< SfrType object
1802  ! -- local
1803  integer(I4B) :: n
1804  integer(I4B) :: ierr
1805  logical(LGP) :: isfound
1806  logical(LGP) :: endOfBlock
1807  integer(I4B) :: i
1808  real(DP) :: rval
1809  integer, allocatable, dimension(:) :: nboundchk
1810  !
1811  ! -- read data
1812  call this%parser%GetBlock('INITIALSTAGES', isfound, ierr, &
1813  supportopenclose=.true., &
1814  blockrequired=.false.)
1815  !
1816  ! -- parse block if detected
1817  if (isfound) then
1818  write (this%iout, '(/1x,a)') &
1819  'PROCESSING '//trim(adjustl(this%text))//' INITIALSTAGES'
1820 
1821  allocate (nboundchk(this%maxbound))
1822  do n = 1, this%maxbound
1823  nboundchk(n) = 0
1824  end do
1825 
1826  do
1827  call this%parser%GetNextLine(endofblock)
1828  if (endofblock) exit
1829 
1830  ! -- read reach number
1831  n = this%parser%GetInteger()
1832 
1833  if (n < 1 .or. n > this%maxbound) then
1834  write (errmsg, '(a,i0,a,1x,i0,a)') &
1835  'Reach number (', n, ') must be greater than 0 and less &
1836  &than or equal to', this%maxbound, '.'
1837  call store_error(errmsg)
1838  cycle
1839  end if
1840 
1841  ! -- increment nboundchk
1842  nboundchk(n) = nboundchk(n) + 1
1843 
1844  rval = this%parser%GetDouble()
1845  this%stage(n) = rval
1846  this%depth(n) = rval - this%strtop(n)
1847 
1848  if (rval < this%strtop(n)) then
1849  write (errmsg, '(a,g0,a,1x,i0,1x,a,g0,a)') &
1850  'Initial stage (', rval, ') for reach', n, &
1851  'is less than the reach top (', this%strtop(n), ').'
1852  call store_error(errmsg)
1853  end if
1854  end do
1855 
1856  write (this%iout, '(1x,a)') &
1857  'END OF '//trim(adjustl(this%text))//' INITIALSTAGES'
1858 
1859  !
1860  ! -- Check to make sure that every reach is specified and that no reach
1861  ! is specified more than once.
1862  do i = 1, this%maxbound
1863  if (nboundchk(i) == 0) then
1864  write (errmsg, '(a,i0,1x,a)') &
1865  'Information for reach ', i, 'not specified in initialstages block.'
1866  call store_error(errmsg)
1867  else if (nboundchk(i) > 1) then
1868  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
1869  'Initial stage information specified', &
1870  nboundchk(i), 'times for reach', i
1871  call store_error(errmsg)
1872  end if
1873  end do
1874  deallocate (nboundchk)
1875  else
1876  ! -- set default initial stage based on a zero depth
1877  if (this%istorage == 1) then
1878  do n = 1, this%maxbound
1879  rval = this%strtop(n)
1880  this%stage(n) = rval
1881  end do
1882  end if
1883  end if
1884  !
1885  ! -- terminate if errors encountered in reach block
1886  if (count_errors() > 0) then
1887  call this%parser%StoreErrorUnit()
1888  end if
subroutine, public read_value_or_time_series_adv(textInput, ii, jj, bndElem, pkgName, auxOrBnd, tsManager, iprpak, varName)
Call this subroutine from advanced packages to define timeseries link for a variable (varName).
Here is the call graph for this function:

◆ sfr_read_packagedata()

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

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

Parameters
[in,out]thisSfrType object

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

934  ! -- modules
936  ! -- dummy
937  class(SfrType), intent(inout) :: this !< SfrType object
938  ! -- local
939  character(len=LINELENGTH) :: text
940  character(len=LINELENGTH) :: cellid
941  character(len=10) :: cnum
942  character(len=LENBOUNDNAME) :: bndName
943  character(len=LENBOUNDNAME) :: bndNameTemp
944  character(len=LENBOUNDNAME) :: hkname
945  character(len=LENBOUNDNAME) :: manningname
946  character(len=LENBOUNDNAME) :: ustrfname
947  character(len=50), dimension(:), allocatable :: caux
948  integer(I4B) :: n, ierr, ival
949  logical(LGP) :: isfound
950  logical(LGP) :: endOfBlock
951  integer(I4B) :: i
952  integer(I4B) :: ii
953  integer(I4B) :: jj
954  integer(I4B) :: iaux
955  integer(I4B) :: nconzero
956  integer(I4B) :: ipos
957  integer, allocatable, dimension(:) :: nboundchk
958  real(DP), pointer :: bndElem => null()
959  !
960  ! -- allocate space for checking sfr reach data
961  allocate (nboundchk(this%maxbound))
962  do i = 1, this%maxbound
963  nboundchk(i) = 0
964  end do
965  nconzero = 0
966  !
967  ! -- allocate local storage for aux variables
968  if (this%naux > 0) then
969  allocate (caux(this%naux))
970  end if
971  !
972  ! -- read reach data
973  call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
974  supportopenclose=.true.)
975  !
976  ! -- parse reaches block if detected
977  if (isfound) then
978  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
979  ' PACKAGEDATA'
980  do
981  call this%parser%GetNextLine(endofblock)
982  if (endofblock) exit
983  ! -- read reach number
984  n = this%parser%GetInteger()
985 
986  if (n < 1 .or. n > this%maxbound) then
987  write (errmsg, '(a,1x,a,1x,i0)') &
988  'Reach number (rno) must be greater than 0 and less', &
989  'than or equal to', this%maxbound
990  call store_error(errmsg)
991  cycle
992  end if
993 
994  ! -- increment nboundchk
995  nboundchk(n) = nboundchk(n) + 1
996  !
997  ! -- get model node number
998  call this%parser%GetCellid(this%dis%ndim, cellid, flag_string=.true.)
999  this%igwfnode(n) = this%dis%noder_from_cellid(cellid, this%inunit, &
1000  this%iout, &
1001  flag_string=.true., &
1002  allow_zero=.true.)
1003  this%igwftopnode(n) = this%igwfnode(n)
1004  !
1005  ! -- read the cellid string and determine if 'none' is specified
1006  if (this%igwfnode(n) < 1) then
1007  this%ianynone = this%ianynone + 1
1008  call upcase(cellid)
1009  if (cellid == 'NONE') then
1010  call this%parser%GetStringCaps(cellid)
1011  !
1012  ! -- create warning message
1013  write (cnum, '(i0)') n
1014  warnmsg = 'CELLID for unconnected reach '//trim(cnum)// &
1015  ' specified to be NONE. Unconnected reaches '// &
1016  'should be specified with a zero for each grid '// &
1017  'dimension. For example, for a DIS grid a CELLID '// &
1018  'of 0 0 0 should be specified for unconnected reaches'
1019  !
1020  ! -- create deprecation warning
1021  call deprecation_warning('PACKAGEDATA', 'CELLID=NONE', '6.4.3', &
1022  warnmsg, this%parser%GetUnit())
1023  else
1024 
1025  end if
1026  end if
1027  ! -- get reach length
1028  this%length(n) = this%parser%GetDouble()
1029  ! -- get reach width
1030  this%width(n) = this%parser%GetDouble()
1031  ! -- get reach slope
1032  this%slope(n) = this%parser%GetDouble()
1033  ! -- get reach stream bottom
1034  this%strtop(n) = this%parser%GetDouble()
1035  ! -- get reach bed thickness
1036  this%bthick(n) = this%parser%GetDouble()
1037  ! -- get reach bed hk
1038  call this%parser%GetStringCaps(hkname)
1039  ! -- get reach roughness
1040  call this%parser%GetStringCaps(manningname)
1041  ! -- get number of connections for reach
1042  ival = this%parser%GetInteger()
1043  this%nconnreach(n) = ival
1044  this%nconn = this%nconn + ival
1045  if (ival < 0) then
1046  write (errmsg, '(a,1x,i0,1x,a,i0,a)') &
1047  'NCON for reach', n, &
1048  'must be greater than or equal to 0 (', ival, ').'
1049  call store_error(errmsg)
1050  else if (ival == 0) then
1051  nconzero = nconzero + 1
1052  end if
1053  ! -- get upstream fraction for reach
1054  call this%parser%GetString(ustrfname)
1055  ! -- get number of diversions for reach
1056  ival = this%parser%GetInteger()
1057  this%ndiv(n) = ival
1058  if (ival > 0) then
1059  this%idiversions = 1
1060  else if (ival < 0) then
1061  ival = 0
1062  end if
1063 
1064  ! -- get aux data
1065  do iaux = 1, this%naux
1066  call this%parser%GetString(caux(iaux))
1067  end do
1068 
1069  ! -- set default bndName
1070  write (cnum, '(i10.10)') n
1071  bndname = 'Reach'//cnum
1072 
1073  ! -- get reach name
1074  if (this%inamedbound /= 0) then
1075  call this%parser%GetStringCaps(bndnametemp)
1076  if (bndnametemp /= '') then
1077  bndname = bndnametemp
1078  end if
1079  !this%boundname(n) = bndName
1080  end if
1081  this%sfrname(n) = bndname
1082  !
1083  ! -- set reach hydraulic conductivity
1084  text = hkname
1085  jj = 1 !for 'BEDK'
1086  bndelem => this%hk(n)
1087  call read_value_or_time_series_adv(text, n, jj, bndelem, &
1088  this%packName, 'BND', &
1089  this%tsManager, this%iprpak, &
1090  'BEDK')
1091  !
1092  ! -- set Mannings
1093  text = manningname
1094  jj = 1 !for 'MANNING'
1095  bndelem => this%rough(n)
1096  call read_value_or_time_series_adv(text, n, jj, bndelem, &
1097  this%packName, 'BND', &
1098  this%tsManager, this%iprpak, &
1099  'MANNING')
1100  !
1101  ! -- set upstream fraction
1102  text = ustrfname
1103  jj = 1 ! For 'USTRF'
1104  bndelem => this%ustrf(n)
1105  call read_value_or_time_series_adv(text, n, jj, bndelem, &
1106  this%packName, 'BND', &
1107  this%tsManager, this%iprpak, 'USTRF')
1108  !
1109  ! -- get aux data
1110  do jj = 1, this%naux
1111  text = caux(jj)
1112  ii = n
1113  bndelem => this%rauxvar(jj, ii)
1114  call read_value_or_time_series_adv(text, ii, jj, bndelem, &
1115  this%packName, 'AUX', &
1116  this%tsManager, this%iprpak, &
1117  this%auxname(jj))
1118  end do
1119  !
1120  ! -- initialize sstage to the top of the reach
1121  ! this value would be used by simple routing reaches
1122  ! on kper = 1 and kstp = 1 if a stage is not specified
1123  ! on the status line for the reach
1124  this%sstage(n) = this%strtop(n)
1125 
1126  end do
1127  write (this%iout, '(1x,a)') &
1128  'END OF '//trim(adjustl(this%text))//' PACKAGEDATA'
1129  else
1130  call store_error('REQUIRED PACKAGEDATA BLOCK NOT FOUND.')
1131  end if
1132  !
1133  ! -- Check to make sure that every reach is specified and that no reach
1134  ! is specified more than once.
1135  do i = 1, this%maxbound
1136  if (nboundchk(i) == 0) then
1137  write (errmsg, '(a,i0,1x,a)') &
1138  'Information for reach ', i, 'not specified in packagedata block.'
1139  call store_error(errmsg)
1140  else if (nboundchk(i) > 1) then
1141  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
1142  'Reach information specified', nboundchk(i), 'times for reach', i
1143  call store_error(errmsg)
1144  end if
1145  end do
1146  deallocate (nboundchk)
1147  !
1148  ! -- Submit warning message if any reach has zero connections
1149  if (nconzero > 0) then
1150  write (warnmsg, '(a,1x,a,1x,a,1x,i0,1x, a)') &
1151  'SFR Package', trim(this%packName), &
1152  'has', nconzero, 'reach(es) with zero connections.'
1153  call store_warning(warnmsg)
1154  end if
1155  !
1156  ! -- terminate if errors encountered in reach block
1157  if (count_errors() > 0) then
1158  call this%parser%StoreErrorUnit()
1159  end if
1160  !
1161  ! -- initialize the cross-section data
1162  ipos = 1
1163  this%iacross(1) = ipos
1164  do i = 1, this%maxbound
1165  this%ncrosspts(i) = 1
1166  this%station(ipos) = this%width(i)
1167  this%xsheight(ipos) = dzero
1168  this%xsrough(ipos) = done
1169  ipos = ipos + 1
1170  this%iacross(i + 1) = ipos
1171  end do
1172  !
1173  ! -- deallocate local storage for aux variables
1174  if (this%naux > 0) then
1175  deallocate (caux)
1176  end if
Here is the call graph for this function:

◆ sfr_rp()

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

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

Parameters
[in,out]thisSfrType object

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

1896  ! -- modules
1897  use tdismodule, only: kper, nper
1900  ! -- dummy
1901  class(SfrType), intent(inout) :: this !< SfrType object
1902  ! -- local
1903  character(len=LINELENGTH) :: title
1904  character(len=LINELENGTH) :: line
1905  character(len=LINELENGTH) :: crossfile
1906  integer(I4B) :: ierr
1907  integer(I4B) :: n
1908  integer(I4B) :: ichkustrm
1909  integer(I4B) :: ichkcross
1910  integer(I4B) :: ncrossptstot
1911  logical(LGP) :: isfound
1912  logical(LGP) :: endOfBlock
1913  type(SfrCrossSection), pointer :: cross_data => null()
1914  ! -- formats
1915  character(len=*), parameter :: fmtblkerr = &
1916  &"('Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
1917  character(len=*), parameter :: fmtlsp = &
1918  &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
1919  character(len=*), parameter :: fmtnbd = &
1920  "(1X,/1X,'The number of active ',A,'S (',I6, &
1921  &') is greater than maximum (',I6,')')"
1922  !
1923  ! -- initialize flags
1924  ichkustrm = 0
1925  ichkcross = 0
1926  if (kper == 1) then
1927  ichkustrm = 1
1928  end if
1929  !
1930  ! -- set nbound to maxbound
1931  this%nbound = this%maxbound
1932  !
1933  ! -- Set ionper to the stress period number for which a new block of data
1934  ! will be read.
1935  if (this%ionper < kper) then
1936  !
1937  ! -- get period block
1938  call this%parser%GetBlock('PERIOD', isfound, ierr, &
1939  supportopenclose=.true., &
1940  blockrequired=.false.)
1941  if (isfound) then
1942  !
1943  ! -- read ionper and check for increasing period numbers
1944  call this%read_check_ionper()
1945  else
1946  !
1947  ! -- PERIOD block not found
1948  if (ierr < 0) then
1949  ! -- End of file found; data applies for remainder of simulation.
1950  this%ionper = nper + 1
1951  else
1952  ! -- Found invalid block
1953  call this%parser%GetCurrentLine(line)
1954  write (errmsg, fmtblkerr) adjustl(trim(line))
1955  call store_error(errmsg)
1956  call this%parser%StoreErrorUnit()
1957  end if
1958  end if
1959  end if
1960  !
1961  ! -- Read data if ionper == kper
1962  if (this%ionper == kper) then
1963  !
1964  ! -- create and initialize cross-section data
1965  call cross_section_cr(cross_data, this%iout, this%iprpak, this%maxbound)
1966  call cross_data%initialize(this%ncrossptstot, this%ncrosspts, &
1967  this%iacross, &
1968  this%station, this%xsheight, &
1969  this%xsrough)
1970  !
1971  ! -- setup table for period data
1972  if (this%iprpak /= 0) then
1973  !
1974  ! -- reset the input table object
1975  title = trim(adjustl(this%text))//' PACKAGE ('// &
1976  trim(adjustl(this%packName))//') DATA FOR PERIOD'
1977  write (title, '(a,1x,i6)') trim(adjustl(title)), kper
1978  call table_cr(this%inputtab, this%packName, title)
1979  call this%inputtab%table_df(1, 4, this%iout, finalize=.false.)
1980  text = 'NUMBER'
1981  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
1982  text = 'KEYWORD'
1983  call this%inputtab%initialize_column(text, 20, alignment=tableft)
1984  do n = 1, 2
1985  write (text, '(a,1x,i6)') 'VALUE', n
1986  call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
1987  end do
1988  end if
1989  !
1990  ! -- read data
1991  do
1992  call this%parser%GetNextLine(endofblock)
1993  if (endofblock) exit
1994  n = this%parser%GetInteger()
1995  if (n < 1 .or. n > this%maxbound) then
1996  write (errmsg, '(a,1x,a,1x,i0,a)') &
1997  'Reach number (RNO) must be greater than 0 and', &
1998  'less than or equal to', this%maxbound, '.'
1999  call store_error(errmsg)
2000  cycle
2001  end if
2002  !
2003  ! -- read data from the rest of the line
2004  call this%sfr_set_stressperiod(n, ichkustrm, crossfile)
2005  !
2006  ! -- write line to table
2007  if (this%iprpak /= 0) then
2008  call this%parser%GetCurrentLine(line)
2009  call this%inputtab%line_to_columns(line)
2010  end if
2011  !
2012  ! -- process cross-section file
2013  if (trim(adjustl(crossfile)) /= 'NONE') then
2014  call cross_data%read_table(n, this%width(n), &
2015  trim(adjustl(crossfile)))
2016  end if
2017  end do
2018  !
2019  ! -- write raw period data
2020  if (this%iprpak /= 0) then
2021  call this%inputtab%finalize_table()
2022  end if
2023  !
2024  ! -- finalize cross-sections
2025 
2026  !
2027  ! -- determine the current size of cross-section data
2028  ncrossptstot = cross_data%get_ncrossptstot()
2029  !
2030  ! -- reallocate sfr package cross-section data
2031  if (ncrossptstot /= this%ncrossptstot) then
2032  this%ncrossptstot = ncrossptstot
2033  call mem_reallocate(this%station, this%ncrossptstot, 'STATION', &
2034  this%memoryPath)
2035  call mem_reallocate(this%xsheight, this%ncrossptstot, 'XSHEIGHT', &
2036  this%memoryPath)
2037  call mem_reallocate(this%xsrough, this%ncrossptstot, 'XSROUGH', &
2038  this%memoryPath)
2039  end if
2040  !
2041  ! -- write cross-section data to the model listing file
2042  call cross_data%output(this%width, this%rough, kstp=1, kper=kper)
2043  !
2044  ! -- pack cross-section data
2045  call cross_data%pack(this%ncrossptstot, this%ncrosspts, &
2046  this%iacross, &
2047  this%station, &
2048  this%xsheight, &
2049  this%xsrough)
2050  !
2051  ! -- deallocate temporary local storage for reach cross-sections
2052  call cross_data%destroy()
2053  deallocate (cross_data)
2054  nullify (cross_data)
2055  !
2056  ! -- Reuse data from last stress period
2057  else
2058  write (this%iout, fmtlsp) trim(this%filtyp)
2059  end if
2060  !
2061  ! -- check upstream fraction values
2062  if (ichkustrm /= 0) then
2063  call this%sfr_check_ustrf()
2064  end if
2065  !
2066  ! -- write summary of package block error messages
2067  if (count_errors() > 0) then
2068  call this%parser%StoreErrorUnit()
2069  end if
integer(i4b), pointer, public nper
number of stress period
Definition: tdis.f90:21
Here is the call graph for this function:

◆ sfr_rp_obs()

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

Method to read and prepare observations for a SFR package.

Parameters
[in,out]thisSfrType object

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

3155  ! -- modules
3156  use tdismodule, only: kper
3157  ! -- dummy
3158  class(SfrType), intent(inout) :: this !< SfrType object
3159  ! -- local
3160  integer(I4B) :: i
3161  integer(I4B) :: j
3162  integer(I4B) :: nn1
3163  character(len=LENBOUNDNAME) :: bname
3164  logical(LGP) :: jfound
3165  class(ObserveType), pointer :: obsrv => null()
3166  ! -- formats
3167 10 format('Boundary "', a, '" for observation "', a, &
3168  '" is invalid in package "', a, '"')
3169 30 format('Boundary name not provided for observation "', a, &
3170  '" in package "', a, '"')
3171  !
3172  ! -- process each package observation
3173  ! only done the first stress period since boundaries are fixed
3174  ! for the simulation
3175  if (kper == 1) then
3176  do i = 1, this%obs%npakobs
3177  obsrv => this%obs%pakobs(i)%obsrv
3178  !
3179  ! -- get node number 1
3180  nn1 = obsrv%NodeNumber
3181  if (nn1 == namedboundflag) then
3182  bname = obsrv%FeatureName
3183  if (bname /= '') then
3184  ! -- Observation location(s) is(are) based on a boundary name.
3185  ! Iterate through all boundaries to identify and store
3186  ! corresponding index(indices) in bound array.
3187  jfound = .false.
3188  do j = 1, this%maxbound
3189  if (this%boundname(j) == bname) then
3190  jfound = .true.
3191  call obsrv%AddObsIndex(j)
3192  end if
3193  end do
3194  if (.not. jfound) then
3195  write (errmsg, 10) &
3196  trim(bname), trim(obsrv%name), trim(this%packName)
3197  call store_error(errmsg)
3198  end if
3199  else
3200  write (errmsg, 30) trim(obsrv%name), trim(this%packName)
3201  call store_error(errmsg)
3202  end if
3203  else if (nn1 < 1 .or. nn1 > this%maxbound) then
3204  write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
3205  trim(adjustl(obsrv%ObsTypeId)), &
3206  'reach must be greater than 0 and less than or equal to', &
3207  this%maxbound, '(specified value is ', nn1, ')'
3208  call store_error(errmsg)
3209  else
3210  if (obsrv%indxbnds_count == 0) then
3211  call obsrv%AddObsIndex(nn1)
3212  else
3213  errmsg = 'Programming error in sfr_rp_obs'
3214  call store_error(errmsg)
3215  end if
3216  end if
3217  !
3218  ! -- catch non-cumulative observation assigned to observation defined
3219  ! by a boundname that is assigned to more than one element
3220  if (obsrv%ObsTypeId == 'STAGE' .or. &
3221  obsrv%ObsTypeId == 'DEPTH' .or. &
3222  obsrv%ObsTypeId == 'WET-PERIMETER' .or. &
3223  obsrv%ObsTypeId == 'WET-AREA' .or. &
3224  obsrv%ObsTypeId == 'WET-WIDTH') then
3225  nn1 = obsrv%NodeNumber
3226  if (nn1 == namedboundflag) then
3227  if (obsrv%indxbnds_count > 1) then
3228  write (errmsg, '(a,3(1x,a))') &
3229  trim(adjustl(obsrv%ObsTypeId)), &
3230  'for observation', trim(adjustl(obsrv%Name)), &
3231  ' must be assigned to a reach with a unique boundname.'
3232  call store_error(errmsg)
3233  end if
3234  end if
3235  end if
3236  !
3237  ! -- check that node number 1 is valid; call store_error if not
3238  do j = 1, obsrv%indxbnds_count
3239  nn1 = obsrv%indxbnds(j)
3240  if (nn1 < 1 .or. nn1 > this%maxbound) then
3241  write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
3242  trim(adjustl(obsrv%ObsTypeId)), &
3243  'reach must be greater than 0 and less than or equal to', &
3244  this%maxbound, '(specified value is ', nn1, ')'
3245  call store_error(errmsg)
3246  end if
3247  end do
3248  end do
3249  !
3250  ! -- evaluate if there are any observation errors
3251  if (count_errors() > 0) then
3252  call this%parser%StoreErrorUnit()
3253  end if
3254  end if
Here is the call graph for this function:

◆ sfr_set_stressperiod()

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

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

Parameters
[in,out]thisSfrType object
[in]nreach number
[in,out]ichkustrmflag indicating if upstream fraction data specified
[in,out]crossfilecross-section file name

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

3305  ! -- modules
3307  ! -- dummy
3308  class(SfrType), intent(inout) :: this !< SfrType object
3309  integer(I4B), intent(in) :: n !< reach number
3310  integer(I4B), intent(inout) :: ichkustrm !< flag indicating if upstream fraction data specified
3311  character(len=LINELENGTH), intent(inout) :: crossfile !< cross-section file name
3312  ! -- local
3313  character(len=10) :: cnum
3314  character(len=LINELENGTH) :: text
3315  character(len=LINELENGTH) :: caux
3316  character(len=LINELENGTH) :: keyword
3317  integer(I4B) :: ival
3318  integer(I4B) :: ii
3319  integer(I4B) :: jj
3320  integer(I4B) :: idiv
3321  integer(I4B) :: ixserror
3322  character(len=10) :: cp
3323  real(DP) :: divq
3324  real(DP), pointer :: bndElem => null()
3325  !
3326  ! -- initialize variables
3327  crossfile = 'NONE'
3328  !
3329  ! -- read line
3330  call this%parser%GetStringCaps(keyword)
3331  select case (keyword)
3332  case ('STATUS')
3333  ichkustrm = 1
3334  call this%parser%GetStringCaps(text)
3335  if (text == 'INACTIVE') then
3336  this%iboundpak(n) = 0
3337  else if (text == 'ACTIVE') then
3338  this%iboundpak(n) = 1
3339  else if (text == 'SIMPLE') then
3340  this%iboundpak(n) = -1
3341  else
3342  write (errmsg, '(2a)') &
3343  'Unknown '//trim(this%text)//' sfr status keyword: ', trim(text)
3344  call store_error(errmsg)
3345  end if
3346  case ('BEDK')
3347  call this%parser%GetString(text)
3348  jj = 1 ! For 'BEDK'
3349  bndelem => this%hk(n)
3350  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3351  this%packName, 'BND', &
3352  this%tsManager, this%iprpak, &
3353  'BEDK')
3354  case ('MANNING')
3355  call this%parser%GetString(text)
3356  jj = 1 ! For 'MANNING'
3357  bndelem => this%rough(n)
3358  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3359  this%packName, 'BND', &
3360  this%tsManager, this%iprpak, &
3361  'MANNING')
3362  case ('STAGE')
3363  call this%parser%GetString(text)
3364  jj = 1 ! For 'STAGE'
3365  bndelem => this%sstage(n)
3366  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3367  this%packName, 'BND', &
3368  this%tsManager, this%iprpak, 'STAGE')
3369  case ('RAINFALL')
3370  call this%parser%GetString(text)
3371  jj = 1 ! For 'RAIN'
3372  bndelem => this%rain(n)
3373  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3374  this%packName, 'BND', &
3375  this%tsManager, this%iprpak, 'RAIN')
3376  case ('EVAPORATION')
3377  call this%parser%GetString(text)
3378  jj = 1 ! For 'EVAP'
3379  bndelem => this%evap(n)
3380  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3381  this%packName, 'BND', &
3382  this%tsManager, this%iprpak, &
3383  'EVAP')
3384  case ('RUNOFF')
3385  call this%parser%GetString(text)
3386  jj = 1 ! For 'RUNOFF'
3387  bndelem => this%runoff(n)
3388  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3389  this%packName, 'BND', &
3390  this%tsManager, this%iprpak, &
3391  'RUNOFF')
3392  case ('INFLOW')
3393  call this%parser%GetString(text)
3394  jj = 1 ! For 'INFLOW'
3395  bndelem => this%inflow(n)
3396  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3397  this%packName, 'BND', &
3398  this%tsManager, this%iprpak, &
3399  'INFLOW')
3400  case ('DIVERSION')
3401  !
3402  ! -- make sure reach has at least one diversion
3403  if (this%ndiv(n) < 1) then
3404  write (cnum, '(i0)') n
3405  errmsg = 'diversions cannot be specified for reach '//trim(cnum)
3406  call store_error(errmsg)
3407  end if
3408  !
3409  ! -- read diversion number
3410  ival = this%parser%GetInteger()
3411  if (ival < 1 .or. ival > this%ndiv(n)) then
3412  write (cnum, '(i0)') n
3413  errmsg = 'Reach '//trim(cnum)
3414  write (cnum, '(i0)') this%ndiv(n)
3415  errmsg = trim(errmsg)//' diversion number should be between 1 '// &
3416  'and '//trim(cnum)//'.'
3417  call store_error(errmsg)
3418  end if
3419  idiv = ival
3420  !
3421  ! -- read value
3422  call this%parser%GetString(text)
3423  ii = this%iadiv(n) + idiv - 1
3424  jj = 1 ! For 'DIVERSION'
3425  bndelem => this%divflow(ii)
3426  call read_value_or_time_series_adv(text, ii, jj, bndelem, &
3427  this%packName, 'BND', &
3428  this%tsManager, this%iprpak, &
3429  'DIVFLOW')
3430  !
3431  ! -- if diversion cprior is 'fraction', ensure that 0.0 <= fraction <= 1.0
3432  cp = this%divcprior(ii)
3433  divq = this%divflow(ii)
3434  if (cp == 'FRACTION' .and. (divq < dzero .or. divq > done)) then
3435  write (errmsg, '(a,1x,i0,a)') &
3436  'cprior is type FRACTION for diversion no.', ii, &
3437  ', but divflow not within the range 0.0 to 1.0'
3438  call store_error(errmsg)
3439  end if
3440  case ('UPSTREAM_FRACTION')
3441  ichkustrm = 1
3442  call this%parser%GetString(text)
3443  jj = 1 ! For 'USTRF'
3444  bndelem => this%ustrf(n)
3445  call read_value_or_time_series_adv(text, n, jj, bndelem, &
3446  this%packName, 'BND', &
3447  this%tsManager, this%iprpak, 'USTRF')
3448 
3449  case ('CROSS_SECTION')
3450  ixserror = 0
3451  !
3452  ! -- read FILE keyword
3453  call this%parser%GetStringCaps(keyword)
3454  select case (keyword)
3455  case ('TAB6')
3456  call this%parser%GetStringCaps(keyword)
3457  if (trim(adjustl(keyword)) /= 'FILEIN') then
3458  errmsg = 'TAB6 keyword must be followed by "FILEIN" '// &
3459  'then by filename.'
3460  call store_error(errmsg)
3461  ixserror = 1
3462  end if
3463  if (ixserror == 0) then
3464  call this%parser%GetString(crossfile)
3465  end if
3466  case default
3467  write (errmsg, '(a,1x,i4,1x,a)') &
3468  'CROSS-SECTION TABLE ENTRY for REACH ', n, &
3469  'MUST INCLUDE TAB6 KEYWORD'
3470  call store_error(errmsg)
3471  end select
3472 
3473  case ('AUXILIARY')
3474  call this%parser%GetStringCaps(caux)
3475  do jj = 1, this%naux
3476  if (trim(adjustl(caux)) /= trim(adjustl(this%auxname(jj)))) cycle
3477  call this%parser%GetString(text)
3478  ii = n
3479  bndelem => this%rauxvar(jj, ii)
3480  call read_value_or_time_series_adv(text, ii, jj, bndelem, &
3481  this%packName, 'AUX', &
3482  this%tsManager, this%iprpak, &
3483  this%auxname(jj))
3484  exit
3485  end do
3486 
3487  case default
3488  write (errmsg, '(a,a)') &
3489  'Unknown '//trim(this%text)//' sfr data keyword: ', &
3490  trim(keyword)//'.'
3491  call store_error(errmsg)
3492  end select
Here is the call graph for this function:

◆ sfr_setup_budobj()

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

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

Parameters
thisSfrType object

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

5031  ! -- dummy
5032  class(SfrType) :: this !< SfrType object
5033  ! -- local
5034  integer(I4B) :: nbudterm
5035  integer(I4B) :: i
5036  integer(I4B) :: n
5037  integer(I4B) :: n1
5038  integer(I4B) :: n2
5039  integer(I4B) :: maxlist
5040  integer(I4B) :: naux
5041  integer(I4B) :: idx
5042  real(DP) :: q
5043  character(len=LENBUDTXT) :: text
5044  character(len=LENBUDTXT), dimension(1) :: auxtxt
5045  !
5046  ! -- Determine the number of sfr budget terms. These are fixed for
5047  ! the simulation and cannot change. This includes FLOW-JA-FACE
5048  ! so they can be written to the binary budget files, but these internal
5049  ! flows are not included as part of the budget table.
5050  nbudterm = 8
5051  if (this%imover == 1) nbudterm = nbudterm + 2
5052  if (this%naux > 0) nbudterm = nbudterm + 1
5053  !
5054  ! -- set up budobj
5055  call budgetobject_cr(this%budobj, this%packName)
5056  call this%budobj%budgetobject_df(this%maxbound, nbudterm, 0, 0, &
5057  ibudcsv=this%ibudcsv)
5058  idx = 0
5059  !
5060  ! -- Go through and set up each budget term
5061  text = ' FLOW-JA-FACE'
5062  idx = idx + 1
5063  maxlist = this%nconn
5064  naux = 1
5065  auxtxt(1) = ' FLOW-AREA'
5066  call this%budobj%budterm(idx)%initialize(text, &
5067  this%name_model, &
5068  this%packName, &
5069  this%name_model, &
5070  this%packName, &
5071  maxlist, .false., .false., &
5072  naux, auxtxt)
5073  !
5074  ! -- store connectivity
5075  call this%budobj%budterm(idx)%reset(this%nconn)
5076  q = dzero
5077  do n = 1, this%maxbound
5078  n1 = n
5079  do i = this%ia(n) + 1, this%ia(n + 1) - 1
5080  n2 = this%ja(i)
5081  call this%budobj%budterm(idx)%update_term(n1, n2, q)
5082  end do
5083  end do
5084  !
5085  ! --
5086  text = ' GWF'
5087  idx = idx + 1
5088  maxlist = this%maxbound - this%ianynone
5089  naux = 1
5090  auxtxt(1) = ' FLOW-AREA'
5091  call this%budobj%budterm(idx)%initialize(text, &
5092  this%name_model, &
5093  this%packName, &
5094  this%name_model, &
5095  this%name_model, &
5096  maxlist, .false., .true., &
5097  naux, auxtxt)
5098  call this%budobj%budterm(idx)%reset(maxlist)
5099  q = dzero
5100  do n = 1, this%maxbound
5101  n2 = this%igwfnode(n)
5102  if (n2 > 0) then
5103  call this%budobj%budterm(idx)%update_term(n, n2, q)
5104  end if
5105  end do
5106  !
5107  ! --
5108  text = ' RAINFALL'
5109  idx = idx + 1
5110  maxlist = this%maxbound
5111  naux = 0
5112  call this%budobj%budterm(idx)%initialize(text, &
5113  this%name_model, &
5114  this%packName, &
5115  this%name_model, &
5116  this%packName, &
5117  maxlist, .false., .false., &
5118  naux)
5119  !
5120  ! --
5121  text = ' EVAPORATION'
5122  idx = idx + 1
5123  maxlist = this%maxbound
5124  naux = 0
5125  call this%budobj%budterm(idx)%initialize(text, &
5126  this%name_model, &
5127  this%packName, &
5128  this%name_model, &
5129  this%packName, &
5130  maxlist, .false., .false., &
5131  naux)
5132  !
5133  ! --
5134  text = ' RUNOFF'
5135  idx = idx + 1
5136  maxlist = this%maxbound
5137  naux = 0
5138  call this%budobj%budterm(idx)%initialize(text, &
5139  this%name_model, &
5140  this%packName, &
5141  this%name_model, &
5142  this%packName, &
5143  maxlist, .false., .false., &
5144  naux)
5145  !
5146  ! --
5147  text = ' EXT-INFLOW'
5148  idx = idx + 1
5149  maxlist = this%maxbound
5150  naux = 0
5151  call this%budobj%budterm(idx)%initialize(text, &
5152  this%name_model, &
5153  this%packName, &
5154  this%name_model, &
5155  this%packName, &
5156  maxlist, .false., .false., &
5157  naux)
5158  !
5159  ! --
5160  text = ' EXT-OUTFLOW'
5161  idx = idx + 1
5162  maxlist = this%maxbound
5163  naux = 0
5164  call this%budobj%budterm(idx)%initialize(text, &
5165  this%name_model, &
5166  this%packName, &
5167  this%name_model, &
5168  this%packName, &
5169  maxlist, .false., .false., &
5170  naux)
5171  !
5172  ! --
5173  text = ' STORAGE'
5174  idx = idx + 1
5175  maxlist = this%maxbound
5176  naux = 1
5177  auxtxt(1) = ' VOLUME'
5178  call this%budobj%budterm(idx)%initialize(text, &
5179  this%name_model, &
5180  this%packName, &
5181  this%name_model, &
5182  this%packName, &
5183  maxlist, .false., .false., &
5184  naux, auxtxt)
5185  !
5186  ! --
5187  if (this%imover == 1) then
5188  !
5189  ! --
5190  text = ' FROM-MVR'
5191  idx = idx + 1
5192  maxlist = this%maxbound
5193  naux = 0
5194  call this%budobj%budterm(idx)%initialize(text, &
5195  this%name_model, &
5196  this%packName, &
5197  this%name_model, &
5198  this%packName, &
5199  maxlist, .false., .false., &
5200  naux)
5201  !
5202  ! --
5203  text = ' TO-MVR'
5204  idx = idx + 1
5205  maxlist = this%maxbound
5206  naux = 0
5207  call this%budobj%budterm(idx)%initialize(text, &
5208  this%name_model, &
5209  this%packName, &
5210  this%name_model, &
5211  this%packName, &
5212  maxlist, .false., .false., &
5213  naux)
5214  end if
5215  !
5216  ! --
5217  naux = this%naux
5218  if (naux > 0) then
5219  !
5220  ! --
5221  text = ' AUXILIARY'
5222  idx = idx + 1
5223  maxlist = this%maxbound
5224  call this%budobj%budterm(idx)%initialize(text, &
5225  this%name_model, &
5226  this%packName, &
5227  this%name_model, &
5228  this%packName, &
5229  maxlist, .false., .false., &
5230  naux, this%auxname)
5231  end if
5232  !
5233  ! -- if sfr flow for each reach are written to the listing file
5234  if (this%iprflow /= 0) then
5235  call this%budobj%flowtable_df(this%iout, cellids='GWF')
5236  end if
Here is the call graph for this function:

◆ sfr_setup_tableobj()

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

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

Parameters
thisSfrType object

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

5478  ! -- dummy
5479  class(SfrType) :: this !< SfrType object
5480  ! -- local
5481  integer(I4B) :: nterms
5482  character(len=LINELENGTH) :: title
5483  character(len=LINELENGTH) :: text
5484  !
5485  ! -- setup stage table
5486  if (this%iprhed > 0) then
5487  !
5488  ! -- Determine the number of sfr budget terms. These are fixed for
5489  ! the simulation and cannot change. This includes FLOW-JA-FACE
5490  ! so they can be written to the binary budget files, but these internal
5491  ! flows are not included as part of the budget table.
5492  nterms = 8
5493  if (this%inamedbound == 1) then
5494  nterms = nterms + 1
5495  end if
5496  !
5497  ! -- set up table title
5498  title = trim(adjustl(this%text))//' PACKAGE ('// &
5499  trim(adjustl(this%packName))//') STAGES FOR EACH CONTROL VOLUME'
5500  !
5501  ! -- set up stage tableobj
5502  call table_cr(this%stagetab, this%packName, title)
5503  call this%stagetab%table_df(this%maxbound, nterms, this%iout, &
5504  transient=.true.)
5505  !
5506  ! -- Go through and set up table budget term
5507  if (this%inamedbound == 1) then
5508  text = 'NAME'
5509  call this%stagetab%initialize_column(text, lenboundname, &
5510  alignment=tableft)
5511  end if
5512  !
5513  ! -- reach number
5514  text = 'NUMBER'
5515  call this%stagetab%initialize_column(text, 10, alignment=tabcenter)
5516  !
5517  ! -- cellids
5518  text = 'CELLID'
5519  call this%stagetab%initialize_column(text, 20, alignment=tableft)
5520  !
5521  ! -- reach stage
5522  text = 'STAGE'
5523  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5524  !
5525  ! -- reach depth
5526  text = 'DEPTH'
5527  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5528  !
5529  ! -- reach width
5530  text = 'WIDTH'
5531  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5532  !
5533  ! -- gwf head
5534  text = 'GWF HEAD'
5535  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5536  !
5537  ! -- streambed conductance
5538  text = 'STREAMBED CONDUCTANCE'
5539  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5540  !
5541  ! -- streambed gradient
5542  text = 'STREAMBED GRADIENT'
5543  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5544  end if
Here is the call graph for this function:

◆ sfr_solve()

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

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

Parameters
thisSfrType object
[in]nreach number
[in]hgroundwater head in cell connected to reach
[in,out]hcofcoefficient term added to the diagonal
[in,out]rhsright-hand side term
[in]updateboolean indicating if the reach depth and stage variables should be updated to current iterate

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

3500  ! -- dummy
3501  class(SfrType) :: this !< SfrType object
3502  integer(I4B), intent(in) :: n !< reach number
3503  real(DP), intent(in) :: h !< groundwater head in cell connected to reach
3504  real(DP), intent(inout) :: hcof !< coefficient term added to the diagonal
3505  real(DP), intent(inout) :: rhs !< right-hand side term
3506  logical(LGP), intent(in), optional :: update !< boolean indicating if the reach depth and stage variables should be updated to current iterate
3507  ! -- local
3508  logical(LGP) :: lupdate
3509  integer(I4B) :: i
3510  integer(I4B) :: ii
3511  integer(I4B) :: n2
3512  real(DP) :: hgwf
3513  real(DP) :: sa
3514  real(DP) :: sa_wet
3515  real(DP) :: qu
3516  real(DP) :: qi
3517  real(DP) :: qr
3518  real(DP) :: qe
3519  real(DP) :: qro
3520  real(DP) :: qsrc
3521  real(DP) :: qfrommvr
3522  real(DP) :: qgwf
3523  real(DP) :: tp
3524  real(DP) :: bt
3525  real(DP) :: hsfr
3526  real(DP) :: qd
3527  real(DP) :: d1
3528  real(DP) :: sumleak
3529  real(DP) :: sumrch
3530  real(DP) :: gwfhcof
3531  real(DP) :: gwfrhs
3532  !
3533  ! -- Process optional dummy variables
3534  if (present(update)) then
3535  lupdate = update
3536  else
3537  lupdate = .true.
3538  end if
3539 
3540  ! -- initialize variables
3541  hcof = dzero
3542  rhs = dzero
3543  !
3544  if (this%iboundpak(n) == 0) then
3545  this%depth(n) = dzero
3546  this%stage(n) = dhnoflo
3547  this%usflow(n) = dzero
3548  this%simevap(n) = dzero
3549  this%simrunoff(n) = dzero
3550  this%dsflow(n) = dzero
3551  this%gwflow(n) = dzero
3552  else
3553  hgwf = h
3554  d1 = dzero
3555  qsrc = dzero
3556  qgwf = dzero
3557  qd = this%dsflow(n)
3558 
3559  ! -- calculate initial depth assuming a wide cross-section and
3560  ! ignore groundwater leakage
3561  ! -- calculate upstream flow
3562  qu = dzero
3563  do i = this%ia(n) + 1, this%ia(n + 1) - 1
3564  if (this%idir(i) < 0) cycle
3565  n2 = this%ja(i)
3566  do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
3567  if (this%idir(ii) > 0) cycle
3568  if (this%ja(ii) /= n) cycle
3569  qu = qu + this%qconn(ii)
3570  end do
3571  end do
3572  this%usflow(n) = qu
3573 
3574  ! -- calculate remaining terms
3575  sa = this%calc_surface_area(n)
3576  sa_wet = this%calc_surface_area_wet(n, this%depth(n))
3577  qi = this%inflow(n)
3578  qr = this%rain(n) * sa
3579  qe = this%evap(n) * sa_wet
3580  qro = this%runoff(n)
3581 
3582  ! -- Water mover term; assume that it goes in at the upstream end of the reach
3583  qfrommvr = dzero
3584  if (this%imover == 1) then
3585  qfrommvr = this%pakmvrobj%get_qfrommvr(n)
3586  end if
3587 
3588  ! -- calculate downstream flow ignoring groundwater leakage
3589  qsrc = qu + qi + qr - qe + qro + qfrommvr
3590 
3591  ! -- adjust runoff or evaporation if sum of sources is negative
3592  call this%sfr_adjust_ro_ev(qsrc, qu, qi, qr, qro, qe, qfrommvr)
3593 
3594  ! -- set simulated evaporation and runoff
3595  this%simevap(n) = qe
3596  this%simrunoff(n) = qro
3597 
3598  ! -- calculate reach flow using appropriate method
3599  if (this%iboundpak(n) < 0) then
3600  call this%sfr_calc_constant(n, d1, hgwf, qgwf, qd)
3601  else
3602  if (this%gwfiss == 0 .and. this%istorage == 1) then
3603  call this%sfr_calc_transient(n, d1, hgwf, qu, qi, &
3604  qfrommvr, qr, qe, qro, &
3605  qgwf, qd)
3606  else
3607  call this%sfr_calc_steady(n, d1, hgwf, qu, qi, &
3608  qfrommvr, qr, qe, qro, &
3609  qgwf, qd)
3610  end if
3611  end if
3612 
3613  ! -- update sfr stage
3614  tp = this%strtop(n)
3615  bt = tp - this%bthick(n)
3616  hsfr = tp + d1
3617 
3618  ! -- update stored values
3619  if (lupdate) then
3620  ! -- save depth and calculate stage
3621  this%depth(n) = d1
3622  this%stage(n) = hsfr
3623  ! -- update flows
3624  call this%sfr_update_flows(n, qd, qgwf)
3625  end if
3626 
3627  ! -- calculate sumleak and sumrch
3628  sumleak = dzero
3629  sumrch = dzero
3630  if (this%gwfiss == 0) then
3631  sumleak = qgwf
3632  else
3633  sumleak = qgwf
3634  end if
3635  if (hgwf < bt) then
3636  sumrch = qgwf
3637  end if
3638 
3639  ! -- make final qgwf calculation and obtain
3640  ! gwfhcof and gwfrhs values
3641  call this%sfr_calc_qgwf(n, d1, hgwf, qgwf, gwfhcof, gwfrhs)
3642 
3643  ! -- update hcof and rhs terms
3644  if (abs(sumleak) > dzero) then
3645  ! -- stream leakage is not head dependent
3646  if (hgwf < bt) then
3647  rhs = rhs - sumrch
3648  ! -- stream leakage is head dependent
3649  else if ((sumleak - qsrc) < -dem30) then
3650  if (this%gwfiss == 0) then
3651  rhs = rhs + gwfrhs - sumrch
3652  else
3653  rhs = rhs + gwfrhs
3654  end if
3655  hcof = gwfhcof
3656  ! -- place holder for UZF
3657  else
3658  if (this%gwfiss == 0) then
3659  rhs = rhs - sumleak - sumrch
3660  else
3661  rhs = rhs - sumleak
3662  end if
3663  end if
3664 
3665  ! -- add groundwater leakage
3666  else if (hgwf < bt) then
3667  rhs = rhs - sumrch
3668  end if
3669  end if

◆ sfr_update_flows()

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

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

Parameters
[in,out]thisSfrType object
[in]nreach number
[in,out]qddownstream reach flow
[in]qgwfgroundwater leakage for reach

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

3678  ! -- dummy
3679  class(SfrType), intent(inout) :: this !< SfrType object
3680  integer(I4B), intent(in) :: n !< reach number
3681  real(DP), intent(inout) :: qd !< downstream reach flow
3682  real(DP), intent(in) :: qgwf !< groundwater leakage for reach
3683  ! -- local
3684  integer(I4B) :: i
3685  integer(I4B) :: n2
3686  integer(I4B) :: idiv
3687  integer(I4B) :: jpos
3688  real(DP) :: qdiv
3689  real(DP) :: f
3690  !
3691  ! -- update reach terms
3692  !
3693  ! -- save final downstream stream flow
3694  this%dsflow(n) = qd
3695  !
3696  ! -- save groundwater leakage
3697  this%gwflow(n) = qgwf
3698  !
3699  ! -- route downstream flow
3700  if (qd > dzero) then
3701  !
3702  ! -- route water to diversions
3703  do i = this%ia(n) + 1, this%ia(n + 1) - 1
3704  if (this%idir(i) > 0) cycle
3705  idiv = this%idiv(i)
3706  if (idiv == 0) cycle
3707  jpos = this%iadiv(n) + idiv - 1
3708  call this%sfr_calc_div(n, idiv, qd, qdiv)
3709  this%qconn(i) = qdiv
3710  this%divq(jpos) = qdiv
3711  end do
3712  !
3713  ! -- Mover terms: store outflow after diversion loss
3714  ! as qformvr and reduce outflow (qd)
3715  ! by how much was actually sent to the mover
3716  if (this%imover == 1) then
3717  call this%pakmvrobj%accumulate_qformvr(n, qd)
3718  qd = max(qd - this%pakmvrobj%get_qtomvr(n), dzero)
3719  end if
3720  !
3721  ! -- route remaining water to downstream reaches
3722  do i = this%ia(n) + 1, this%ia(n + 1) - 1
3723  if (this%idir(i) > 0) cycle
3724  if (this%idiv(i) > 0) cycle
3725  n2 = this%ja(i)
3726  if (this%iboundpak(n2) == 0) cycle
3727  f = this%ustrf(n2) / this%ftotnd(n)
3728  this%qconn(i) = qd * f
3729  end do
3730  else
3731  do i = this%ia(n) + 1, this%ia(n + 1) - 1
3732  if (this%idir(i) > 0) cycle
3733  this%qconn(i) = dzero
3734  idiv = this%idiv(i)
3735  if (idiv == 0) cycle
3736  jpos = this%iadiv(n) + idiv - 1
3737  this%divq(jpos) = dzero
3738  end do
3739  end if

Variable Documentation

◆ ftype

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

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

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

◆ text

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

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

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