gwfbuymodule Module Reference

Data Types
type	concentrationpointer
type	gwfbuytype

Functions/Subroutines
real(dp) function	calcdens (denseref, drhodc, crhoref, conc)
	Generic function to calculate fluid density from concentration.
subroutine, public	buy_cr (buyobj, name_model, input_mempath, inunit, iout)
	Create a new BUY object.
subroutine	buy_df (this, dis, buy_input)
	Read options and package data, or set from argument.
subroutine	buy_ar (this, npf, ibound)
	Allocate and Read.
subroutine	buy_ar_bnd (this, packobj, hnew)
	Buoyancy ar_bnd routine to activate density in packages.
subroutine	buy_rp (this)
	Check for new buy period data.
subroutine	buy_ad (this)
	Advance.
subroutine	buy_cf (this, kiter)
	Fill coefficients.
subroutine	buy_cf_bnd (this, packobj, hnew)
	Fill coefficients.
real(dp) function	get_bnd_density (n, locdense, locconc, denseref, drhodc, crhoref, ctemp, auxvar)
	Return the density of the boundary package using one of several different options in the following order of priority:
subroutine	buy_cf_ghb (packobj, dense, elev, denseref, locelev, locdense, locconc, drhodc, crhoref, ctemp)
	Fill ghb coefficients.
subroutine	calc_ghb_hcof_rhs_terms (denseref, denseghb, densenode, elevghb, elevnode, hghb, cond, rhsterm, hcofterm)
	Calculate density hcof and rhs terms for ghb conditions.
subroutine	buy_cf_riv (packobj, hnew, dense, elev, denseref, locelev, locdense, locconc, drhodc, crhoref, ctemp)
	Fill riv coefficients.
subroutine	buy_cf_drn (packobj, dense, denseref)
	Fill drn coefficients.
subroutine	buy_cf_lak (packobj, hnew, dense, elev, denseref, locdense, locconc, drhodc, crhoref, ctemp, iform)
	Pass density information into lak package; density terms are calculated in the lake package as part of lak_calculate_density_exchange method.
subroutine	buy_cf_sfr (packobj, hnew, dense, elev, denseref, locdense, locconc, drhodc, crhoref, ctemp, iform)
	Pass density information into sfr package; density terms are calculated in the sfr package as part of sfr_calculate_density_exchange method.
subroutine	buy_cf_maw (packobj, hnew, dense, elev, denseref, locdense, locconc, drhodc, crhoref, ctemp, iform)
	Pass density information into maw package; density terms are calculated in the maw package as part of maw_calculate_density_exchange method.
subroutine	buy_fc (this, kiter, matrix_sln, idxglo, rhs, hnew)
	Fill coefficients.
subroutine	buy_ot_dv (this, idvfl)
	Save density array to binary file.
subroutine	buy_cq (this, hnew, flowja)
	Add buy term to flowja.
subroutine	buy_da (this)
	Deallocate.
subroutine	source_dimensions (this)
	@ brief Source dimensions for package
subroutine	source_packagedata (this)
	@ brief source packagedata for package
subroutine	set_packagedata (this, input_data)
	Sets package data instead of reading from file.
subroutine	calcbuy (this, n, m, icon, hn, hm, buy)
	Calculate buyancy term for this connection.
subroutine	calchhterms (this, n, m, icon, hn, hm, rhsterm, amatnn, amatnm)
	Calculate hydraulic head term for this connection.
subroutine	buy_calcdens (this)
	calculate fluid density from concentration
subroutine	buy_calcelev (this)
	Calculate cell elevations to use in density flow equations.
subroutine	allocate_scalars (this)
	Allocate scalars used by the package.
subroutine	allocate_arrays (this, nodes)
	Allocate arrays used by the package.
subroutine	source_options (this)
	@ brief Source input options
subroutine	log_options (this, found, densityfile)
	@ brief log input options
subroutine	set_options (this, input_data)
	Sets options as opposed to reading them from a file.
subroutine	set_concentration_pointer (this, modelname, conc, icbund)
	Pass in a gwt model name, concentration array and ibound, and store a pointer to these in the BUY package so that density can be calculated from them.

Function/Subroutine Documentation

allocate_arrays()

subroutine gwfbuymodule::allocate_arrays ( class(gwfbuytype) this, integer(i4b), intent(in) nodes )

private

Definition at line 1326 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    integer(I4B), intent(in) :: nodes
    ! -- local
    integer(I4B) :: i
    !
    ! -- Allocate
    call mem_allocate(this%dense, nodes, 'DENSE', this%memoryPath)
    call mem_allocate(this%concbuy, 0, 'CONCBUY', this%memoryPath)
    call mem_allocate(this%elev, nodes, 'ELEV', this%memoryPath)
    call mem_allocate(this%drhodc, this%nrhospecies, 'DRHODC', this%memoryPath)
    call mem_allocate(this%crhoref, this%nrhospecies, 'CRHOREF', this%memoryPath)
    call mem_allocate(this%ctemp, this%nrhospecies, 'CTEMP', this%memoryPath)
    allocate (this%cmodelname(this%nrhospecies))
    allocate (this%cauxspeciesname(this%nrhospecies))
    allocate (this%modelconc(this%nrhospecies))
    !
    ! -- Initialize
    do i = 1, nodes
      this%dense(i) = this%denseref
      this%elev(i) = dzero
    end do
    !
    ! -- Initialize nrhospecies arrays
    do i = 1, this%nrhospecies
      this%drhodc(i) = dzero
      this%crhoref(i) = dzero
      this%ctemp(i) = dzero
      this%cmodelname(i) = ''
      this%cauxspeciesname(i) = ''
    end do

allocate_scalars()

subroutine gwfbuymodule::allocate_scalars ( class(gwfbuytype) this )

private

Definition at line 1290 of file gwf-buy.f90.

    ! -- modules
    use constantsmodule, only: dzero
    ! -- dummy
    class(GwfBuyType) :: this
    ! -- local
    !
    ! -- allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- Allocate
    call mem_allocate(this%ioutdense, 'IOUTDENSE', this%memoryPath)
    call mem_allocate(this%iform, 'IFORM', this%memoryPath)
    call mem_allocate(this%ireadelev, 'IREADELEV', this%memoryPath)
    call mem_allocate(this%ireadconcbuy, 'IREADCONCBUY', this%memoryPath)
    call mem_allocate(this%iconcset, 'ICONCSET', this%memoryPath)
    call mem_allocate(this%denseref, 'DENSEREF', this%memoryPath)
    !
    call mem_allocate(this%nrhospecies, 'NRHOSPECIES', this%memoryPath)
    !
    ! -- Initialize
    this%ioutdense = 0
    this%ireadelev = 0
    this%iconcset = 0
    this%ireadconcbuy = 0
    this%denseref = 1000.d0
    !
    this%nrhospecies = 0
    !
    ! -- Initialize default to LHS implementation of hydraulic head formulation
    this%iform = 2
    this%iasym = 1

buy_ad()

subroutine gwfbuymodule::buy_ad ( class(gwfbuytype) this )

private

Definition at line 275 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    !
    ! -- update density using the last concentration
    call this%buy_calcdens()

buy_ar()

subroutine gwfbuymodule::buy_ar ( class(gwfbuytype) this, type(gwfnpftype), intent(in), pointer npf, integer(i4b), dimension(:), pointer ibound )

private

