gwtmstmodule Module Reference

– @ brief Mobile Storage and Transfer (MST) Module More...

Data Types
type	gwtmsttype
	@ brief Mobile storage and transfer More...

Functions/Subroutines
subroutine, public	mst_cr (mstobj, name_model, inunit, iout, fmi)
	@ brief Create a new package object More...
subroutine	mst_ar (this, dis, ibound)
	@ brief Allocate and read method for package More...
subroutine	mst_fc (this, nodes, cold, nja, matrix_sln, idxglo, cnew, rhs, kiter)
	@ brief Fill coefficient method for package More...
subroutine	mst_fc_sto (this, nodes, cold, nja, matrix_sln, idxglo, rhs)
	@ brief Fill storage coefficient method for package More...
subroutine	mst_fc_dcy (this, nodes, cold, cnew, nja, matrix_sln, idxglo, rhs, kiter)
	@ brief Fill decay coefficient method for package More...
subroutine	mst_fc_srb (this, nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
	@ brief Fill sorption coefficient method for package More...
subroutine	mst_srb_term (isrb, volfracm, rhobm, vcell, tled, cnew, cold, swnew, swold, const1, const2, rate, hcofval, rhsval)
	@ brief Calculate sorption terms More...
subroutine	mst_fc_dcy_srb (this, nodes, cold, nja, matrix_sln, idxglo, rhs, cnew, kiter)
	@ brief Fill sorption-decay coefficient method for package More...
subroutine	mst_cq (this, nodes, cnew, cold, flowja)
	@ brief Calculate flows for package More...
subroutine	mst_cq_sto (this, nodes, cnew, cold, flowja)
	@ brief Calculate storage terms for package More...
subroutine	mst_cq_dcy (this, nodes, cnew, cold, flowja)
	@ brief Calculate decay terms for package More...
subroutine	mst_cq_srb (this, nodes, cnew, cold, flowja)
	@ brief Calculate sorption terms for package More...
subroutine	mst_cq_dcy_srb (this, nodes, cnew, cold, flowja)
	@ brief Calculate decay-sorption terms for package More...
subroutine	mst_bd (this, isuppress_output, model_budget)
	@ brief Calculate budget terms for package More...
subroutine	mst_ot_flow (this, icbcfl, icbcun)
	@ brief Output flow terms for package More...
subroutine	mst_da (this)
	@ brief Deallocate More...
subroutine	allocate_scalars (this)
	@ brief Allocate scalar variables for package More...
subroutine	allocate_arrays (this, nodes)
	@ brief Allocate arrays for package More...
subroutine	read_options (this)
	@ brief Read options for package More...
subroutine	read_data (this)
	@ brief Read data for package More...
subroutine	addto_volfracim (this, volfracim)
	@ brief Add volfrac values to volfracim More...
real(dp) function	get_volfracm (this, node)
	@ brief Return mobile domain volume fraction More...
real(dp) function	get_freundlich_conc (conc, kf, a)
	@ brief Calculate sorption concentration using Freundlich More...
real(dp) function	get_langmuir_conc (conc, kl, sbar)
	@ brief Calculate sorption concentration using Langmuir More...
real(dp) function	get_freundlich_derivative (conc, kf, a)
	@ brief Calculate sorption derivative using Freundlich More...
real(dp) function	get_langmuir_derivative (conc, kl, sbar)
	@ brief Calculate sorption derivative using Langmuir More...
real(dp) function	get_zero_order_decay (decay_rate_usr, decay_rate_last, kiter, cold, cnew, delt)
	@ brief Calculate zero-order decay rate and constrain if necessary More...

Variables
integer(i4b), parameter	nbditems = 4
character(len=lenbudtxt), dimension(nbditems)	budtxt

Detailed Description

The GwtMstModule contains the GwtMstType, which is the derived type responsible for adding the effects of

1. Changes in dissolved solute mass
2. Decay of dissolved solute mass
3. Sorption
4. Decay of sorbed solute mass

Function/Subroutine Documentation

addto_volfracim()

subroutine gwtmstmodule::addto_volfracim	(	class(gwtmsttype)	this,
		real(dp), dimension(:), intent(in)	volfracim
	)

Method to add immobile domain volume fracions, which are stored as a cumulative value in volfracim.

Parameters

	this	GwtMstType object
[in]	volfracim	immobile domain volume fraction that contributes to total immobile volume fraction

Definition at line 1419 of file gwt-mst.f90.

    ! -- modules
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    real(DP), dimension(:), intent(in) :: volfracim !< immobile domain volume fraction that contributes to total immobile volume fraction
    ! -- local
    integer(I4B) :: n
    !
    ! -- Add to volfracim
    do n = 1, this%dis%nodes
      this%volfracim(n) = this%volfracim(n) + volfracim(n)
    end do
    !
    ! -- An immobile domain is adding a volume fraction, so update thetam
    !    accordingly.
    do n = 1, this%dis%nodes
      this%thetam(n) = this%get_volfracm(n) * this%porosity(n)
    end do
    !
    ! -- Return
    return

allocate_arrays()

subroutine gwtmstmodule::allocate_arrays	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	nodes
	)

Method to allocate arrays for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes

Definition at line 1051 of file gwt-mst.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    use constantsmodule, only: dzero
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: nodes !< number of nodes
    ! -- local
    integer(I4B) :: n
    !
    ! -- Allocate
    ! -- sto
    call mem_allocate(this%porosity, nodes, 'POROSITY', this%memoryPath)
    call mem_allocate(this%thetam, nodes, 'THETAM', this%memoryPath)
    call mem_allocate(this%volfracim, nodes, 'VOLFRACIM', this%memoryPath)
    call mem_allocate(this%ratesto, nodes, 'RATESTO', this%memoryPath)
    !
    ! -- dcy
    if (this%idcy == 0) then
      call mem_allocate(this%ratedcy, 1, 'RATEDCY', this%memoryPath)
      call mem_allocate(this%decay, 1, 'DECAY', this%memoryPath)
      call mem_allocate(this%decaylast, 1, 'DECAYLAST', this%memoryPath)
    else
      call mem_allocate(this%ratedcy, this%dis%nodes, 'RATEDCY', this%memoryPath)
      call mem_allocate(this%decay, nodes, 'DECAY', this%memoryPath)
      call mem_allocate(this%decaylast, nodes, 'DECAYLAST', this%memoryPath)
    end if
    if (this%idcy /= 0 .and. this%isrb /= 0) then
      call mem_allocate(this%ratedcys, this%dis%nodes, 'RATEDCYS', &
                        this%memoryPath)
      call mem_allocate(this%decayslast, this%dis%nodes, 'DECAYSLAST', &
                        this%memoryPath)
    else
      call mem_allocate(this%ratedcys, 1, 'RATEDCYS', this%memoryPath)
      call mem_allocate(this%decayslast, 1, 'DECAYSLAST', this%memoryPath)
    end if
    call mem_allocate(this%decay_sorbed, 1, 'DECAY_SORBED', &
                      this%memoryPath)
    !
    ! -- srb
    if (this%isrb == 0) then
      call mem_allocate(this%bulk_density, 1, 'BULK_DENSITY', this%memoryPath)
      call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
      call mem_allocate(this%distcoef, 1, 'DISTCOEF', this%memoryPath)
      call mem_allocate(this%ratesrb, 1, 'RATESRB', this%memoryPath)
    else
      call mem_allocate(this%bulk_density, nodes, 'BULK_DENSITY', this%memoryPath)
      call mem_allocate(this%distcoef, nodes, 'DISTCOEF', this%memoryPath)
      call mem_allocate(this%ratesrb, nodes, 'RATESRB', this%memoryPath)
      if (this%isrb == 1) then
        call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
      else
        call mem_allocate(this%sp2, nodes, 'SP2', this%memoryPath)
      end if
    end if
    !
    ! -- Initialize
    do n = 1, nodes
      this%porosity(n) = dzero
      this%thetam(n) = dzero
      this%volfracim(n) = dzero
      this%ratesto(n) = dzero
    end do
    do n = 1, size(this%decay)
      this%decay(n) = dzero
      this%ratedcy(n) = dzero
      this%decaylast(n) = dzero
    end do
    do n = 1, size(this%bulk_density)
      this%bulk_density(n) = dzero
      this%distcoef(n) = dzero
      this%ratesrb(n) = dzero
    end do
    do n = 1, size(this%sp2)
      this%sp2(n) = dzero
    end do
    do n = 1, size(this%ratedcys)
      this%ratedcys(n) = dzero
      this%decayslast(n) = dzero
    end do
    !
    ! -- Return
    return

