gwtmstmodule Module Reference

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

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

Enumerations
enum	{ decay_off = 0 , decay_first_order = 1 , decay_zero_order = 2 }
	Enumerator that defines the decay options. More...

Functions/Subroutines
subroutine, public	mst_cr (mstobj, name_model, input_mempath, 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_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_calc_csrb (this, cnew)
	@ brief Calculate sorbed concentration 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_ot_dv (this, idvsave)
	Save sorbate concentration array to binary file. 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	source_options (this)
	@ brief Source options for package More...
subroutine	log_options (this, found, sorbate_fname)
	Log user options to list file. More...
subroutine	source_data (this)
	@ brief Source data for package More...
subroutine	log_data (this, found)
	Log user data to list file. 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_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

Enumeration Type Documentation

anonymous enum

anonymous enum

Enumerator
decay_off	Decay (or production) of mass inactive (default)
decay_first_order	First-order decay.
decay_zero_order	Zeroth-order decay.

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

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 1408 of file gwt-mst.f90.

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

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 1015 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 == decay_off) 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 /= decay_off .and. this%isrb /= sorption_off) 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 == sorption_off) 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)
      call mem_allocate(this%csrb, 1, 'CSRB', 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)
      call mem_allocate(this%csrb, nodes, 'CSRB', this%memoryPath)
      if (this%isrb == sorption_linear) 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
      this%csrb(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

allocate_scalars()

subroutine gwtmstmodule::allocate_scalars

(

class(gwtmsttype)

this

)

Method to allocate scalar variables for the package.

Parameters

this	GwtMstType object

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

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

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 1431 of file gwt-mst.f90.

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

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 1449 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

Here is the caller graph for this function:

log_data()

subroutine gwtmstmodule::log_data	(	class(gwtmsttype)	this,
		type(gwtmstparamfoundtype), intent(in)	found
	)

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

    use gwtmstinputmodule, only: gwtmstparamfoundtype
    class(GwTMstType) :: this
    type(GwtMstParamFoundType), intent(in) :: found
 
    write (this%iout, '(1x,a)') 'PROCESSING GRIDDATA'
    if (found%porosity) then
      write (this%iout, '(4x,a)') 'MOBILE DOMAIN POROSITY set from input file'
    end if
    if (found%decay) then
      write (this%iout, '(4x,a)') 'DECAY RATE set from input file'
    end if
    if (found%decay_sorbed) then
      write (this%iout, '(4x,a)') 'DECAY SORBED RATE set from input file'
    end if
    if (found%bulk_density) then
      write (this%iout, '(4x,a)') 'BULK DENSITY set from input file'
    end if
    if (found%distcoef) then
      write (this%iout, '(4x,a)') 'DISTRIBUTION COEFFICIENT set from input file'
    end if
    if (found%sp2) then
      write (this%iout, '(4x,a)') 'SECOND SORPTION PARAM set from input file'
    end if
    write (this%iout, '(1x,a)') 'END PROCESSING GRIDDATA'

log_options()

subroutine gwtmstmodule::log_options	(	class(gwtmsttype)	this,
		type(gwtmstparamfoundtype), intent(in)	found,
		character(len=*), intent(in)	sorbate_fname
	)

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

    use gwtmstinputmodule, only: gwtmstparamfoundtype
    class(GwTMstType) :: this
    type(GwtMstParamFoundType), intent(in) :: found
    character(len=*), intent(in) :: sorbate_fname
    ! -- 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 :: fmtlinear = &
      &"(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. ')"
    character(len=*), parameter :: fmtfileout = &
      "(4x,'MST ',1x,a,1x,' WILL BE SAVED TO FILE: ',a,/4x,&
      &'OPENED ON UNIT: ',I7)"
 
    write (this%iout, '(1x,a)') 'PROCESSING MOBILE STORAGE AND TRANSFER OPTIONS'
    if (found%save_flows) then
      write (this%iout, fmtisvflow)
    end if
    if (found%order1_decay) then
      write (this%iout, fmtidcy1)
    end if
    if (found%order0_decay) then
      write (this%iout, fmtidcy2)
    end if
    if (found%sorption) then
      select case (this%isrb)
      case (sorption_linear)
        write (this%iout, fmtlinear)
      case (sorption_freund)
        write (this%iout, fmtfreundlich)
      case (sorption_lang)
        write (this%iout, fmtlangmuir)
      end select
    end if
    if (found%sorbatefile) then
      write (this%iout, fmtfileout) &
        'SORBATE', sorbate_fname, this%ioutsorbate
    end if
    write (this%iout, '(1x,a)') 'END OF MOBILE STORAGE AND TRANSFER OPTIONS'

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 142 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 MEMPATH: ', A, //)"
    !
    ! --print a message identifying the mobile storage and transfer package.
    write (this%iout, fmtmst) this%input_mempath
    !
    ! -- Source options
    call this%source_options()
    !
    ! -- store pointers to arguments that were passed in
    this%dis => dis
    this%ibound => ibound
    !
    ! -- Allocate arrays
    call this%allocate_arrays(dis%nodes)
    !
    ! -- source the data block
    call this%source_data()
    !
    ! -- Create isotherm object if sorption is active
    this%isotherm => create_isotherm(this%isrb, this%distcoef, this%sp2)