Definition at line 174 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    type(GwfNpfType), pointer, intent(in) :: npf
    integer(I4B), dimension(:), pointer :: ibound
    !
    ! -- store pointers to arguments that were passed in
    this%npf => npf
    this%ibound => ibound
    !
    ! -- Ensure NPF XT3D is not on
    if (this%npf%ixt3d /= 0) then
      call store_error('Error in model '//trim(this%name_model)// &
                       '.  The XT3D option cannot be used with the BUY Package.')
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- Calculate cell elevations
    call this%buy_calcelev()

Here is the call graph for this function:

buy_ar_bnd()

subroutine gwfbuymodule::buy_ar_bnd ( class(gwfbuytype) this, class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) hnew )

private

This routine is called from gwf_ar() as it goes through each package

Definition at line 199 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    use lakmodule, only: laktype
    use sfrmodule, only: sfrtype
    use mawmodule, only: mawtype
    ! -- dummy
    class(GwfBuyType) :: this
    class(BndType), pointer :: packobj
    real(DP), intent(in), dimension(:) :: hnew
    !
    ! -- Add density terms based on boundary package type
    select case (packobj%filtyp)
    case ('LAK')
      !
      ! -- activate density for lake package
      select type (packobj)
      type is (laktype)
        call packobj%lak_activate_density()
      end select
      !
    case ('SFR')
      !
      ! -- activate density for sfr package
      select type (packobj)
      type is (sfrtype)
        call packobj%sfr_activate_density()
      end select
      !
    case ('MAW')
      !
      ! -- activate density for maw package
      select type (packobj)
      type is (mawtype)
        call packobj%maw_activate_density()
      end select
      !
    case default
      !
      ! -- nothing
    end select

buy_calcdens()

subroutine gwfbuymodule::buy_calcdens ( class(gwfbuytype) this )

private

Definition at line 1248 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
 
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: i
    !
    ! -- Calculate the density using the specified concentration array
    do n = 1, this%dis%nodes
      do i = 1, this%nrhospecies
        if (this%modelconc(i)%icbund(n) == 0) then
          this%ctemp(i) = dzero
        else
          this%ctemp(i) = this%modelconc(i)%conc(n)
        end if
      end do
      this%dense(n) = calcdens(this%denseref, this%drhodc, this%crhoref, &
                               this%ctemp)
    end do

Here is the call graph for this function:

buy_calcelev()

subroutine gwfbuymodule::buy_calcelev ( class(gwfbuytype) this )

private

Definition at line 1272 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    ! -- local
    integer(I4B) :: n
    real(DP) :: tp, bt, frac
    !
    ! -- Calculate the elev array
    do n = 1, this%dis%nodes
      tp = this%dis%top(n)
      bt = this%dis%bot(n)
      frac = this%npf%sat(n)
      this%elev(n) = bt + dhalf * frac * (tp - bt)
    end do

buy_cf()

subroutine gwfbuymodule::buy_cf ( class(gwfbuytype) this, integer(i4b) kiter )

private

Definition at line 285 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    integer(I4B) :: kiter
    !
    ! -- Recalculate the elev array for this iteration
    if (this%ireadelev == 0) then
      if (this%iform == 1 .or. this%iform == 2) then
        call this%buy_calcelev()
      end if
    end if

buy_cf_bnd()

subroutine gwfbuymodule::buy_cf_bnd ( class(gwfbuytype) this, class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) hnew )

private

Definition at line 300 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    ! -- dummy
    class(GwfBuyType) :: this
    class(BndType), pointer :: packobj
    real(DP), intent(in), dimension(:) :: hnew
    ! -- local
    integer(I4B) :: i, j
    integer(I4B) :: n, locdense, locelev
    integer(I4B), dimension(:), allocatable :: locconc
    !
    ! -- Return if freshwater head formulation; all boundary heads must be
    !    entered as freshwater equivalents
    if (this%iform == 0) return
    !
    ! -- initialize
    locdense = 0
    locelev = 0
    allocate (locconc(this%nrhospecies))
    locconc(:) = 0
    !
    ! -- Add buoyancy terms for head-dependent boundaries
    do n = 1, packobj%naux
      if (packobj%auxname(n) == 'DENSITY') then
        locdense = n
      else if (packobj%auxname(n) == 'ELEVATION') then
        locelev = n
      end if
    end do
    !
    ! -- find aux columns for concentrations that affect density
    do i = 1, this%nrhospecies
      locconc(i) = 0
      do j = 1, packobj%naux
        if (this%cauxspeciesname(i) == packobj%auxname(j)) then
          locconc(i) = j
          exit
        end if
      end do
      if (locconc(i) == 0) then
        ! -- one not found, so don't use and mark all as 0
        locconc(:) = 0
        exit
      end if
    end do
    !
    ! -- Add density terms based on boundary package type
    select case (packobj%filtyp)
    case ('GHB')
      !
      ! -- general head boundary
      call buy_cf_ghb(packobj, this%dense, this%elev, this%denseref, &
                      locelev, locdense, locconc, this%drhodc, this%crhoref, &
                      this%ctemp)
    case ('RIV')
      !
      ! -- river
      call buy_cf_riv(packobj, hnew, this%dense, this%elev, this%denseref, &
                      locelev, locdense, locconc, this%drhodc, this%crhoref, &
                      this%ctemp)
    case ('DRN')
      !
      ! -- drain
      call buy_cf_drn(packobj, this%dense, this%denseref)
    case ('LAK')
      !
      ! -- lake
      call buy_cf_lak(packobj, hnew, this%dense, this%elev, this%denseref, &
                      locdense, locconc, this%drhodc, this%crhoref, &
                      this%ctemp, this%iform)
    case ('SFR')
      !
      ! -- sfr
      call buy_cf_sfr(packobj, hnew, this%dense, this%elev, this%denseref, &
                      locdense, locconc, this%drhodc, this%crhoref, &
                      this%ctemp, this%iform)
    case ('MAW')
      !
      ! -- maw
      call buy_cf_maw(packobj, hnew, this%dense, this%elev, this%denseref, &
                      locdense, locconc, this%drhodc, this%crhoref, &
                      this%ctemp, this%iform)
    case default
      !
      ! -- nothing
    end select
    !
    ! -- deallocate
    deallocate (locconc)

Here is the call graph for this function:

buy_cf_drn()

subroutine gwfbuymodule::buy_cf_drn ( class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) dense, real(dp), intent(in) denseref )

private

Definition at line 601 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    use drnmodule, only: drntype
    class(BndType), pointer :: packobj
    ! -- dummy
    real(DP), intent(in), dimension(:) :: dense
    real(DP), intent(in) :: denseref
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: node
    real(DP) :: rho
    !
    ! -- Process density terms for each DRN
    select type (packobj)
    type is (drntype)
      do n = 1, packobj%nbound
        node = packobj%nodelist(n)
        if (packobj%ibound(node) <= 0) cycle
        ! -- DRN uses the density of the discharging water. Scale the
        !    assembled coefficients so AUXMULTNAME, AUXDEPTHNAME, and the
        !    corresponding discharge elevation apply to the entire flux.
        !    The Newton term for AUXDEPTHNAME in drn_fn is not scaled, so
        !    the Jacobian is approximate, as for the other density terms.
        rho = dense(node) / denseref
        packobj%hcof(n) = rho * packobj%hcof(n)
        packobj%rhs(n) = rho * packobj%rhs(n)
      end do
    end select