allocate_scalars()

subroutine gwtmstmodule::allocate_scalars

(

class(gwtmsttype)

this

)

Method to allocate scalar variables for the package.

Parameters

this	GwtMstType object

Definition at line 1024 of file gwt-mst.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    ! -- local
    !
    ! -- Allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- Allocate
    call mem_allocate(this%isrb, 'ISRB', this%memoryPath)
    call mem_allocate(this%idcy, 'IDCY', this%memoryPath)
    !
    ! -- Initialize
    this%isrb = 0
    this%idcy = 0
    !
    ! -- Return
    return

get_freundlich_conc()

real(dp) function gwtmstmodule::get_freundlich_conc	(	real(dp), intent(in)	conc,
		real(dp), intent(in)	kf,
		real(dp), intent(in)	a
	)

Function to calculate sorption concentration using Freundlich

Parameters

[in]	conc	solute concentration
[in]	kf	freundlich constant
[in]	a	freundlich exponent

Definition at line 1466 of file gwt-mst.f90.

    ! -- dummy
    real(DP), intent(in) :: conc !< solute concentration
    real(DP), intent(in) :: kf !< freundlich constant
    real(DP), intent(in) :: a !< freundlich exponent
    ! -- return
    real(DP) :: cbar
    !
    if (conc > dzero) then
      cbar = kf * conc**a
    else
      cbar = dzero
    end if
    return

Here is the caller graph for this function:

get_freundlich_derivative()

real(dp) function gwtmstmodule::get_freundlich_derivative	(	real(dp), intent(in)	conc,
		real(dp), intent(in)	kf,
		real(dp), intent(in)	a
	)

Function to calculate sorption derivative using Freundlich

Parameters

[in]	conc	solute concentration
[in]	kf	freundlich constant
[in]	a	freundlich exponent

Definition at line 1508 of file gwt-mst.f90.

    ! -- dummy
    real(DP), intent(in) :: conc !< solute concentration
    real(DP), intent(in) :: kf !< freundlich constant
    real(DP), intent(in) :: a !< freundlich exponent
    ! -- return
    real(DP) :: derv
    !
    if (conc > dzero) then
      derv = kf * a * conc**(a - done)
    else
      derv = dzero
    end if
    return

Here is the caller graph for this function:

get_langmuir_conc()

real(dp) function gwtmstmodule::get_langmuir_conc	(	real(dp), intent(in)	conc,
		real(dp), intent(in)	kl,
		real(dp), intent(in)	sbar
	)

Function to calculate sorption concentration using Langmuir

Parameters

[in]	conc	solute concentration
[in]	kl	langmuir constant
[in]	sbar	langmuir sorption sites

Definition at line 1487 of file gwt-mst.f90.

    ! -- dummy
    real(DP), intent(in) :: conc !< solute concentration
    real(DP), intent(in) :: kl !< langmuir constant
    real(DP), intent(in) :: sbar !< langmuir sorption sites
    ! -- return
    real(DP) :: cbar
    !
    if (conc > dzero) then
      cbar = (kl * sbar * conc) / (done + kl * conc)
    else
      cbar = dzero
    end if
    return

Here is the caller graph for this function:

get_langmuir_derivative()

real(dp) function gwtmstmodule::get_langmuir_derivative	(	real(dp), intent(in)	conc,
		real(dp), intent(in)	kl,
		real(dp), intent(in)	sbar
	)

Function to calculate sorption derivative using Langmuir

Parameters

[in]	conc	solute concentration
[in]	kl	langmuir constant
[in]	sbar	langmuir sorption sites

Definition at line 1529 of file gwt-mst.f90.

    ! -- dummy
    real(DP), intent(in) :: conc !< solute concentration
    real(DP), intent(in) :: kl !< langmuir constant
    real(DP), intent(in) :: sbar !< langmuir sorption sites
    ! -- return
    real(DP) :: derv
    !
    if (conc > dzero) then
      derv = (kl * sbar) / (done + kl * conc)**dtwo
    else
      derv = dzero
    end if
    return

Here is the caller graph for this function:

get_volfracm()

real(dp) function gwtmstmodule::get_volfracm	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	node
	)

Calculate and return the volume fraction of the aquifer that is mobile

Parameters

	this	GwtMstType object
[in]	node	node number

Definition at line 1447 of file gwt-mst.f90.

    ! -- modules
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: node !< node number
    ! -- return
    real(DP) :: volfracm
    !
    volfracm = done - this%volfracim(node)
    !
    ! -- Return
    return

get_zero_order_decay()

real(dp) function gwtmstmodule::get_zero_order_decay	(	real(dp), intent(in)	decay_rate_usr,
		real(dp), intent(in)	decay_rate_last,
		integer(i4b), intent(in)	kiter,
		real(dp), intent(in)	cold,
		real(dp), intent(in)	cnew,
		real(dp), intent(in)	delt
	)

Function to calculate the zero-order decay rate from the user specified decay rate. If the decay rate is positive, then the decay rate must be constrained so that more mass is not removed than is available. Without this constraint, negative concentrations could result from zero-order decay.

Parameters

[in]	decay_rate_usr	user-entered decay rate
[in]	decay_rate_last	decay rate used for last iteration
[in]	kiter	Picard iteration counter
[in]	cold	concentration at end of last time step
[in]	cnew	concentration at end of this time step
[in]	delt	length of time step

Returns

returned value for decay rate