Here is the call graph for this function:

mst_bd()

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

Parameters

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

Definition at line 811 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 /= decay_off) 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 /= sorption_off) 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 /= sorption_off .and. this%idcy /= decay_off) then
      call rate_accumulator(this%ratedcys, rin, rout)
      call model_budget%addentry(rin, rout, delt, budtxt(4), &
                                 isuppress_output, rowlabel=this%packName)
    end if

Here is the call graph for this function:

mst_calc_csrb()

subroutine gwtmstmodule::mst_calc_csrb	(	class(gwtmsttype)	this,
		real(dp), dimension(:), intent(in)	cnew
	)

Parameters

	this	GwtMstType object
[in]	cnew	concentration at end of this time step

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

    ! -- dummy
    class(GwtMstType) :: this !< GwtMstType object
    real(DP), intent(in), dimension(:) :: cnew !< concentration at end of this time step
    ! -- local
    integer(I4B) :: n
    real(DP) :: csrb
 
    ! Calculate sorbed concentration
    do n = 1, size(cnew)
      csrb = dzero
      if (this%ibound(n) > 0 .and. this%isrb /= sorption_off) then
        csrb = this%isotherm%value(cnew, n)
      end if
      this%csrb(n) = csrb
    end do
 

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 485 of file gwt-mst.f90.

    ! -- 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
    !
    ! - storage
    call this%mst_cq_sto(nodes, cnew, cold, flowja)
    !
    ! -- decay
    if (this%idcy /= decay_off) then
      call this%mst_cq_dcy(nodes, cnew, cold, flowja)
    end if
    !
    ! -- sorption
    if (this%isrb /= sorption_off) then
      call this%mst_cq_srb(nodes, cnew, cold, flowja)
    end if
    !
    ! -- decay sorbed
    if (this%isrb /= sorption_off .and. this%idcy /= decay_off) then
      call this%mst_cq_dcy_srb(nodes, cnew, cold, flowja)
    end if
    !
    ! -- calculate csrb
    if (this%isrb /= sorption_off) then
      call this%mst_calc_csrb(cnew)
    end if

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 569 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 == decay_first_order) then
        if (cnew(n) > dzero) then
          hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
        end if
      elseif (this%idcy == decay_zero_order) 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

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 695 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
      select case (this%idcy)
      case (decay_first_order)
 
        !
        select case (this%isrb)
        case (sorption_linear)
          !
          ! -- first order decay rate is a function of concentration, so add
          !    to left hand side
          if (cnew(n) > dzero) then
            hhcof = -term * distcoef
          end if
        case (sorption_freund)
          !
          ! -- nonlinear Freundlich sorption, so add to RHS
          csrb = this%isotherm%value(cnew, n)
          rrhs = term * csrb
        case (sorption_lang)
          !
          ! -- nonlinear Lanmuir sorption, so add to RHS
          csrb = this%isotherm%value(cnew, n)
          rrhs = term * csrb
        end select
      case (decay_zero_order)
        !
        ! -- 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
          csrbold = this%isotherm%value(cold, n)
          csrbnew = this%isotherm%value(cnew, n)
 
          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 select
      !
      ! -- calculate rate
      rate = hhcof * cnew(n) - rrhs
      this%ratedcys(n) = rate
      idiag = this%dis%con%ia(n)
      flowja(idiag) = flowja(idiag) + rate
      !
    end do

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 625 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) :: vcell
    real(DP) :: volfracm
    real(DP) :: rhobm
    real(DP), dimension(nodes) :: c_half
    real(DP) :: cbar_derv_half
    real(DP) :: cbar_new, cbar_half, cbar_old
    real(DP) :: sat_new, sat_half, sat_old
    !
    ! -- initialize
    tled = done / delt
    c_half = 0.5_dp * (cold + cnew)
    !
    ! -- Calculate sorption change
    do n = 1, nodes
      !
      ! -- initialize rates
      this%ratesrb(n) = dzero
      rate = 0.0_dp
      !
      ! -- 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))
 
      rhobm = this%bulk_density(n)
      volfracm = this%get_volfracm(n)
      sat_new = this%fmi%gwfsat(n)
      sat_old = this%fmi%gwfsatold(n, delt)
 
      ! -- Midpoint formulation using average values
      cbar_new = this%isotherm%value(cnew, n)
      cbar_old = this%isotherm%value(cold, n)
      cbar_half = 0.5 * (cbar_new + cbar_old)
 
      sat_half = 0.5 * (sat_new + sat_old)
      cbar_derv_half = this%isotherm%derivative(c_half, n)
 
      rate = -volfracm * rhobm * cbar_derv_half * sat_half * cnew(n) &
             * vcell * tled
      rate = rate + volfracm * rhobm * cbar_derv_half * sat_half * cold(n) &
             * vcell * tled
      rate = rate - volfracm * rhobm * cbar_half * (sat_new - sat_old) &
             * vcell * tled
 
      this%ratesrb(n) = rate
      idiag = this%dis%con%ia(n)
      flowja(idiag) = flowja(idiag) + rate
      !
    end do

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 521 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