Here is the caller graph for this function:

buy_cf_ghb()

subroutine gwfbuymodule::buy_cf_ghb ( class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) dense, real(dp), dimension(:), intent(in) elev, real(dp), intent(in) denseref, integer(i4b), intent(in) locelev, integer(i4b), intent(in) locdense, integer(i4b), dimension(:), intent(in) locconc, real(dp), dimension(:), intent(in) drhodc, real(dp), dimension(:), intent(in) crhoref, real(dp), dimension(:), intent(inout) ctemp )

private

Definition at line 435 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    use ghbmodule, only: ghbtype
    class(BndType), pointer :: packobj
    ! -- dummy
    real(DP), intent(in), dimension(:) :: dense
    real(DP), intent(in), dimension(:) :: elev
    real(DP), intent(in) :: denseref
    integer(I4B), intent(in) :: locelev
    integer(I4B), intent(in) :: locdense
    integer(I4B), dimension(:), intent(in) :: locconc
    real(DP), dimension(:), intent(in) :: drhodc
    real(DP), dimension(:), intent(in) :: crhoref
    real(DP), dimension(:), intent(inout) :: ctemp
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: node
    real(DP) :: denseghb
    real(DP) :: elevghb
    real(DP) :: hghb
    real(DP) :: cond
    real(DP) :: hcofterm, rhsterm
    !
    ! -- Process density terms for each GHB
    select type (packobj)
    type is (ghbtype)
      do n = 1, packobj%nbound
        node = packobj%nodelist(n)
        if (packobj%ibound(node) <= 0) cycle
        !
        ! -- density
        denseghb = get_bnd_density(n, locdense, locconc, denseref, &
                                   drhodc, crhoref, ctemp, packobj%auxvar)
        !
        ! -- elevation
        elevghb = elev(node)
        if (locelev > 0) elevghb = packobj%auxvar(locelev, n)
        !
        ! -- boundary head and conductance
        hghb = packobj%bhead(n)
        cond = packobj%cond_mult(n)
        !
        ! -- calculate HCOF and RHS terms
        call calc_ghb_hcof_rhs_terms(denseref, denseghb, dense(node), &
                                     elevghb, elev(node), hghb, cond, &
                                     rhsterm, hcofterm)
        packobj%hcof(n) = packobj%hcof(n) + hcofterm
        packobj%rhs(n) = packobj%rhs(n) - rhsterm
        !
      end do
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

buy_cf_lak()

subroutine gwfbuymodule::buy_cf_lak ( class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) hnew, real(dp), dimension(:), intent(in) dense, real(dp), dimension(:), intent(in) elev, real(dp), intent(in) denseref, integer(i4b), intent(in) locdense, integer(i4b), dimension(:), intent(in) locconc, real(dp), dimension(:), intent(in) drhodc, real(dp), dimension(:), intent(in) crhoref, real(dp), dimension(:), intent(inout) ctemp, integer(i4b), intent(in) iform )

private

Definition at line 636 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    use lakmodule, only: laktype
    class(BndType), pointer :: packobj
    ! -- dummy
    real(DP), intent(in), dimension(:) :: hnew
    real(DP), intent(in), dimension(:) :: dense
    real(DP), intent(in), dimension(:) :: elev
    real(DP), intent(in) :: denseref
    integer(I4B), intent(in) :: locdense
    integer(I4B), dimension(:), intent(in) :: locconc
    real(DP), dimension(:), intent(in) :: drhodc
    real(DP), dimension(:), intent(in) :: crhoref
    real(DP), dimension(:), intent(inout) :: ctemp
    integer(I4B), intent(in) :: iform
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: node
    real(DP) :: denselak
    !
    ! -- Insert the lake and gwf relative densities into col 1 and 2 and the
    !    gwf elevation into col 3 of the lake package denseterms array
    select type (packobj)
    type is (laktype)
      do n = 1, packobj%nbound
        !
        ! -- get gwf node number
        node = packobj%nodelist(n)
        if (packobj%ibound(node) <= 0) cycle
        !
        ! -- Determine lak density
        denselak = get_bnd_density(n, locdense, locconc, denseref, &
                                   drhodc, crhoref, ctemp, packobj%auxvar)
        !
        ! -- fill lak relative density into column 1 of denseterms
        packobj%denseterms(1, n) = denselak / denseref
        !
        ! -- fill gwf relative density into column 2 of denseterms
        packobj%denseterms(2, n) = dense(node) / denseref
        !
        ! -- fill gwf elevation into column 3 of denseterms
        packobj%denseterms(3, n) = elev(node)
        !
      end do
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

buy_cf_maw()

subroutine gwfbuymodule::buy_cf_maw ( class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) hnew, real(dp), dimension(:), intent(in) dense, real(dp), dimension(:), intent(in) elev, real(dp), intent(in) denseref, integer(i4b), intent(in) locdense, integer(i4b), dimension(:), intent(in) locconc, real(dp), dimension(:), intent(in) drhodc, real(dp), dimension(:), intent(in) crhoref, real(dp), dimension(:), intent(inout) ctemp, integer(i4b), intent(in) iform )

private

Definition at line 742 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    use mawmodule, only: mawtype
    class(BndType), pointer :: packobj
    ! -- dummy
    real(DP), intent(in), dimension(:) :: hnew
    real(DP), intent(in), dimension(:) :: dense
    real(DP), intent(in), dimension(:) :: elev
    real(DP), intent(in) :: denseref
    integer(I4B), intent(in) :: locdense
    integer(I4B), dimension(:), intent(in) :: locconc
    real(DP), dimension(:), intent(in) :: drhodc
    real(DP), dimension(:), intent(in) :: crhoref
    real(DP), dimension(:), intent(inout) :: ctemp
    integer(I4B), intent(in) :: iform
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: node
    real(DP) :: densemaw
    !
    ! -- Insert the maw and gwf relative densities into col 1 and 2 and the
    !    gwf elevation into col 3 of the maw package denseterms array
    select type (packobj)
    type is (mawtype)
      do n = 1, packobj%nbound
        !
        ! -- get gwf node number
        node = packobj%nodelist(n)
        if (packobj%ibound(node) <= 0) cycle
        !
        ! -- Determine maw density
        densemaw = get_bnd_density(n, locdense, locconc, denseref, &
                                   drhodc, crhoref, ctemp, packobj%auxvar)
        !
        ! -- fill maw relative density into column 1 of denseterms
        packobj%denseterms(1, n) = densemaw / denseref
        !
        ! -- fill gwf relative density into column 2 of denseterms
        packobj%denseterms(2, n) = dense(node) / denseref
        !
        ! -- fill gwf elevation into column 3 of denseterms
        packobj%denseterms(3, n) = elev(node)
        !
      end do
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

buy_cf_riv()

subroutine gwfbuymodule::buy_cf_riv ( class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) hnew, real(dp), dimension(:), intent(in) dense, real(dp), dimension(:), intent(in) elev, real(dp), intent(in) denseref, integer(i4b), intent(in) locelev, integer(i4b), intent(in) locdense, integer(i4b), dimension(:), intent(in) locconc, real(dp), dimension(:), intent(in) drhodc, real(dp), dimension(:), intent(in) crhoref, real(dp), dimension(:), intent(inout) ctemp )

private