Definition at line 1553 of file gwt-mst.f90.

    ! -- dummy
    real(DP), intent(in) :: decay_rate_usr !< user-entered decay rate
    real(DP), intent(in) :: decay_rate_last !< decay rate used for last iteration
    integer(I4B), intent(in) :: kiter !< Picard iteration counter
    real(DP), intent(in) :: cold !< concentration at end of last time step
    real(DP), intent(in) :: cnew !< concentration at end of this time step
    real(DP), intent(in) :: delt !< length of time step
    ! -- return
    real(DP) :: decay_rate !< returned value for decay rate
    !
    ! -- Return user rate if production, otherwise constrain, if necessary
    if (decay_rate_usr < dzero) then
      !
      ! -- Production, no need to limit rate
      decay_rate = decay_rate_usr
    else
      !
      ! -- Need to ensure decay does not result in negative
      !    concentration, so reduce the rate if it would result in
      !    removing more mass than is in the cell.
      if (kiter == 1) then
        decay_rate = min(decay_rate_usr, cold / delt)
      else
        decay_rate = decay_rate_last
        if (cnew < dzero) then
          decay_rate = decay_rate_last + cnew / delt
        else if (cnew > cold) then
          decay_rate = decay_rate_last + cnew / delt
        end if
        decay_rate = min(decay_rate_usr, decay_rate)
      end if
      decay_rate = max(decay_rate, dzero)
    end if
    return

Here is the caller graph for this function:

mst_ar()

subroutine gwtmstmodule::mst_ar	(	class(gwtmsttype), intent(inout)	this,
		class(disbasetype), intent(in), pointer	dis,
		integer(i4b), dimension(:), pointer, contiguous	ibound
	)

Method to allocate and read static data for the package.

Parameters

[in,out]	this	GwtMstType object
[in]	dis	pointer to dis package
	ibound	pointer to GWT ibound array

Definition at line 131 of file gwt-mst.f90.

    ! -- modules
    ! -- dummy
    class(GwtMstType), intent(inout) :: this !< GwtMstType object
    class(DisBaseType), pointer, intent(in) :: dis !< pointer to dis package
    integer(I4B), dimension(:), pointer, contiguous :: ibound !< pointer to GWT ibound array
    ! -- local
    ! -- formats
    character(len=*), parameter :: fmtmst = &
      "(1x,/1x,'MST -- MOBILE STORAGE AND TRANSFER PACKAGE, VERSION 1, &
      &7/29/2020 INPUT READ FROM UNIT ', i0, //)"
    !
    ! --print a message identifying the mobile storage and transfer package.
    write (this%iout, fmtmst) this%inunit
    !
    ! -- Read options
    call this%read_options()
    !
    ! -- store pointers to arguments that were passed in
    this%dis => dis
    this%ibound => ibound
    !
    ! -- Allocate arrays
    call this%allocate_arrays(dis%nodes)
    !
    ! -- read the data block
    call this%read_data()
    !
    ! -- Return
    return

mst_bd()

subroutine gwtmstmodule::mst_bd	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	isuppress_output,
		type(budgettype), intent(inout)	model_budget
	)

Method to calculate budget terms for the package.

Parameters

	this	GwtMstType object
[in]	isuppress_output	flag to suppress output
[in,out]	model_budget	model budget object

Definition at line 875 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    use budgetmodule, only: budgettype, rate_accumulator
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: isuppress_output !< flag to suppress output
    type(BudgetType), intent(inout) :: model_budget !< model budget object
    ! -- local
    real(DP) :: rin
    real(DP) :: rout
    !
    ! -- sto
    call rate_accumulator(this%ratesto, rin, rout)
    call model_budget%addentry(rin, rout, delt, budtxt(1), &
                               isuppress_output, rowlabel=this%packName)
    !
    ! -- dcy
    if (this%idcy /= 0) then
      call rate_accumulator(this%ratedcy, rin, rout)
      call model_budget%addentry(rin, rout, delt, budtxt(2), &
                                 isuppress_output, rowlabel=this%packName)
    end if
    !
    ! -- srb
    if (this%isrb /= 0) then
      call rate_accumulator(this%ratesrb, rin, rout)
      call model_budget%addentry(rin, rout, delt, budtxt(3), &
                                 isuppress_output, rowlabel=this%packName)
    end if
    !
    ! -- srb dcy
    if (this%isrb /= 0 .and. this%idcy /= 0) then
      call rate_accumulator(this%ratedcys, rin, rout)
      call model_budget%addentry(rin, rout, delt, budtxt(4), &
                                 isuppress_output, rowlabel=this%packName)
    end if
    !
    ! -- Return
    return

Here is the call graph for this function:

mst_cq()

subroutine gwtmstmodule::mst_cq	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cnew,
		real(dp), dimension(nodes), intent(in)	cold,
		real(dp), dimension(:), intent(inout), contiguous	flowja
	)

Method to calculate flows for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cnew	concentration at end of this time step
[in]	cold	concentration at end of last time step
[in,out]	flowja	flow between two connected control volumes

Definition at line 565 of file gwt-mst.f90.

    ! -- modules
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
    ! -- local
    !
    ! - storage
    call this%mst_cq_sto(nodes, cnew, cold, flowja)
    !
    ! -- decay
    if (this%idcy /= 0) then
      call this%mst_cq_dcy(nodes, cnew, cold, flowja)
    end if
    !
    ! -- sorption
    if (this%isrb /= 0) then
      call this%mst_cq_srb(nodes, cnew, cold, flowja)
    end if
    !
    ! -- decay sorbed
    if (this%isrb /= 0 .and. this%idcy /= 0) then
      call this%mst_cq_dcy_srb(nodes, cnew, cold, flowja)
    end if
    !
    ! -- Return
    return

mst_cq_dcy()

subroutine gwtmstmodule::mst_cq_dcy	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cnew,
		real(dp), dimension(nodes), intent(in)	cold,
		real(dp), dimension(:), intent(inout), contiguous	flowja
	)

Method to calculate decay terms for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cnew	concentration at end of this time step
[in]	cold	concentration at end of last time step
[in,out]	flowja	flow between two connected control volumes

Definition at line 654 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: rate
    real(DP) :: swtpdt
    real(DP) :: hhcof, rrhs
    real(DP) :: vcell
    real(DP) :: decay_rate
    !
    ! -- initialize
    !
    ! -- Calculate decay change
    do n = 1, nodes
      !
      ! -- skip if transport inactive
      this%ratedcy(n) = dzero
      if (this%ibound(n) <= 0) cycle
      !
      ! -- calculate new and old water volumes
      vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
      swtpdt = this%fmi%gwfsat(n)
      !
      ! -- calculate decay gains and losses
      rate = dzero
      hhcof = dzero
      rrhs = dzero
      if (this%idcy == 1) then
        hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
      elseif (this%idcy == 2) then
        decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), &
                                          0, cold(n), cnew(n), delt)
        rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
      end if
      rate = hhcof * cnew(n) - rrhs
      this%ratedcy(n) = rate
      idiag = this%dis%con%ia(n)
      flowja(idiag) = flowja(idiag) + rate
      !
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

mst_cq_dcy_srb()

subroutine gwtmstmodule::mst_cq_dcy_srb	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cnew,
		real(dp), dimension(nodes), intent(in)	cold,
		real(dp), dimension(:), intent(inout), contiguous	flowja
	)