mst_cr()

subroutine, public gwtmstmodule::mst_cr	(	type(gwtmsttype), pointer	mstobj,
		character(len=*), intent(in)	name_model,
		character(len=*), intent(in)	input_mempath,
		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]	input_mempath	memory path of input
[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 114 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
    character(len=*), intent(in) :: input_mempath !< memory path of input
    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', input_mempath)
    !
    ! -- Allocate scalars
    call mstobj%allocate_scalars()
    !
    ! -- Set variables
    mstobj%inunit = inunit
    mstobj%iout = iout
    mstobj%fmi => fmi

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 945 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%ioutsorbate)
      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%csrb)
      call mem_deallocate(this%ratedcys)
      this%ibound => null()
      this%fmi => null()
    end if
    !
    ! -- Scalars
    !
    ! -- Objects
    if (associated(this%isotherm)) then
      deallocate (this%isotherm)
      nullify (this%isotherm)
    end if
    !
    ! -- deallocate parent
    call this%NumericalPackageType%da()

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 178 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
    !
    ! -- storage contribution
    call this%mst_fc_sto(nodes, cold, nja, matrix_sln, idxglo, rhs)
    !
    ! -- decay contribution
    if (this%idcy /= decay_off) then
      call this%mst_fc_dcy(nodes, cold, cnew, nja, matrix_sln, idxglo, &
                           rhs, kiter)
    end if
    !
    ! -- sorption contribution
    if (this%isrb /= sorption_off) then
      call this%mst_fc_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
    end if
    !
    ! -- decay sorbed contribution
    if (this%isrb /= sorption_off .and. this%idcy /= decay_off) then
      call this%mst_fc_dcy_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, &
                               cnew, kiter)
    end if

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 263 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)
      select case (this%idcy)
      case (decay_first_order)
        !
        ! -- first order decay rate is a function of concentration, so add
        !    to left hand side
        if (cnew(n) > dzero) then
          hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
          call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
        end if
      case (decay_zero_order)
        !
        ! -- 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 select
      !
    end do

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 391 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
      select case (this%idcy)
      case (decay_first_order)
        select case (this%isrb)
        case (sorption_linear)
          !
          ! -- first order decay rate is a function of concentration, so add
          !    to left hand side
          if (cnew(n) > dzero) then
            hhcof = -term * distcoef
          end if
        case (sorption_freund)
          !
          ! -- nonlinear Freundlich sorption, so add to RHS
          csrb = this%isotherm%value(cnew, n)
          rrhs = term * csrb
        case (sorption_lang)
          !
          ! -- nonlinear Lanmuir sorption, so add to RHS
          csrb = this%isotherm%value(cnew, n)
          rrhs = term * csrb
        end select
      case (decay_zero_order)
        !
        ! -- 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
          csrbold = this%isotherm%value(cold, n)
          csrbnew = this%isotherm%value(cnew, n)
          !
          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 select
      !
      ! -- 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

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) :: vcell
    real(DP) :: volfracm
    real(DP) :: rhobm
    real(DP), dimension(nodes) :: c_half
    real(DP) :: cbar_derv_half
    real(DP) :: cbar_new, cbar_half, cbar_old
    real(DP) :: sat_new, sat_half, sat_old
    !
    ! -- set variables
    tled = done / delt
    c_half = 0.5_dp * (cold + cnew)
    !
    ! -- 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))
      volfracm = this%get_volfracm(n)
      rhobm = this%bulk_density(n)
      sat_new = this%fmi%gwfsat(n)
      sat_old = this%fmi%gwfsatold(n, delt)
 
      ! -- Midpoint formulation using average values
      cbar_new = this%isotherm%value(cnew, n)
      cbar_old = this%isotherm%value(cold, n)
      cbar_half = 0.5 * (cbar_new + cbar_old)
 
      sat_half = 0.5 * (sat_new + sat_old)
      cbar_derv_half = this%isotherm%derivative(c_half, n)
 
      hhcof = -volfracm * rhobm * cbar_derv_half * sat_half * vcell * tled
      idiag = this%dis%con%ia(n)
      call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
 
      rrhs = -volfracm * rhobm * cbar_derv_half * sat_half * cold(n) &
             * vcell * tled
      rhs(n) = rhs(n) + rrhs
 
      rrhs = volfracm * rhobm * cbar_half * (sat_new - sat_old) * vcell * tled
      rhs(n) = rhs(n) + rrhs
 
    end do

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 217 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