Definition at line 532 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    use rivmodule, only: rivtype
    class(BndType), pointer :: packobj
    ! -- dummy
    real(DP), intent(in), dimension(:) :: hnew
    real(DP), intent(in), dimension(:) :: dense
    real(DP), intent(in), dimension(:) :: elev
    real(DP), intent(in) :: denseref
    integer(I4B), intent(in) :: locelev
    integer(I4B), intent(in) :: locdense
    integer(I4B), dimension(:), intent(in) :: locconc
    real(DP), dimension(:), intent(in) :: drhodc
    real(DP), dimension(:), intent(in) :: crhoref
    real(DP), dimension(:), intent(inout) :: ctemp
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: node
    real(DP) :: denseriv
    real(DP) :: elevriv
    real(DP) :: hriv
    real(DP) :: rbot
    real(DP) :: cond
    real(DP) :: hcofterm
    real(DP) :: rhsterm
    !
    ! -- Process density terms for each RIV
    select type (packobj)
    type is (rivtype)
      do n = 1, packobj%nbound
        node = packobj%nodelist(n)
        if (packobj%ibound(node) <= 0) cycle
        !
        ! -- density
        denseriv = get_bnd_density(n, locdense, locconc, denseref, &
                                   drhodc, crhoref, ctemp, packobj%auxvar)
        !
        ! -- elevation
        elevriv = elev(node)
        if (locelev > 0) elevriv = packobj%auxvar(locelev, n)
        !
        ! -- boundary head and conductance
        hriv = packobj%stage(n)
        cond = packobj%cond_mult(n)
        rbot = packobj%rbot(n)
        !
        ! -- calculate and add terms depending on whether head is above rbot
        if (hnew(node) > rbot) then
          !
          ! --calculate HCOF and RHS terms, similar to GHB in this case
          call calc_ghb_hcof_rhs_terms(denseref, denseriv, dense(node), &
                                       elevriv, elev(node), hriv, cond, &
                                       rhsterm, hcofterm)
        else
          hcofterm = dzero
          rhsterm = cond * (denseriv / denseref - done) * (hriv - rbot)
        end if
        !
        ! -- Add terms to package hcof and rhs accumulators
        packobj%hcof(n) = packobj%hcof(n) + hcofterm
        packobj%rhs(n) = packobj%rhs(n) - rhsterm
      end do
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

buy_cf_sfr()

subroutine gwfbuymodule::buy_cf_sfr ( class(bndtype), pointer packobj, real(dp), dimension(:), intent(in) hnew, real(dp), dimension(:), intent(in) dense, real(dp), dimension(:), intent(in) elev, real(dp), intent(in) denseref, integer(i4b), intent(in) locdense, integer(i4b), dimension(:), intent(in) locconc, real(dp), dimension(:), intent(in) drhodc, real(dp), dimension(:), intent(in) crhoref, real(dp), dimension(:), intent(inout) ctemp, integer(i4b), intent(in) iform )

private

Definition at line 689 of file gwf-buy.f90.

    ! -- modules
    use bndmodule, only: bndtype
    use sfrmodule, only: sfrtype
    class(BndType), pointer :: packobj
    ! -- dummy
    real(DP), intent(in), dimension(:) :: hnew
    real(DP), intent(in), dimension(:) :: dense
    real(DP), intent(in), dimension(:) :: elev
    real(DP), intent(in) :: denseref
    integer(I4B), intent(in) :: locdense
    integer(I4B), dimension(:), intent(in) :: locconc
    real(DP), dimension(:), intent(in) :: drhodc
    real(DP), dimension(:), intent(in) :: crhoref
    real(DP), dimension(:), intent(inout) :: ctemp
    integer(I4B), intent(in) :: iform
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: node
    real(DP) :: densesfr
    !
    ! -- Insert the sfr and gwf relative densities into col 1 and 2 and the
    !    gwf elevation into col 3 of the sfr package denseterms array
    select type (packobj)
    type is (sfrtype)
      do n = 1, packobj%nbound
        !
        ! -- get gwf node number
        node = packobj%nodelist(n)
        if (packobj%ibound(node) <= 0) cycle
        !
        ! -- Determine sfr density
        densesfr = get_bnd_density(n, locdense, locconc, denseref, &
                                   drhodc, crhoref, ctemp, packobj%auxvar)
        !
        ! -- fill sfr relative density into column 1 of denseterms
        packobj%denseterms(1, n) = densesfr / denseref
        !
        ! -- fill gwf relative density into column 2 of denseterms
        packobj%denseterms(2, n) = dense(node) / denseref
        !
        ! -- fill gwf elevation into column 3 of denseterms
        packobj%denseterms(3, n) = elev(node)
        !
      end do
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

buy_cq()

subroutine gwfbuymodule::buy_cq ( class(gwfbuytype) this, real(dp), dimension(:), intent(in) hnew, real(dp), dimension(:), intent(inout) flowja )

private

Definition at line 870 of file gwf-buy.f90.

    implicit none
    class(GwfBuyType) :: this
    real(DP), intent(in), dimension(:) :: hnew
    real(DP), intent(inout), dimension(:) :: flowja
    integer(I4B) :: n, m, ipos
    real(DP) :: deltaQ
    real(DP) :: rhsterm, amatnn, amatnm
    !
    ! -- Calculate the flow across each cell face and store in flowja
    do n = 1, this%dis%nodes
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        m = this%dis%con%ja(ipos)
        if (m < n) cycle
        if (this%iform == 0) then
          ! -- equivalent freshwater head formulation
          call this%calcbuy(n, m, ipos, hnew(n), hnew(m), deltaq)
        else
          ! -- hydraulic head formulation
          call this%calchhterms(n, m, ipos, hnew(n), hnew(m), rhsterm, &
                                amatnn, amatnm)
          deltaq = amatnm * hnew(m) - amatnn * hnew(n) + rhsterm
        end if
        flowja(ipos) = flowja(ipos) + deltaq
        flowja(this%dis%con%isym(ipos)) = flowja(this%dis%con%isym(ipos)) - &
                                          deltaq
      end do
    end do

buy_cr()

subroutine, public gwfbuymodule::buy_cr	(	type(gwfbuytype), pointer	buyobj,
		character(len=*), intent(in)	name_model,
		character(len=*), intent(in)	input_mempath,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout )

Definition at line 103 of file gwf-buy.f90.

    ! -- dummy
    type(GwfBuyType), pointer :: buyobj
    character(len=*), intent(in) :: name_model
    character(len=*), intent(in) :: input_mempath
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    !
    ! -- Create the object
    allocate (buyobj)
    !
    ! -- create name and memory path
    call buyobj%set_names(1, name_model, 'BUY', 'BUY', input_mempath)
    !
    ! -- Allocate scalars
    call buyobj%allocate_scalars()
    !
    ! -- Set variables
    buyobj%inunit = inunit
    buyobj%iout = iout

Here is the caller graph for this function:

buy_da()

subroutine gwfbuymodule::buy_da ( class(gwfbuytype) this )

private

Definition at line 902 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    !
    ! -- Deallocate arrays if package was active
    if (this%inunit > 0) then
      call mem_deallocate(this%elev)
      call mem_deallocate(this%dense)
      call mem_deallocate(this%concbuy)
      call mem_deallocate(this%drhodc)
      call mem_deallocate(this%crhoref)
      call mem_deallocate(this%ctemp)
      deallocate (this%cmodelname)
      deallocate (this%cauxspeciesname)
      deallocate (this%modelconc)
    end if
    !
    ! -- Scalars
    call mem_deallocate(this%ioutdense)
    call mem_deallocate(this%iform)
    call mem_deallocate(this%ireadelev)
    call mem_deallocate(this%ireadconcbuy)
    call mem_deallocate(this%iconcset)
    call mem_deallocate(this%denseref)
    !
    call mem_deallocate(this%nrhospecies)
    !
    ! -- deallocate parent
    call this%NumericalPackageType%da()