Method to calculate decay-sorption terms for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cnew	concentration at end of this time step
[in]	cold	concentration at end of last time step
[in,out]	flowja	flow between two connected control volumes

Definition at line 772 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: rate
    real(DP) :: hhcof, rrhs
    real(DP) :: vcell
    real(DP) :: swnew
    real(DP) :: distcoef
    real(DP) :: volfracm
    real(DP) :: rhobm
    real(DP) :: term
    real(DP) :: csrb
    real(DP) :: csrbnew
    real(DP) :: csrbold
    real(DP) :: decay_rate
    !
    ! -- Calculate sorbed decay change
    !    This routine will only be called if sorption and decay are active
    do n = 1, nodes
      !
      ! -- initialize rates
      this%ratedcys(n) = dzero
      !
      ! -- skip if transport inactive
      if (this%ibound(n) <= 0) cycle
      !
      ! -- set variables
      hhcof = dzero
      rrhs = dzero
      vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
      swnew = this%fmi%gwfsat(n)
      distcoef = this%distcoef(n)
      volfracm = this%get_volfracm(n)
      rhobm = this%bulk_density(n)
      term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
      !
      ! -- add sorbed mass decay rate terms to accumulators
      if (this%idcy == 1) then
        !
        if (this%isrb == 1) then
          !
          ! -- first order decay rate is a function of concentration, so add
          !    to left hand side
          hhcof = -term * distcoef
        else if (this%isrb == 2) then
          !
          ! -- nonlinear Freundlich sorption, so add to RHS
          csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
          rrhs = term * csrb
        else if (this%isrb == 3) then
          !
          ! -- nonlinear Lanmuir sorption, so add to RHS
          csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
          rrhs = term * csrb
        end if
      elseif (this%idcy == 2) then
        !
        ! -- Call function to get zero-order decay rate, which may be changed
        !    from the user-specified rate to prevent negative concentrations
        if (distcoef > dzero) then
          if (this%isrb == 1) then
            csrbold = cold(n) * distcoef
            csrbnew = cnew(n) * distcoef
          else if (this%isrb == 2) then
            csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n))
            csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
          else if (this%isrb == 3) then
            csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n))
            csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
          end if
          decay_rate = get_zero_order_decay(this%decay_sorbed(n), &
                                            this%decayslast(n), &
                                            0, csrbold, csrbnew, delt)
          rrhs = decay_rate * volfracm * rhobm * swnew * vcell
        end if
      end if
      !
      ! -- calculate rate
      rate = hhcof * cnew(n) - rrhs
      this%ratedcys(n) = rate
      idiag = this%dis%con%ia(n)
      flowja(idiag) = flowja(idiag) + rate
      !
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

mst_cq_srb()

subroutine gwtmstmodule::mst_cq_srb	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cnew,
		real(dp), dimension(nodes), intent(in)	cold,
		real(dp), dimension(:), intent(inout), contiguous	flowja
	)

Method to calculate sorption terms for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cnew	concentration at end of this time step
[in]	cold	concentration at end of last time step
[in,out]	flowja	flow between two connected control volumes

Definition at line 712 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: rate
    real(DP) :: tled
    real(DP) :: swt, swtpdt
    real(DP) :: vcell
    real(DP) :: volfracm
    real(DP) :: rhobm
    real(DP) :: const1
    real(DP) :: const2
    !
    ! -- initialize
    tled = done / delt
    !
    ! -- Calculate sorption change
    do n = 1, nodes
      !
      ! -- initialize rates
      this%ratesrb(n) = dzero
      !
      ! -- skip if transport inactive
      if (this%ibound(n) <= 0) cycle
      !
      ! -- assign variables
      vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
      swtpdt = this%fmi%gwfsat(n)
      swt = this%fmi%gwfsatold(n, delt)
      volfracm = this%get_volfracm(n)
      rhobm = this%bulk_density(n)
      const1 = this%distcoef(n)
      const2 = 0.
      if (this%isrb > 1) const2 = this%sp2(n)
      call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
                        cold(n), swtpdt, swt, const1, const2, &
                        rate=rate)
      this%ratesrb(n) = rate
      idiag = this%dis%con%ia(n)
      flowja(idiag) = flowja(idiag) + rate
      !
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

mst_cq_sto()

subroutine gwtmstmodule::mst_cq_sto	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cnew,
		real(dp), dimension(nodes), intent(in)	cold,
		real(dp), dimension(:), intent(inout), contiguous	flowja
	)

Method to calculate storage terms for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cnew	concentration at end of this time step
[in]	cold	concentration at end of last time step
[in,out]	flowja	flow between two connected control volumes

Definition at line 602 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: rate
    real(DP) :: tled
    real(DP) :: vnew, vold
    real(DP) :: hhcof, rrhs
    !
    ! -- initialize
    tled = done / delt
    !
    ! -- Calculate storage change
    do n = 1, nodes
      this%ratesto(n) = dzero
      !
      ! -- skip if transport inactive
      if (this%ibound(n) <= 0) cycle
      !
      ! -- calculate new and old water volumes
      vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
             this%fmi%gwfsat(n) * this%thetam(n)
      vold = vnew
      if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) * delt
      if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) * delt
      !
      ! -- calculate rate
      hhcof = -vnew * tled
      rrhs = -vold * tled * cold(n)
      rate = hhcof * cnew(n) - rrhs
      this%ratesto(n) = rate
      idiag = this%dis%con%ia(n)
      flowja(idiag) = flowja(idiag) + rate
    end do
    !
    ! -- Return
    return

mst_cr()

subroutine, public gwtmstmodule::mst_cr	(	type(gwtmsttype), pointer	mstobj,
		character(len=*), intent(in)	name_model,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout,
		type(tspfmitype), intent(in), target	fmi
	)

Create a new MST object

Parameters

	mstobj	unallocated new mst object to create
[in]	name_model	name of the model
[in]	inunit	unit number of WEL package input file
[in]	iout	unit number of model listing file
[in]	fmi	fmi package for this GWT model

Definition at line 97 of file gwt-mst.f90.

    ! -- dummy
    type(GwtMstType), pointer :: mstobj !< unallocated new mst object to create
    character(len=*), intent(in) :: name_model !< name of the model
    integer(I4B), intent(in) :: inunit !< unit number of WEL package input file
    integer(I4B), intent(in) :: iout !< unit number of model listing file
    type(TspFmiType), intent(in), target :: fmi !< fmi package for this GWT model
    !
    ! -- Create the object
    allocate (mstobj)
    !
    ! -- create name and memory path
    call mstobj%set_names(1, name_model, 'MST', 'MST')
    !
    ! -- Allocate scalars
    call mstobj%allocate_scalars()
    !
    ! -- Set variables
    mstobj%inunit = inunit
    mstobj%iout = iout
    mstobj%fmi => fmi
    !
    ! -- Initialize block parser
    call mstobj%parser%Initialize(mstobj%inunit, mstobj%iout)
    !
    ! -- Return
    return

Here is the caller graph for this function:

mst_da()

subroutine gwtmstmodule::mst_da

(

class(gwtmsttype)

this

)

Method to deallocate memory for the package.

Parameters

this	GwtMstType object