mst_ot_dv()

subroutine gwtmstmodule::mst_ot_dv	(	class(gwtmsttype)	this,
		integer(i4b), intent(in)	idvsave
	)

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

    ! -- dummy
    class(GwtMstType) :: this
    integer(I4B), intent(in) :: idvsave
    ! -- 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 sorbate output
    if (this%ioutsorbate /= 0) then
      ibinun = 1
    else
      ibinun = 0
    end if
    if (idvsave == 0) ibinun = 0
 
    ! save sorbate concentration array
    if (ibinun /= 0) then
      iprint = 0
      dinact = dhnoflo
      if (this%ioutsorbate /= 0) then
        ibinun = this%ioutsorbate
        call this%dis%record_array(this%csrb, this%iout, iprint, ibinun, &
                                   '         SORBATE', cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      end if
    end if
 

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 854 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
    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 /= decay_off) &
        call this%dis%record_array(this%ratedcy, this%iout, iprint, -ibinun, &
                                   budtxt(2), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      !
      ! -- srb
      if (this%isrb /= sorption_off) &
        call this%dis%record_array(this%ratesrb, this%iout, iprint, -ibinun, &
                                   budtxt(3), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      !
      ! -- dcy srb
      if (this%isrb /= sorption_off .and. this%idcy /= decay_off) &
        call this%dis%record_array(this%ratedcys, this%iout, iprint, -ibinun, &
                                   budtxt(4), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
    end if

source_data()

subroutine gwtmstmodule::source_data

(

class(gwtmsttype)

this

)

Method to source data for the package.

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

    ! -- modules
    use memorymanagermodule, only: get_isize, mem_reallocate
    use memorymanagerextmodule, only: mem_set_value
    use gwtmstinputmodule, only: gwtmstparamfoundtype
    ! -- dummy
    class(GwtMsttype) :: this
    ! -- locals
    character(len=LINELENGTH) :: errmsg
    type(GwtMstParamFoundType) :: found
    integer(I4B) :: n, asize
    integer(I4B), dimension(:), pointer, contiguous :: map
    !
    ! -- set map to convert user input data into reduced data
    map => null()
    if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
    !
    ! -- reallocate
    if (this%isrb == sorption_off) then
      call get_isize('BULK_DENSITY', this%input_mempath, asize)
      if (asize > 0) &
        call mem_reallocate(this%bulk_density, this%dis%nodes, &
                            'BULK_DENSITY', this%memoryPath)
      call get_isize('DISTCOEF', this%input_mempath, asize)
      if (asize > 0) &
        call mem_reallocate(this%distcoef, this%dis%nodes, 'DISTCOEF', &
                            this%memoryPath)
    end if
    if (this%idcy == decay_off) then
      call get_isize('DECAY', this%input_mempath, asize)
      if (asize > 0) &
        call mem_reallocate(this%decay, this%dis%nodes, 'DECAY', this%memoryPath)
    end if
    call get_isize('DECAY_SORBED', this%input_mempath, asize)
    if (asize > 0) then
      call mem_reallocate(this%decay_sorbed, this%dis%nodes, &
                          'DECAY_SORBED', this%memoryPath)
    end if
    if (this%isrb == sorption_off .or. this%isrb == sorption_linear) then
      call get_isize('SP2', this%input_mempath, asize)
      if (asize > 0) &
        call mem_reallocate(this%sp2, this%dis%nodes, 'SP2', this%memoryPath)
    end if
    !
    ! -- update defaults with memory sourced values
    call mem_set_value(this%porosity, 'POROSITY', this%input_mempath, map, &
                       found%porosity)
    call mem_set_value(this%decay, 'DECAY', this%input_mempath, map, &
                       found%decay)
    call mem_set_value(this%decay_sorbed, 'DECAY_SORBED', this%input_mempath, &
                       map, found%decay_sorbed)
    call mem_set_value(this%bulk_density, 'BULK_DENSITY', this%input_mempath, &
                       map, found%bulk_density)
    call mem_set_value(this%distcoef, 'DISTCOEF', this%input_mempath, map, &
                       found%distcoef)
    call mem_set_value(this%sp2, 'SP2', this%input_mempath, map, &
                       found%sp2)
 
    ! -- log options
    if (this%iout > 0) then
      call this%log_data(found)
    end if
 
    ! -- Check for required porosity
    if (.not. found%porosity) then
      write (errmsg, '(a)') 'POROSITY not specified in GRIDDATA block.'
      call store_error(errmsg)
    end if
 
    ! -- Check for required sorption variables
    if (this%isrb == sorption_linear .or. this%isrb == sorption_freund .or. &
        this%isrb == sorption_lang) then
      if (.not. found%bulk_density) 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. found%distcoef) 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 == sorption_freund .or. this%isrb == sorption_lang) then
        if (.not. found%sp2) 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 (found%bulk_density) 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 (found%distcoef) 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 (found%sp2) 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 /= decay_off) then
      if (.not. found%decay) 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. found%decay_sorbed) then
        !
        ! -- If DECAY_SORBED not specified and sorption is active, then
        !    terminate with an error
        if (this%isrb == sorption_linear .or. this%isrb == sorption_freund .or. &
            this%isrb == sorption_lang) 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 (found%decay) 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 (found%decay_sorbed) 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 store_error_filename(this%input_fname)
    end if
 
    ! -- initialize thetam from porosity
    do n = 1, size(this%porosity)
      this%thetam(n) = this%porosity(n)
    end do