buy_df()

subroutine gwfbuymodule::buy_df ( class(gwfbuytype) this, class(disbasetype), intent(in), pointer dis, type(gwfbuyinputdatatype), intent(in), optional buy_input )

private

Parameters

	this	this buoyancy package
[in]	dis	pointer to discretization
[in]	buy_input	optional buy input data, otherwise read from file

Definition at line 127 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this !< this buoyancy package
    class(DisBaseType), pointer, intent(in) :: dis !< pointer to discretization
    type(GwfBuyInputDataType), optional, intent(in) :: buy_input !< optional buy input data, otherwise read from file
    ! -- local
    ! -- formats
    character(len=*), parameter :: fmtbuy = &
      "(1x,/1x,'BUY -- Buoyancy Package, Version 1, 5/16/2018', &
      &' input read from mempath: ', a, /)"
    !
    ! --print a message identifying the buoyancy package.
    if (.not. present(buy_input)) then
      write (this%iout, fmtbuy) this%input_mempath
    end if
    !
    ! -- store pointers to arguments that were passed in
    this%dis => dis
 
    if (.not. present(buy_input)) then
      !
      ! -- Read buoyancy options
      call this%source_options()
      !
      ! -- Read buoyancy dimensions
      call this%source_dimensions()
    else
      ! set from input data instead
      call this%set_options(buy_input)
      this%nrhospecies = buy_input%nrhospecies
    end if
    !
    ! -- Allocate arrays
    call this%allocate_arrays(dis%nodes)
 
    if (.not. present(buy_input)) then
      !
      ! -- Read buoyancy packagedata
      call this%source_packagedata()
    else
      ! set from input data instead
      call this%set_packagedata(buy_input)
    end if

buy_fc()

subroutine gwfbuymodule::buy_fc ( class(gwfbuytype) this, integer(i4b) kiter, class(matrixbasetype), pointer matrix_sln, integer(i4b), dimension(:), intent(in) idxglo, real(dp), dimension(:), intent(inout) rhs, real(dp), dimension(:), intent(inout) hnew )

private

Definition at line 793 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    integer(I4B) :: kiter
    class(MatrixBaseType), pointer :: matrix_sln
    integer(I4B), intent(in), dimension(:) :: idxglo
    real(DP), dimension(:), intent(inout) :: rhs
    real(DP), intent(inout), dimension(:) :: hnew
    ! -- local
    integer(I4B) :: n, m, ipos, idiag
    real(DP) :: rhsterm, amatnn, amatnm
    !
    ! -- initialize
    amatnn = dzero
    amatnm = dzero
    !
    ! -- fill buoyancy flow term
    do n = 1, this%dis%nodes
      if (this%ibound(n) == 0) cycle
      idiag = this%dis%con%ia(n)
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) == 0) cycle
        if (this%iform == 0) then
          call this%calcbuy(n, m, ipos, hnew(n), hnew(m), rhsterm)
        else
          call this%calchhterms(n, m, ipos, hnew(n), hnew(m), rhsterm, &
                                amatnn, amatnm)
        end if
        !
        ! -- Add terms to rhs, diagonal, and off diagonal
        rhs(n) = rhs(n) - rhsterm
        call matrix_sln%add_value_pos(idxglo(idiag), -amatnn)
        call matrix_sln%add_value_pos(idxglo(ipos), amatnm)
      end do
    end do

buy_ot_dv()

subroutine gwfbuymodule::buy_ot_dv ( class(gwfbuytype) this, integer(i4b), intent(in) idvfl )

private