Definition at line 982 of file gwt-mst.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    !
    ! -- Deallocate arrays if package was active
    if (this%inunit > 0) then
      call mem_deallocate(this%porosity)
      call mem_deallocate(this%thetam)
      call mem_deallocate(this%volfracim)
      call mem_deallocate(this%ratesto)
      call mem_deallocate(this%idcy)
      call mem_deallocate(this%decay)
      call mem_deallocate(this%decay_sorbed)
      call mem_deallocate(this%ratedcy)
      call mem_deallocate(this%decaylast)
      call mem_deallocate(this%decayslast)
      call mem_deallocate(this%isrb)
      call mem_deallocate(this%bulk_density)
      call mem_deallocate(this%distcoef)
      call mem_deallocate(this%sp2)
      call mem_deallocate(this%ratesrb)
      call mem_deallocate(this%ratedcys)
      this%ibound => null()
      this%fmi => null()
    end if
    !
    ! -- Scalars
    !
    ! -- deallocate parent
    call this%NumericalPackageType%da()
    !
    ! -- Return
    return

mst_fc()

subroutine gwtmstmodule::mst_fc	(	class(gwtmsttype)	this,
		integer, intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cold,
		integer(i4b), intent(in)	nja,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(nja), intent(in)	idxglo,
		real(dp), dimension(nodes), intent(in)	cnew,
		real(dp), dimension(nodes), intent(inout)	rhs,
		integer(i4b), intent(in)	kiter
	)

Method to calculate and fill coefficients for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cold	concentration at end of last time step
[in]	nja	number of GWT connections
	matrix_sln	solution matrix
[in]	idxglo	mapping vector for model (local) to solution (global)
[in,out]	rhs	right-hand side vector for model
[in]	cnew	concentration at end of this time step
[in]	kiter	solution outer iteration number

Definition at line 168 of file gwt-mst.f90.

    ! -- modules
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer, intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    integer(I4B), intent(in) :: nja !< number of GWT connections
    class(MatrixBaseType), pointer :: matrix_sln !< solution matrix
    integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
    real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    integer(I4B), intent(in) :: kiter !< solution outer iteration number
    ! -- local
    !
    ! -- storage contribution
    call this%mst_fc_sto(nodes, cold, nja, matrix_sln, idxglo, rhs)
    !
    ! -- decay contribution
    if (this%idcy /= 0) then
      call this%mst_fc_dcy(nodes, cold, cnew, nja, matrix_sln, idxglo, &
                           rhs, kiter)
    end if
    !
    ! -- sorption contribution
    if (this%isrb /= 0) then
      call this%mst_fc_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
    end if
    !
    ! -- decay sorbed contribution
    if (this%isrb /= 0 .and. this%idcy /= 0) then
      call this%mst_fc_dcy_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, &
                               cnew, kiter)
    end if
    !
    ! -- Return
    return

mst_fc_dcy()

subroutine gwtmstmodule::mst_fc_dcy	(	class(gwtmsttype)	this,
		integer, intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cold,
		real(dp), dimension(nodes), intent(in)	cnew,
		integer(i4b), intent(in)	nja,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(nja), intent(in)	idxglo,
		real(dp), dimension(nodes), intent(inout)	rhs,
		integer(i4b), intent(in)	kiter
	)

Method to calculate and fill decay coefficients for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cold	concentration at end of last time step
[in]	cnew	concentration at end of this time step
[in]	nja	number of GWT connections
	matrix_sln	solution coefficient matrix
[in]	idxglo	mapping vector for model (local) to solution (global)
[in,out]	rhs	right-hand side vector for model
[in]	kiter	solution outer iteration number

Definition at line 262 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer, intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    integer(I4B), intent(in) :: nja !< number of GWT connections
    class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
    integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
    real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
    integer(I4B), intent(in) :: kiter !< solution outer iteration number
    ! -- local
    integer(I4B) :: n, idiag
    real(DP) :: hhcof, rrhs
    real(DP) :: swtpdt
    real(DP) :: vcell
    real(DP) :: decay_rate
    !
    ! -- loop through and calculate decay contribution to hcof and rhs
    do n = 1, this%dis%nodes
      !
      ! -- skip if transport inactive
      if (this%ibound(n) <= 0) cycle
      !
      ! -- calculate new and old water volumes
      vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
      swtpdt = this%fmi%gwfsat(n)
      !
      ! -- add decay rate terms to accumulators
      idiag = this%dis%con%ia(n)
      if (this%idcy == 1) then
        !
        ! -- first order decay rate is a function of concentration, so add
        !    to left hand side
        hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
        call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
      elseif (this%idcy == 2) then
        !
        ! -- Call function to get zero-order decay rate, which may be changed
        !    from the user-specified rate to prevent negative concentrations
        decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), &
                                          kiter, cold(n), cnew(n), delt)
        this%decaylast(n) = decay_rate
        rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
        rhs(n) = rhs(n) + rrhs
      end if
      !
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

mst_fc_dcy_srb()

subroutine gwtmstmodule::mst_fc_dcy_srb	(	class(gwtmsttype)	this,
		integer, intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cold,
		integer(i4b), intent(in)	nja,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(nja), intent(in)	idxglo,
		real(dp), dimension(nodes), intent(inout)	rhs,
		real(dp), dimension(nodes), intent(in)	cnew,
		integer(i4b), intent(in)	kiter
	)

Method to calculate and fill sorption-decay coefficients for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cold	concentration at end of last time step
[in]	nja	number of GWT connections
	matrix_sln	solution coefficient matrix
[in]	idxglo	mapping vector for model (local) to solution (global)
[in,out]	rhs	right-hand side vector for model
[in]	cnew	concentration at end of this time step
[in]	kiter	solution outer iteration number

Definition at line 460 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer, intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    integer(I4B), intent(in) :: nja !< number of GWT connections
    class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
    integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
    real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    integer(I4B), intent(in) :: kiter !< solution outer iteration number
    ! -- local
    integer(I4B) :: n, idiag
    real(DP) :: hhcof, rrhs
    real(DP) :: vcell
    real(DP) :: swnew
    real(DP) :: distcoef
    real(DP) :: volfracm
    real(DP) :: rhobm
    real(DP) :: term
    real(DP) :: csrb
    real(DP) :: decay_rate
    real(DP) :: csrbold
    real(DP) :: csrbnew
    !
    ! -- loop through and calculate sorption contribution to hcof and rhs
    do n = 1, this%dis%nodes
      !
      ! -- skip if transport inactive
      if (this%ibound(n) <= 0) cycle
      !
      ! -- set variables
      hhcof = dzero
      rrhs = dzero
      vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
      swnew = this%fmi%gwfsat(n)
      distcoef = this%distcoef(n)
      idiag = this%dis%con%ia(n)
      volfracm = this%get_volfracm(n)
      rhobm = this%bulk_density(n)
      term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
      !
      ! -- add sorbed mass decay rate terms to accumulators
      if (this%idcy == 1) then
        !
        if (this%isrb == 1) then
          !
          ! -- first order decay rate is a function of concentration, so add
          !    to left hand side
          hhcof = -term * distcoef
        else if (this%isrb == 2) then
          !
          ! -- nonlinear Freundlich sorption, so add to RHS
          csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
          rrhs = term * csrb
        else if (this%isrb == 3) then
          !
          ! -- nonlinear Lanmuir sorption, so add to RHS
          csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
          rrhs = term * csrb
        end if
      elseif (this%idcy == 2) then
        !
        ! -- Call function to get zero-order decay rate, which may be changed
        !    from the user-specified rate to prevent negative concentrations
        if (distcoef > dzero) then
 
          if (this%isrb == 1) then
            csrbold = cold(n) * distcoef
            csrbnew = cnew(n) * distcoef
          else if (this%isrb == 2) then
            csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n))
            csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
          else if (this%isrb == 3) then
            csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n))
            csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
          end if
 
          decay_rate = get_zero_order_decay(this%decay_sorbed(n), &
                                            this%decayslast(n), &
                                            kiter, csrbold, csrbnew, delt)
          this%decayslast(n) = decay_rate
          rrhs = decay_rate * volfracm * rhobm * swnew * vcell
        end if
 
      end if
      !
      ! -- Add hhcof to diagonal and rrhs to right-hand side
      call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
      rhs(n) = rhs(n) + rrhs
      !
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