Here is the call graph for this function:

source_options()

subroutine gwtmstmodule::source_options

(

class(gwtmsttype)

this

)

Method to source options for the package.

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

    ! -- modules
    use constantsmodule, only: lenvarname
    use openspecmodule, only: access, form
    use inputoutputmodule, only: getunit, openfile
    use memorymanagerextmodule, only: mem_set_value
    use gwtmstinputmodule, only: gwtmstparamfoundtype
    ! -- dummy
    class(GwtMstType) :: this
    ! -- locals
    type(GwtMstParamFoundType) :: found
    character(len=LENVARNAME), dimension(3) :: sorption_method = &
      &[character(len=LENVARNAME) :: 'LINEAR', 'FREUNDLICH', 'LANGMUIR']
    character(len=linelength) :: fname
    !
    ! -- update defaults with memory sourced values
    call mem_set_value(this%ipakcb, 'SAVE_FLOWS', this%input_mempath, &
                       found%save_flows)
    call mem_set_value(this%idcy, 'ORDER1_DECAY', this%input_mempath, &
                       found%order1_decay)
    call mem_set_value(this%idcy, 'ORDER0_DECAY', this%input_mempath, &
                       found%order0_decay)
    call mem_set_value(this%isrb, 'SORPTION', this%input_mempath, &
                       sorption_method, found%sorption)
    call mem_set_value(fname, 'SORBATEFILE', this%input_mempath, &
                       found%sorbatefile)
 
    ! -- found side effects
    if (found%save_flows) this%ipakcb = -1
    if (found%order1_decay) this%idcy = decay_first_order
    if (found%order0_decay) this%idcy = decay_zero_order
    if (found%sorption) then
      if (this%isrb == izero) then
        call store_error('Unknown sorption type was specified. &
                         &Sorption must be specified as LINEAR, &
                         &FREUNDLICH, or LANGMUIR.')
        call store_error_filename(this%input_fname)
      end if
    end if
    if (found%sorbatefile) then
      this%ioutsorbate = getunit()
      call openfile(this%ioutsorbate, this%iout, fname, 'DATA(BINARY)', &
                    form, access, 'REPLACE', mode_opt=mnormal)
    end if
    !
    ! -- log options
    if (this%iout > 0) then
      call this%log_options(found, fname)
    end if

Here is the call graph for this function:

Variable Documentation

budtxt

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

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

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

nbditems

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

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

  integer(I4B), parameter :: NBDITEMS = 4