Definition at line 833 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    integer(I4B), intent(in) :: idvfl
    ! -- local
    character(len=1) :: cdatafmp = ' ', editdesc = ' '
    integer(I4B) :: ibinun
    integer(I4B) :: iprint
    integer(I4B) :: nvaluesp
    integer(I4B) :: nwidthp
    real(DP) :: dinact
    !
    ! -- Set unit number for density output
    if (this%ioutdense /= 0) then
      ibinun = 1
    else
      ibinun = 0
    end if
    if (idvfl == 0) ibinun = 0
    !
    ! -- save density array
    if (ibinun /= 0) then
      iprint = 0
      dinact = dhnoflo
      !
      ! -- write density to binary file
      if (this%ioutdense /= 0) then
        ibinun = this%ioutdense
        call this%dis%record_array(this%dense, this%iout, iprint, ibinun, &
                                   '         DENSITY', cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      end if
    end if

buy_rp()

subroutine gwfbuymodule::buy_rp ( class(gwfbuytype) this )

private

Definition at line 244 of file gwf-buy.f90.

    ! -- modules
    use tdismodule, only: kstp, kper
    ! -- dummy
    class(GwfBuyType) :: this
    ! -- local
    character(len=LINELENGTH) :: errmsg
    integer(I4B) :: i
    ! -- formats
    character(len=*), parameter :: fmtc = &
      "('Buoyancy Package does not have a concentration set &
       &for species ',i0,'. One or more model names may be specified &
       &incorrectly in the PACKAGEDATA block or a gwf-gwt exchange may need &
       &to be activated.')"
    !
    ! -- Check to make sure all concentration pointers have been set
    if (kstp * kper == 1) then
      do i = 1, this%nrhospecies
        if (.not. associated(this%modelconc(i)%conc)) then
          write (errmsg, fmtc) i
          call store_error(errmsg)
        end if
      end do
      if (count_errors() > 0) then
        call store_error_filename(this%input_fname)
      end if
    end if

Here is the call graph for this function:

calc_ghb_hcof_rhs_terms()

subroutine gwfbuymodule::calc_ghb_hcof_rhs_terms ( real(dp), intent(in) denseref, real(dp), intent(in) denseghb, real(dp), intent(in) densenode, real(dp), intent(in) elevghb, real(dp), intent(in) elevnode, real(dp), intent(in) hghb, real(dp), intent(in) cond, real(dp), intent(inout) rhsterm, real(dp), intent(inout) hcofterm )

private

These terms only affect the diagonal and cannot make the matrix asymmetric, so they are added to the matrix for both hydraulic head formulations. Lagging them on the right-hand side for HHFORMULATION_RHS caused outer iterations to cycle when a RIV switched across rbot.

Definition at line 497 of file gwf-buy.f90.

    ! -- dummy
    real(DP), intent(in) :: denseref
    real(DP), intent(in) :: denseghb
    real(DP), intent(in) :: densenode
    real(DP), intent(in) :: elevghb
    real(DP), intent(in) :: elevnode
    real(DP), intent(in) :: hghb
    real(DP), intent(in) :: cond
    real(DP), intent(inout) :: rhsterm
    real(DP), intent(inout) :: hcofterm
    ! -- local
    real(DP) :: t1, t2
    real(DP) :: avgdense, avgelev
    !
    ! -- Calculate common terms
    avgdense = dhalf * denseghb + dhalf * densenode
    avgelev = dhalf * elevghb + dhalf * elevnode
    t1 = avgdense / denseref - done
    t2 = (denseghb - densenode) / denseref
    !
    ! -- Add hcof terms
    hcofterm = -cond * t1
    hcofterm = hcofterm + dhalf * cond * t2
    !
    ! -- Add rhs terms
    rhsterm = cond * t1 * hghb
    rhsterm = rhsterm - cond * t2 * avgelev
    rhsterm = rhsterm + dhalf * cond * t2 * hghb

Here is the caller graph for this function:

calcbuy()

subroutine gwfbuymodule::calcbuy ( class(gwfbuytype) this, integer(i4b), intent(in) n, integer(i4b), intent(in) m, integer(i4b), intent(in) icon, real(dp), intent(in) hn, real(dp), intent(in) hm, real(dp), intent(inout) buy )

private

Definition at line 1073 of file gwf-buy.f90.

    ! -- modules
    use gwfconductanceutilsmodule, only: hcond, vcond
    ! -- dummy
    class(GwfBuyType) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: m
    integer(I4B), intent(in) :: icon
    real(DP), intent(in) :: hn
    real(DP), intent(in) :: hm
    real(DP), intent(inout) :: buy
    ! -- local
    integer(I4B) :: ihc
    real(DP) :: densen, densem, cl1, cl2, avgdense, wt, elevn, elevm, &
                cond, tp, bt
    real(DP) :: hyn
    real(DP) :: hym
    !
    ! -- Average density
    densen = this%dense(n)
    densem = this%dense(m)
    if (m > n) then
      cl1 = this%dis%con%cl1(this%dis%con%jas(icon))
      cl2 = this%dis%con%cl2(this%dis%con%jas(icon))
    else
      cl1 = this%dis%con%cl2(this%dis%con%jas(icon))
      cl2 = this%dis%con%cl1(this%dis%con%jas(icon))
    end if
    wt = cl1 / (cl1 + cl2)
    avgdense = wt * densen + (done - wt) * densem
    !
    ! -- Elevations
    if (this%ireadelev == 0) then
      tp = this%dis%top(n)
      bt = this%dis%bot(n)
      elevn = bt + dhalf * this%npf%sat(n) * (tp - bt)
      tp = this%dis%top(m)
      bt = this%dis%bot(m)
      elevm = bt + dhalf * this%npf%sat(m) * (tp - bt)
    else
      elevn = this%elev(n)
      elevm = this%elev(m)
    end if
    !
    ihc = this%dis%con%ihc(this%dis%con%jas(icon))
    hyn = this%npf%hy_eff(n, m, ihc, ipos=icon)
    hym = this%npf%hy_eff(m, n, ihc, ipos=icon)
    !
    ! -- Conductance
    if (this%dis%con%ihc(this%dis%con%jas(icon)) == 0) then
      cond = vcond(this%ibound(n), this%ibound(m), &
                   this%npf%icelltype(n), this%npf%icelltype(m), &
                   this%npf%inewton, &
                   this%npf%ivarcv, this%npf%idewatcv, &
                   this%npf%condsat(this%dis%con%jas(icon)), hn, hm, &
                   hyn, hym, &
                   this%npf%sat(n), this%npf%sat(m), &
                   this%dis%top(n), this%dis%top(m), &
                   this%dis%bot(n), this%dis%bot(m), &
                   this%dis%con%hwva(this%dis%con%jas(icon)))
    else
      cond = hcond(this%ibound(n), this%ibound(m), &
                   this%npf%icelltype(n), this%npf%icelltype(m), &
                   this%npf%inewton, &
                   this%dis%con%ihc(this%dis%con%jas(icon)), &
                   this%npf%icellavg, &
                   this%npf%condsat(this%dis%con%jas(icon)), &
                   hn, hm, this%npf%sat(n), this%npf%sat(m), &
                   hyn, hym, &
                   this%dis%top(n), this%dis%top(m), &
                   this%dis%bot(n), this%dis%bot(m), &
                   this%dis%con%cl1(this%dis%con%jas(icon)), &
                   this%dis%con%cl2(this%dis%con%jas(icon)), &
                   this%dis%con%hwva(this%dis%con%jas(icon)))
    end if
    !
    ! -- Calculate buoyancy term
    buy = cond * (avgdense - this%denseref) / this%denseref * (elevm - elevn)

Here is the call graph for this function:

calcdens()

real(dp) function gwfbuymodule::calcdens	(	real(dp), intent(in)	denseref,
		real(dp), dimension(:), intent(in)	drhodc,
		real(dp), dimension(:), intent(in)	crhoref,
		real(dp), dimension(:), intent(in)	conc )

Definition at line 82 of file gwf-buy.f90.

    ! -- dummy
    real(DP), intent(in) :: denseref
    real(DP), dimension(:), intent(in) :: drhodc
    real(DP), dimension(:), intent(in) :: crhoref
    real(DP), dimension(:), intent(in) :: conc
    ! -- return
    real(DP) :: dense
    ! -- local
    integer(I4B) :: nrhospec
    integer(I4B) :: i
    !
    nrhospec = size(drhodc)
    dense = denseref
    do i = 1, nrhospec
      dense = dense + drhodc(i) * (conc(i) - crhoref(i))
    end do

Here is the caller graph for this function:

calchhterms()

subroutine gwfbuymodule::calchhterms ( class(gwfbuytype) this, integer(i4b), intent(in) n, integer(i4b), intent(in) m, integer(i4b), intent(in) icon, real(dp), intent(in) hn, real(dp), intent(in) hm, real(dp), intent(inout) rhsterm, real(dp), intent(inout) amatnn, real(dp), intent(inout) amatnm )

private

Definition at line 1155 of file gwf-buy.f90.

    ! -- modules
    use gwfconductanceutilsmodule, only: hcond, vcond
    ! -- dummy
    class(GwfBuyType) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: m
    integer(I4B), intent(in) :: icon
    real(DP), intent(in) :: hn
    real(DP), intent(in) :: hm
    real(DP), intent(inout) :: rhsterm
    real(DP), intent(inout) :: amatnn
    real(DP), intent(inout) :: amatnm
    ! -- local
    integer(I4B) :: ihc
    real(DP) :: densen, densem, cl1, cl2, avgdense, wt, elevn, elevm, cond
    real(DP) :: rhonormn, rhonormm
    real(DP) :: rhoterm
    real(DP) :: elevnm
    real(DP) :: hphi
    real(DP) :: hyn
    real(DP) :: hym
    !
    ! -- Average density
    densen = this%dense(n)
    densem = this%dense(m)
    if (m > n) then
      cl1 = this%dis%con%cl1(this%dis%con%jas(icon))
      cl2 = this%dis%con%cl2(this%dis%con%jas(icon))
    else
      cl1 = this%dis%con%cl2(this%dis%con%jas(icon))
      cl2 = this%dis%con%cl1(this%dis%con%jas(icon))
    end if
    wt = cl1 / (cl1 + cl2)
    avgdense = wt * densen + (1.0 - wt) * densem
    !
    ! -- Elevations
    elevn = this%elev(n)
    elevm = this%elev(m)
    elevnm = (done - wt) * elevn + wt * elevm
    !
    ihc = this%dis%con%ihc(this%dis%con%jas(icon))
    hyn = this%npf%hy_eff(n, m, ihc, ipos=icon)
    hym = this%npf%hy_eff(m, n, ihc, ipos=icon)
    !
    ! -- Conductance
    if (ihc == 0) then
      cond = vcond(this%ibound(n), this%ibound(m), &
                   this%npf%icelltype(n), this%npf%icelltype(m), &
                   this%npf%inewton, &
                   this%npf%ivarcv, this%npf%idewatcv, &
                   this%npf%condsat(this%dis%con%jas(icon)), hn, hm, &
                   hyn, hym, &
                   this%npf%sat(n), this%npf%sat(m), &
                   this%dis%top(n), this%dis%top(m), &
                   this%dis%bot(n), this%dis%bot(m), &
                   this%dis%con%hwva(this%dis%con%jas(icon)))
    else
      cond = hcond(this%ibound(n), this%ibound(m), &
                   this%npf%icelltype(n), this%npf%icelltype(m), &
                   this%npf%inewton, &
                   this%dis%con%ihc(this%dis%con%jas(icon)), &
                   this%npf%icellavg, &
                   this%npf%condsat(this%dis%con%jas(icon)), &
                   hn, hm, this%npf%sat(n), this%npf%sat(m), &
                   hyn, hym, &
                   this%dis%top(n), this%dis%top(m), &
                   this%dis%bot(n), this%dis%bot(m), &
                   this%dis%con%cl1(this%dis%con%jas(icon)), &
                   this%dis%con%cl2(this%dis%con%jas(icon)), &
                   this%dis%con%hwva(this%dis%con%jas(icon)))
    end if
    !
    ! -- Calculate terms
    rhonormn = densen / this%denseref
    rhonormm = densem / this%denseref
    rhoterm = wt * rhonormn + (done - wt) * rhonormm
    amatnn = cond * (rhoterm - done)
    amatnm = amatnn
    rhsterm = -cond * (rhonormm - rhonormn) * elevnm
    if (this%iform == 1) then
      ! -- rhs (lag the h terms and keep matrix symmetric)
      hphi = (done - wt) * hn + wt * hm
      rhsterm = rhsterm + cond * hphi * (rhonormm - rhonormn)
    else
      ! -- lhs, results in asymmetric matrix due to weight term
      amatnn = amatnn - cond * (done - wt) * (rhonormm - rhonormn)
      amatnm = amatnm + cond * wt * (rhonormm - rhonormn)
    end if

Here is the call graph for this function:

get_bnd_density()

real(dp) function gwfbuymodule::get_bnd_density ( integer(i4b), intent(in) n, integer(i4b), intent(in) locdense, integer(i4b), dimension(:), intent(in) locconc, real(dp), intent(in) denseref, real(dp), dimension(:), intent(in) drhodc, real(dp), dimension(:), intent(in) crhoref, real(dp), dimension(:), intent(inout) ctemp, real(dp), dimension(:, :), intent(in) auxvar )

private

1. Assign as aux variable in column with name 'DENSITY'
2. Calculate using equation of state and nrhospecies aux columns
3. If neither of those, then assign as denseref

Definition at line 398 of file gwf-buy.f90.

    ! -- dummy
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: locdense
    integer(I4B), dimension(:), intent(in) :: locconc
    real(DP), intent(in) :: denseref
    real(DP), dimension(:), intent(in) :: drhodc
    real(DP), dimension(:), intent(in) :: crhoref
    real(DP), dimension(:), intent(inout) :: ctemp
    real(DP), dimension(:, :), intent(in) :: auxvar
    ! -- return
    real(DP) :: densebnd
    ! -- local
    integer(I4B) :: i
    !
    ! -- assign boundary density based on one of three options
    if (locdense > 0) then
      ! -- assign density to an aux column named 'DENSITY'
      densebnd = auxvar(locdense, n)
    else if (locconc(1) > 0) then
      ! -- calculate density using one or more concentration auxcolumns
      do i = 1, size(locconc)
        ctemp(i) = dzero
        if (locconc(i) > 0) then
          ctemp(i) = auxvar(locconc(i), n)
        end if
      end do
      densebnd = calcdens(denseref, drhodc, crhoref, ctemp)
    else
      ! -- neither of the above, so assign as denseref
      densebnd = denseref
    end if

Here is the call graph for this function:

Here is the caller graph for this function:

log_options()

subroutine gwfbuymodule::log_options ( class(gwfbuytype), intent(inout) this, type(gwfbuyparamfoundtype), intent(in) found, character(len=*), intent(in) densityfile )

private

Definition at line 1407 of file gwf-buy.f90.

    ! -- modules
    use gwfbuyinputmodule, only: gwfbuyparamfoundtype
    ! -- dummy variables
    class(GwfBuyType), intent(inout) :: this
    type(GwfBuyParamFoundType), intent(in) :: found
    character(len=*), intent(in) :: densityfile
    ! -- local variables
    ! -- formats
    character(len=*), parameter :: fmtfileout = &
      "(4x, 'BUY ', 1x, a, 1x, ' will be saved to file: ', &
      &a, ' opened on unit: ', I7)"
    !
    write (this%iout, '(1x,a)') 'Processing BUY OPTIONS block'
 
    if (found%hhform_rhs) then
      write (this%iout, '(4x,a)') &
        'Hydraulic head formulation set to right-hand side'
    end if
    if (found%denseref) then
      write (this%iout, '(4x,a,1pg15.6)') &
        'Reference density has been set to: ', this%denseref
    end if
    if (found%dev_efh_form) then
      write (this%iout, '(4x,a)') &
        'Formulation set to equivalent freshwater head'
    end if
    if (found%densityfile) then
      write (this%iout, fmtfileout) &
        'DENSITY', trim(densityfile), this%ioutdense
    end if
 
    write (this%iout, '(1x,a)') 'End of BUY OPTIONS block'