mst_fc_srb()

subroutine gwtmstmodule::mst_fc_srb	(	class(gwtmsttype)	this,
		integer, intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cold,
		integer(i4b), intent(in)	nja,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(nja), intent(in)	idxglo,
		real(dp), dimension(nodes), intent(inout)	rhs,
		real(dp), dimension(nodes), intent(in)	cnew
	)

Method to calculate and fill sorption coefficients for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cold	concentration at end of last time step
[in]	nja	number of GWT connections
	matrix_sln	solution coefficient matrix
[in]	idxglo	mapping vector for model (local) to solution (global)
[in,out]	rhs	right-hand side vector for model
[in]	cnew	concentration at end of this time step

Definition at line 323 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer, intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    integer(I4B), intent(in) :: nja !< number of GWT connections
    class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
    integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
    real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
    real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
    ! -- local
    integer(I4B) :: n, idiag
    real(DP) :: tled
    real(DP) :: hhcof, rrhs
    real(DP) :: swt, swtpdt
    real(DP) :: vcell
    real(DP) :: const1
    real(DP) :: const2
    real(DP) :: volfracm
    real(DP) :: rhobm
    !
    ! -- set variables
    tled = done / delt
    !
    ! -- loop through and calculate sorption contribution to hcof and rhs
    do n = 1, this%dis%nodes
      !
      ! -- skip if transport inactive
      if (this%ibound(n) <= 0) cycle
      !
      ! -- assign variables
      vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
      swtpdt = this%fmi%gwfsat(n)
      swt = this%fmi%gwfsatold(n, delt)
      idiag = this%dis%con%ia(n)
      const1 = this%distcoef(n)
      const2 = 0.
      if (this%isrb > 1) const2 = this%sp2(n)
      volfracm = this%get_volfracm(n)
      rhobm = this%bulk_density(n)
      call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
                        cold(n), swtpdt, swt, const1, const2, &
                        hcofval=hhcof, rhsval=rrhs)
      !
      ! -- Add hhcof to diagonal and rrhs to right-hand side
      call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
      rhs(n) = rhs(n) + rrhs
      !
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

mst_fc_sto()

subroutine gwtmstmodule::mst_fc_sto	(	class(gwtmsttype)	this,
		integer, intent(in)	nodes,
		real(dp), dimension(nodes), intent(in)	cold,
		integer(i4b), intent(in)	nja,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(nja), intent(in)	idxglo,
		real(dp), dimension(nodes), intent(inout)	rhs
	)

Method to calculate and fill storage coefficients for the package.

Parameters

	this	GwtMstType object
[in]	nodes	number of nodes
[in]	cold	concentration at end of last time step
[in]	nja	number of GWT connections
	matrix_sln	solution coefficient matrix
[in]	idxglo	mapping vector for model (local) to solution (global)
[in,out]	rhs	right-hand side vector for model

Definition at line 212 of file gwt-mst.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer, intent(in) :: nodes !< number of nodes
    real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
    integer(I4B), intent(in) :: nja !< number of GWT connections
    class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
    integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
    real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
    ! -- local
    integer(I4B) :: n, idiag
    real(DP) :: tled
    real(DP) :: hhcof, rrhs
    real(DP) :: vnew, vold
    !
    ! -- set variables
    tled = done / delt
    !
    ! -- loop through and calculate storage contribution to hcof and rhs
    do n = 1, this%dis%nodes
      !
      ! -- skip if transport inactive
      if (this%ibound(n) <= 0) cycle
      !
      ! -- calculate new and old water volumes
      vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
             this%fmi%gwfsat(n) * this%thetam(n)
      vold = vnew
      if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) * delt
      if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) * delt
      !
      ! -- add terms to diagonal and rhs accumulators
      hhcof = -vnew * tled
      rrhs = -vold * tled * cold(n)
      idiag = this%dis%con%ia(n)
      call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
      rhs(n) = rhs(n) + rrhs
    end do
    !
    ! -- Return
    return

mst_ot_flow()

subroutine gwtmstmodule::mst_ot_flow	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	icbcfl,
		integer(i4b), intent(in)	icbcun
	)

Method to output terms for the package.

Parameters

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

Definition at line 922 of file gwt-mst.f90.

    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    integer(I4B), intent(in) :: icbcfl !< flag and unit number for cell-by-cell output
    integer(I4B), intent(in) :: icbcun !< flag indication if cell-by-cell data should be saved
    ! -- local
    integer(I4B) :: ibinun
    !character(len=16), dimension(2) :: aname
    integer(I4B) :: iprint, nvaluesp, nwidthp
    character(len=1) :: cdatafmp = ' ', editdesc = ' '
    real(DP) :: dinact
    !
    ! -- Set unit number for binary output
    if (this%ipakcb < 0) then
      ibinun = icbcun
    elseif (this%ipakcb == 0) then
      ibinun = 0
    else
      ibinun = this%ipakcb
    end if
    if (icbcfl == 0) ibinun = 0
    !
    ! -- Record the storage rate if requested
    if (ibinun /= 0) then
      iprint = 0
      dinact = dzero
      !
      ! -- sto
      call this%dis%record_array(this%ratesto, this%iout, iprint, -ibinun, &
                                 budtxt(1), cdatafmp, nvaluesp, &
                                 nwidthp, editdesc, dinact)
      !
      ! -- dcy
      if (this%idcy /= 0) &
        call this%dis%record_array(this%ratedcy, this%iout, iprint, -ibinun, &
                                   budtxt(2), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      !
      ! -- srb
      if (this%isrb /= 0) &
        call this%dis%record_array(this%ratesrb, this%iout, iprint, -ibinun, &
                                   budtxt(3), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      !
      ! -- dcy srb
      if (this%isrb /= 0 .and. this%idcy /= 0) &
        call this%dis%record_array(this%ratedcys, this%iout, iprint, -ibinun, &
                                   budtxt(4), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
    end if
    !
    ! -- Return
    return

mst_srb_term()

subroutine gwtmstmodule::mst_srb_term	(	integer(i4b), intent(in)	isrb,
		real(dp), intent(in)	volfracm,
		real(dp), intent(in)	rhobm,
		real(dp), intent(in)	vcell,
		real(dp), intent(in)	tled,
		real(dp), intent(in)	cnew,
		real(dp), intent(in)	cold,
		real(dp), intent(in)	swnew,
		real(dp), intent(in)	swold,
		real(dp), intent(in)	const1,
		real(dp), intent(in)	const2,
		real(dp), intent(out), optional	rate,
		real(dp), intent(out), optional	hcofval,
		real(dp), intent(out), optional	rhsval
	)

Subroutine to calculate sorption terms

Parameters

[in]	isrb	sorption flag 1, 2, 3 are linear, freundlich, and langmuir
[in]	volfracm	volume fraction of mobile domain (fhat_m)
[in]	rhobm	bulk density of mobile domain (rhob_m)
[in]	vcell	volume of cell
[in]	tled	one over time step length
[in]	cnew	concentration at end of this time step
[in]	cold	concentration at end of last time step
[in]	swnew	cell saturation at end of this time step
[in]	swold	cell saturation at end of last time step
[in]	const1	distribution coefficient or freundlich or langmuir constant
[in]	const2	zero, freundlich exponent, or langmuir sorption sites
[out]	rate	calculated sorption rate
[out]	hcofval	diagonal contribution to solution coefficient matrix
[out]	rhsval	contribution to solution right-hand-side

Definition at line 385 of file gwt-mst.f90.

    ! -- modules
    ! -- dummy
    integer(I4B), intent(in) :: isrb !< sorption flag 1, 2, 3 are linear, freundlich, and langmuir
    real(DP), intent(in) :: volfracm !< volume fraction of mobile domain (fhat_m)
    real(DP), intent(in) :: rhobm !< bulk density of mobile domain (rhob_m)
    real(DP), intent(in) :: vcell !< volume of cell
    real(DP), intent(in) :: tled !< one over time step length
    real(DP), intent(in) :: cnew !< concentration at end of this time step
    real(DP), intent(in) :: cold !< concentration at end of last time step
    real(DP), intent(in) :: swnew !< cell saturation at end of this time step
    real(DP), intent(in) :: swold !< cell saturation at end of last time step
    real(DP), intent(in) :: const1 !< distribution coefficient or freundlich or langmuir constant
    real(DP), intent(in) :: const2 !< zero, freundlich exponent, or langmuir sorption sites
    real(DP), intent(out), optional :: rate !< calculated sorption rate
    real(DP), intent(out), optional :: hcofval !< diagonal contribution to solution coefficient matrix
    real(DP), intent(out), optional :: rhsval !< contribution to solution right-hand-side
    ! -- local
    real(DP) :: term
    real(DP) :: derv
    real(DP) :: cbarnew
    real(DP) :: cbarold
    real(DP) :: cavg
    real(DP) :: cbaravg
    real(DP) :: swavg
    !
    ! -- Calculate based on type of sorption
    if (isrb == 1) then
      ! -- linear
      term = -volfracm * rhobm * vcell * tled * const1
      if (present(hcofval)) hcofval = term * swnew
      if (present(rhsval)) rhsval = term * swold * cold
      if (present(rate)) rate = term * swnew * cnew - term * swold * cold
    else
      !
      ! -- calculate average aqueous concentration
      cavg = dhalf * (cold + cnew)
      !
      ! -- set values based on isotherm
      if (isrb == 2) then
        ! -- freundlich
        cbarnew = get_freundlich_conc(cnew, const1, const2)
        cbarold = get_freundlich_conc(cold, const1, const2)
        derv = get_freundlich_derivative(cavg, const1, const2)
      else if (isrb == 3) then
        ! -- langmuir
        cbarnew = get_langmuir_conc(cnew, const1, const2)
        cbarold = get_langmuir_conc(cold, const1, const2)
        derv = get_langmuir_derivative(cavg, const1, const2)
      end if
      !
      ! -- calculate hcof, rhs, and rate for freundlich and langmuir
      term = -volfracm * rhobm * vcell * tled
      cbaravg = (cbarold + cbarnew) * dhalf
      swavg = (swnew + swold) * dhalf
      if (present(hcofval)) then
        hcofval = term * derv * swavg
      end if
      if (present(rhsval)) then
        rhsval = term * derv * swavg * cold - term * cbaravg * (swnew - swold)
      end if
      if (present(rate)) then
        rate = term * derv * swavg * (cnew - cold) &
               + term * cbaravg * (swnew - swold)
      end if
    end if
    return

Here is the call graph for this function:

Here is the caller graph for this function:

read_data()

subroutine gwtmstmodule::read_data

(

class(gwtmsttype)

this

)

Method to read data for the package.

Parameters

this	GwtMstType object

Definition at line 1218 of file gwt-mst.f90.

    ! -- modules
    use constantsmodule, only: linelength
    use memorymanagermodule, only: mem_reallocate, mem_reassignptr
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    ! -- local
    character(len=LINELENGTH) :: keyword
    character(len=:), allocatable :: line
    integer(I4B) :: istart, istop, lloc, ierr, n
    logical :: isfound, endOfBlock
    logical, dimension(6) :: lname
    character(len=24), dimension(6) :: aname
    ! -- formats
    ! -- data
    data aname(1)/'  MOBILE DOMAIN POROSITY'/
    data aname(2)/'            BULK DENSITY'/
    data aname(3)/'DISTRIBUTION COEFFICIENT'/
    data aname(4)/'              DECAY RATE'/
    data aname(5)/'       DECAY SORBED RATE'/
    data aname(6)/'   SECOND SORPTION PARAM'/
    !
    ! -- initialize
    isfound = .false.
    lname(:) = .false.
    !
    ! -- get griddata block
    call this%parser%GetBlock('GRIDDATA', isfound, ierr)
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING GRIDDATA'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        call this%parser%GetRemainingLine(line)
        lloc = 1
        select case (keyword)
        case ('POROSITY')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%porosity, &
                                        aname(1))
          lname(1) = .true.
        case ('BULK_DENSITY')
          if (this%isrb == 0) &
            call mem_reallocate(this%bulk_density, this%dis%nodes, &
                                'BULK_DENSITY', trim(this%memoryPath))
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, &
                                        this%bulk_density, aname(2))
          lname(2) = .true.
        case ('DISTCOEF')
          if (this%isrb == 0) &
            call mem_reallocate(this%distcoef, this%dis%nodes, 'DISTCOEF', &
                                trim(this%memoryPath))
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%distcoef, &
                                        aname(3))
          lname(3) = .true.
        case ('DECAY')
          if (this%idcy == 0) &
            call mem_reallocate(this%decay, this%dis%nodes, 'DECAY', &
                                trim(this%memoryPath))
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%decay, &
                                        aname(4))
          lname(4) = .true.
        case ('DECAY_SORBED')
          call mem_reallocate(this%decay_sorbed, this%dis%nodes, &
                              'DECAY_SORBED', trim(this%memoryPath))
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, &
                                        this%decay_sorbed, aname(5))
          lname(5) = .true.
        case ('SP2')
          if (this%isrb < 2) &
            call mem_reallocate(this%sp2, this%dis%nodes, 'SP2', &
                                trim(this%memoryPath))
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%sp2, &
                                        aname(6))
          lname(6) = .true.
        case default
          write (errmsg, '(a,a)') 'Unknown GRIDDATA tag: ', trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END PROCESSING GRIDDATA'
    else
      write (errmsg, '(a)') 'Required GRIDDATA block not found.'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- Check for required porosity
    if (.not. lname(1)) then
      write (errmsg, '(a)') 'POROSITY not specified in GRIDDATA block.'
      call store_error(errmsg)
    end if
    !
    ! -- Check for required sorption variables
    if (this%isrb > 0) then
      if (.not. lname(2)) then
        write (errmsg, '(a)') 'Sorption is active but BULK_DENSITY &
          &not specified.  BULK_DENSITY must be specified in GRIDDATA block.'
        call store_error(errmsg)
      end if
      if (.not. lname(3)) then
        write (errmsg, '(a)') 'Sorption is active but distribution &
          &coefficient not specified.  DISTCOEF must be specified in &
          &GRIDDATA block.'
        call store_error(errmsg)
      end if
      if (this%isrb > 1) then
        if (.not. lname(6)) then
          write (errmsg, '(a)') 'Freundlich or langmuir sorption is active &
            &but SP2 not specified.  SP2 must be specified in &
            &GRIDDATA block.'
          call store_error(errmsg)
        end if
      end if
    else
      if (lname(2)) then
        write (warnmsg, '(a)') 'Sorption is not active but &
          &BULK_DENSITY was specified.  BULK_DENSITY will have no affect on &
          &simulation results.'
        call store_warning(warnmsg)
        write (this%iout, '(1x,a)') 'WARNING.  '//warnmsg
      end if
      if (lname(3)) then
        write (warnmsg, '(a)') 'Sorption is not active but &
          &distribution coefficient was specified.  DISTCOEF will have &
          &no affect on simulation results.'
        call store_warning(warnmsg)
        write (this%iout, '(1x,a)') 'WARNING.  '//warnmsg
      end if
      if (lname(6)) then
        write (warnmsg, '(a)') 'Sorption is not active but &
          &SP2 was specified.  SP2 will have &
          &no affect on simulation results.'
        call store_warning(warnmsg)
        write (this%iout, '(1x,a)') 'WARNING.  '//warnmsg
      end if
    end if
    !
    ! -- Check for required decay/production rate coefficients
    if (this%idcy > 0) then
      if (.not. lname(4)) then
        write (errmsg, '(a)') 'First or zero order decay is &
          &active but the first rate coefficient is not specified.  DECAY &
          &must be specified in GRIDDATA block.'
        call store_error(errmsg)
      end if
      if (.not. lname(5)) then
        !
        ! -- If DECAY_SORBED not specified and sorption is active, then
        !    terminate with an error
        if (this%isrb > 0) then
          write (errmsg, '(a)') 'DECAY_SORBED not provided in GRIDDATA &
            &block but decay and sorption are active.  Specify DECAY_SORBED &
            &in GRIDDATA block.'
          call store_error(errmsg)
        end if
      end if
    else
      if (lname(4)) then
        write (warnmsg, '(a)') 'First- or zero-order decay &
          &is not active but decay was specified.  DECAY will &
          &have no affect on simulation results.'
        call store_warning(warnmsg)
        write (this%iout, '(1x,a)') 'WARNING.  '//warnmsg
      end if
      if (lname(5)) then
        write (warnmsg, '(a)') 'First- or zero-order decay &
          &is not active but DECAY_SORBED was specified.  &
          &DECAY_SORBED will have no affect on simulation results.'
        call store_warning(warnmsg)
        write (this%iout, '(1x,a)') 'WARNING.  '//warnmsg
      end if
    end if
    !
    ! -- terminate if errors
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- initialize thetam from porosity
    do n = 1, size(this%porosity)
      this%thetam(n) = this%porosity(n)
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