set_concentration_pointer()

subroutine gwfbuymodule::set_concentration_pointer ( class(gwfbuytype) this, character(len=lenmodelname), intent(in) modelname, real(dp), dimension(:), pointer conc, integer(i4b), dimension(:), pointer icbund )

private

This routine is called from the gwfgwt exchange in the exg_ar() method

Definition at line 1465 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    character(len=LENMODELNAME), intent(in) :: modelname
    real(DP), dimension(:), pointer :: conc
    integer(I4B), dimension(:), pointer :: icbund
    ! -- local
    integer(I4B) :: i
    logical :: found
    !
    this%iconcset = 1
    found = .false.
    do i = 1, this%nrhospecies
      if (this%cmodelname(i) == modelname) then
        this%modelconc(i)%conc => conc
        this%modelconc(i)%icbund => icbund
        found = .true.
        exit
      end if
    end do

set_options()

subroutine gwfbuymodule::set_options ( class(gwfbuytype) this, type(gwfbuyinputdatatype), intent(in) input_data )

private

Parameters

[in]

input_data

the input data to be set

Definition at line 1444 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this
    type(GwfBuyInputDataType), intent(in) :: input_data !< the input data to be set
    !
    this%iform = input_data%iform
    this%denseref = input_data%denseref
    !
    ! derived option:
    ! if not iform==2, there is no asymmetry
    if (this%iform == 0 .or. this%iform == 1) then
      this%iasym = 0
    end if

set_packagedata()

subroutine gwfbuymodule::set_packagedata ( class(gwfbuytype) this, type(gwfbuyinputdatatype), intent(in) input_data )

private

Parameters

	this	this buyoancy pkg
[in]	input_data	the input data to be set