read_options()

subroutine gwtmstmodule::read_options

(

class(gwtmsttype)

this

)

Method to read options for the package.

Parameters

this	GwtMstType object

Definition at line 1141 of file gwt-mst.f90.

    ! -- modules
    use constantsmodule, only: linelength
    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    ! -- local
    character(len=LINELENGTH) :: keyword, keyword2
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    ! -- formats
    character(len=*), parameter :: fmtisvflow = &
      "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
      &WHENEVER ICBCFL IS NOT ZERO.')"
    character(len=*), parameter :: fmtisrb = &
      &"(4x,'LINEAR SORPTION IS ACTIVE. ')"
    character(len=*), parameter :: fmtfreundlich = &
      &"(4x,'FREUNDLICH SORPTION IS ACTIVE. ')"
    character(len=*), parameter :: fmtlangmuir = &
      &"(4x,'LANGMUIR SORPTION IS ACTIVE. ')"
    character(len=*), parameter :: fmtidcy1 = &
      &"(4x,'FIRST-ORDER DECAY IS ACTIVE. ')"
    character(len=*), parameter :: fmtidcy2 = &
      &"(4x,'ZERO-ORDER DECAY IS ACTIVE. ')"
    !
    ! -- get options block
    call this%parser%GetBlock('OPTIONS', isfound, ierr, blockrequired=.false., &
                              supportopenclose=.true.)
    !
    ! -- parse options block if detected
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING MOBILE STORAGE AND TRANSFER OPTIONS'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('SAVE_FLOWS')
          this%ipakcb = -1
          write (this%iout, fmtisvflow)
        case ('SORBTION', 'SORPTION')
          this%isrb = 1
          call this%parser%GetStringCaps(keyword2)
          if (trim(adjustl(keyword2)) == 'LINEAR') this%isrb = 1
          if (trim(adjustl(keyword2)) == 'FREUNDLICH') this%isrb = 2
          if (trim(adjustl(keyword2)) == 'LANGMUIR') this%isrb = 3
          select case (this%isrb)
          case (1)
            write (this%iout, fmtisrb)
          case (2)
            write (this%iout, fmtfreundlich)
          case (3)
            write (this%iout, fmtlangmuir)
          end select
        case ('FIRST_ORDER_DECAY')
          this%idcy = 1
          write (this%iout, fmtidcy1)
        case ('ZERO_ORDER_DECAY')
          this%idcy = 2
          write (this%iout, fmtidcy2)
        case default
          write (errmsg, '(a,a)') 'Unknown MST option: ', trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END OF MOBILE STORAGE AND TRANSFER OPTIONS'
    end if
    !
    ! -- Return
    return

Here is the call graph for this function:

Variable Documentation

budtxt

character(len=lenbudtxt), dimension(nbditems) gwtmstmodule::budtxt

Definition at line 27 of file gwt-mst.f90.

  character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt

nbditems

integer(i4b), parameter gwtmstmodule::nbditems = 4

Definition at line 26 of file gwt-mst.f90.

  integer(I4B), parameter :: NBDITEMS = 4