Definition at line 1056 of file gwf-buy.f90.

    ! -- dummy
    class(GwfBuyType) :: this !< this buyoancy pkg
    type(GwfBuyInputDataType), intent(in) :: input_data !< the input data to be set
    ! -- local
    integer(I4B) :: ispec
 
    do ispec = 1, this%nrhospecies
      this%drhodc(ispec) = input_data%drhodc(ispec)
      this%crhoref(ispec) = input_data%crhoref(ispec)
      this%cmodelname(ispec) = input_data%cmodelname(ispec)
      this%cauxspeciesname(ispec) = input_data%cauxspeciesname(ispec)
    end do

source_dimensions()

subroutine gwfbuymodule::source_dimensions ( class(gwfbuytype), intent(inout) this )

private

Definition at line 935 of file gwf-buy.f90.

    ! -- modules
    use memorymanagerextmodule, only: mem_set_value
    use gwfbuyinputmodule, only: gwfbuyparamfoundtype
    ! -- dummy
    class(GwfBuyType), intent(inout) :: this
    ! -- local
    type(GwfBuyParamFoundType) :: found
 
    ! update defaults from input context
    call mem_set_value(this%nrhospecies, 'NRHOSPECIES', this%input_mempath, &
                       found%nrhospecies)
 
    write (this%iout, '(/1x,a)') 'Processing BUY DIMENSIONS block'
    write (this%iout, '(4x,a,i0)') 'NRHOSPECIES = ', this%nrhospecies
    write (this%iout, '(1x,a)') 'End of BUY DIMENSIONS block'
 
    ! check dimension
    if (this%nrhospecies < 1) then
      call store_error('NRHOSPECIES must be greater than zero.')
      call store_error_filename(this%input_fname)
    end if

Here is the call graph for this function:

source_options()

subroutine gwfbuymodule::source_options ( class(gwfbuytype), intent(inout) this )

private

Definition at line 1362 of file gwf-buy.f90.

    ! -- modules
    use memorymanagerextmodule, only: mem_set_value
    use openspecmodule, only: access, form
    use inputoutputmodule, only: getunit, openfile
    use gwfbuyinputmodule, only: gwfbuyparamfoundtype
    ! -- dummy
    class(GwfBuyType), intent(inout) :: this
    ! -- local
    character(len=LINELENGTH) :: densityfile
    type(GwfBuyParamFoundType) :: found
 
    ! initialize variables
    densityfile = ''
 
    ! update defaults from input context
    call mem_set_value(this%iform, 'HHFORM_RHS', this%input_mempath, &
                       found%hhform_rhs)
    call mem_set_value(this%denseref, 'DENSEREF', this%input_mempath, &
                       found%denseref)
    call mem_set_value(this%iform, 'DEV_EFH_FORM', this%input_mempath, &
                       found%dev_efh_form)
    call mem_set_value(densityfile, 'DENSITYFILE', this%input_mempath, &
                       found%densityfile)
 
    ! update input dependent internal state
    if (found%hhform_rhs) this%iasym = 0
    if (found%dev_efh_form) then
      this%iform = 0
      this%iasym = 0
    end if
 
    ! fileout options
    if (found%densityfile) then
      this%ioutdense = getunit()
      call openfile(this%ioutdense, this%iout, densityfile, 'DATA(BINARY)', &
                    form, access, 'REPLACE')
    end if
 
    ! log options
    call this%log_options(found, densityfile)

Here is the call graph for this function:

source_packagedata()

subroutine gwfbuymodule::source_packagedata ( class(gwfbuytype), intent(inout) this )

private

Definition at line 961 of file gwf-buy.f90.

    ! -- modules
    use memorymanagermodule, only: mem_setptr, mem_allocate
    use memorymanagerextmodule, only: mem_set_value, memorystore_release
    use characterstringmodule, only: characterstringtype
    ! -- dummy
    class(GwfBuyType), intent(inout) :: this
    ! -- local
    integer(I4B), dimension(:), pointer, contiguous :: irhospec
    type(CharacterStringType), dimension(:), pointer, &
      contiguous :: modelnames, auxspeciesnames
    real(DP), dimension(:), pointer, contiguous :: drhodc, crhoref
    integer(I4B), dimension(:), allocatable :: itemp
    character(len=LINELENGTH) :: modelname, auxspeciesname, line, errmsg
    character(len=10) :: c10
    character(len=16) :: c16
    integer(I4B) :: n
    ! -- format
    character(len=*), parameter :: fmterr = &
      "('Invalid value for IRHOSPEC (',i0,') detected in BUY Package. &
      &IRHOSPEC must be > 0 and <= NRHOSPECIES, and duplicate values &
      &are not allowed.')"
 
    ! initialize
    allocate (itemp(this%nrhospecies))
    itemp(:) = 0
 
    ! set input context pointers
    call mem_setptr(irhospec, 'IRHOSPEC', this%input_mempath)
    call mem_setptr(drhodc, 'DRHODC', this%input_mempath)
    call mem_setptr(crhoref, 'CRHOREF', this%input_mempath)
    call mem_setptr(modelnames, 'MODELNAME', this%input_mempath)
    call mem_setptr(auxspeciesnames, 'AUXSPECIESNAME', this%input_mempath)
 
    ! process package data
    do n = 1, size(irhospec)
      modelname = modelnames(n)
      auxspeciesname = auxspeciesnames(n)
 
      if (irhospec(n) < 1 .or. irhospec(n) > this%nrhospecies) then
        write (errmsg, fmterr) irhospec(n)
        call store_error(errmsg)
      end if
      if (itemp(irhospec(n)) /= 0) then
        write (errmsg, fmterr) irhospec(n)
        call store_error(errmsg)
      end if
      itemp(irhospec(n)) = 1
 
      this%drhodc(irhospec(n)) = drhodc(n)
      this%crhoref(irhospec(n)) = crhoref(n)
      this%cmodelname(irhospec(n)) = trim(modelname)
      this%cauxspeciesname(irhospec(n)) = trim(auxspeciesname)
    end do
 
    ! Check for errors.
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
 
    ! log package data
    write (this%iout, '(/,1x,a)') 'Processing BUY PACKAGEDATA block'
 
    ! write packagedata information
    write (this%iout, '(1x,a)') 'Summary of species information in BUY Package'
    write (this%iout, '(1a11, 4a17)') &
      'SPECIES', 'DRHODC', 'CRHOREF', 'MODEL', &
      'AUXSPECIESNAME'
    do n = 1, this%nrhospecies
      write (c10, '(i0)') n
      line = ' '//adjustr(c10)
      write (c16, '(g15.6)') this%drhodc(n)
      line = trim(line)//' '//adjustr(c16)
      write (c16, '(g15.6)') this%crhoref(n)
      line = trim(line)//' '//adjustr(c16)
      write (c16, '(a)') this%cmodelname(n)
      line = trim(line)//' '//adjustr(c16)
      write (c16, '(a)') this%cauxspeciesname(n)
      line = trim(line)//' '//adjustr(c16)
      write (this%iout, '(a)') trim(line)
    end do
 
    write (this%iout, '(1x,a)') 'End of BUY PACKAGEDATA block'
 
    ! cleanup
    deallocate (itemp)
    call memorystore_release('IRHOSPEC', this%input_mempath)
    call memorystore_release('DRHODC', this%input_mempath)
    call memorystore_release('CRHOREF', this%input_mempath)
    call memorystore_release('MODELNAME', this%input_mempath)
    call memorystore_release('AUXSPECIESNAME', this%input_mempath)

Here is the call graph for this function: