gwtistmodule Module Reference

– @ brief Immobile Storage and Transfer (IST) Module More...

Data Types
type	gwtisttype
	@ brief Immobile storage and transfer More...

Functions/Subroutines
subroutine, public	ist_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname, fmi, mst)
	@ brief Create a new package object More...
subroutine	ist_ar (this)
	@ brief Allocate and read method for package More...
subroutine	ist_rp (this)
	@ brief Read and prepare method for package More...
subroutine	ist_ad (this)
	@ brief Advance the ist package More...
subroutine	ist_fc (this, rhs, ia, idxglo, matrix_sln)
	@ brief Fill coefficient method for package More...
subroutine	ist_cq (this, x, flowja, iadv)
	@ brief Calculate package flows. More...
subroutine	ist_calc_csrb (this, cim)
	@ brief Calculate immobile sorbed concentration More...
subroutine	ist_bd (this, model_budget)
	@ brief Add package flows to model budget. More...
subroutine	ist_ot_model_flows (this, icbcfl, ibudfl, icbcun, imap)
	@ brief Output model flow terms. More...
subroutine	ist_ot_dv (this, idvsave, idvprint)
	@ brief Output dependent variables. More...
subroutine	output_immobile_concentration (this, idvsave, idvprint)
	@ brief Output immobile domain aqueous concentration. More...
subroutine	output_immobile_sorbate_concentration (this, idvsave, idvprint)
	@ brief Output immobile domain sorbate concentration. More...
subroutine	ist_ot_bdsummary (this, kstp, kper, iout, ibudfl)
	@ brief Output IST package budget summary. More...
subroutine	ist_da (this)
	@ brief Deallocate package memory More...
subroutine	allocate_scalars (this)
	@ brief Allocate package scalars More...
subroutine	ist_allocate_arrays (this)
	@ brief Allocate package arrays More...
subroutine	read_options (this)
	@ brief Read options for package More...
subroutine	ist_read_dimensions (this)
	@ brief Read dimensions for package More...
subroutine	read_data (this)
	@ brief Read data for package More...
real(dp) function	get_thetaim (this, node)
	@ brief Return thetaim More...
subroutine	get_ddterm (thetaim, vcell, delt, swtpdt, volfracim, rhobim, kdnew, kdold, lambda1im, lambda2im, gamma1im, gamma2im, zetaim, ddterm, f)
	@ brief Calculate immobile domain equation terms More...
subroutine	get_hcofrhs (ddterm, f, cimold, hcof, rhs)
	@ brief Calculate the hcof and rhs terms for immobile domain More...
real(dp) function	get_ddconc (ddterm, f, cimold, cnew)
	@ brief Calculate the immobile domain concentration More...
subroutine	accumulate_budterm (budterm, ddterm, cimnew, cimold, cnew, idcy)
	@ brief Calculate the immobile domain budget terms More...

Variables
character(len=lenftype)	ftype = 'IST'
character(len=lenpackagename)	text = ' IMMOBILE DOMAIN'
integer(i4b), parameter	nbditems = 5
character(len=lenbudtxt), dimension(nbditems)	budtxt

Detailed Description

The GwtIstModule is contains the GwtIstType, which is the derived type responsible for adding the effects of an immobile domain. In addition to representing transfer of mass between the mobile and immobile domain, the IST Package also represents the following processes within the immobile domain

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

Function/Subroutine Documentation

accumulate_budterm()

subroutine gwtistmodule::accumulate_budterm ( real(dp), dimension(:, :), intent(inout)  budterm, real(dp), dimension(:), intent(in)  ddterm, real(dp), intent(in)  cimnew, real(dp), intent(in)  cimold, real(dp), intent(in)  cnew, integer(i4b), intent(in)  idcy  )

private

This subroutine calculates and accumulates the immobile domain budget terms into the budterm accumulator

Parameters

[in]	ddterm	terms comprising the balance equation of the immobile domain
[in]	cimnew	immobile domain concenration at the end of this time step
[in]	cimold	immobile domain concentration at end of last time step
[in]	cnew	mobile domain concentration at the end of this time step
[in]	idcy	order of decay rate (0:none, 1:first, 2:zero)

Definition at line 1467 of file gwt-ist.f90.

    ! -- modules
    ! -- dummy
    real(DP), dimension(:, :), intent(inout) :: budterm !<
    real(DP), dimension(:), intent(in) :: ddterm !< terms comprising the balance equation of the immobile domain
    real(DP), intent(in) :: cimnew !< immobile domain concenration at the end of this time step
    real(DP), intent(in) :: cimold !< immobile domain concentration at end of last time step
    real(DP), intent(in) :: cnew !< mobile domain concentration at the end of this time step
    integer(I4B), intent(in) :: idcy !< order of decay rate (0:none, 1:first, 2:zero)
    ! -- local
    real(DP) :: rate
    integer(I4B) :: i
    !
    ! -- calculate STORAGE-AQUEOUS
    i = 1
    rate = -ddterm(1) * cimnew + ddterm(2) * cimold
    if (rate > dzero) then
      budterm(1, i) = budterm(1, i) + rate
    else
      budterm(2, i) = budterm(2, i) - rate
    end if
    !
    ! -- calculate STORAGE-SORBED
    i = 2
    rate = -ddterm(3) * cimnew + ddterm(4) * cimold
    if (rate > dzero) then
      budterm(1, i) = budterm(1, i) + rate
    else
      budterm(2, i) = budterm(2, i) - rate
    end if
    !
    ! -- calculate DECAY-AQUEOUS
    i = 3
    rate = dzero
    if (idcy == 1) then
      rate = -ddterm(5) * cimnew
    else if (idcy == 2) then
      rate = -ddterm(7)
    else
      rate = dzero
    end if
    if (rate > dzero) then
      budterm(1, i) = budterm(1, i) + rate
    else
      budterm(2, i) = budterm(2, i) - rate
    end if
    !
    ! -- calculate DECAY-SORBED
    i = 4
    if (idcy == 1) then
      rate = -ddterm(6) * cimnew
    else if (idcy == 2) then
      rate = -ddterm(8)
    else
      rate = dzero
    end if
    if (rate > dzero) then
      budterm(1, i) = budterm(1, i) + rate
    else
      budterm(2, i) = budterm(2, i) - rate
    end if
    !
    ! -- calculate MOBILE-DOMAIN
    i = 5
    rate = ddterm(9) * cnew - ddterm(9) * cimnew
    if (rate > dzero) then
      budterm(1, i) = budterm(1, i) + rate
    else
      budterm(2, i) = budterm(2, i) - rate
    end if
    !

Here is the caller graph for this function:

allocate_scalars()

subroutine gwtistmodule::allocate_scalars

(

class(gwtisttype)

this

)

Allocate and initialize package scalars.

Parameters

this	GwtIstType object

Definition at line 845 of file gwt-ist.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    use outputcontroldatamodule, only: ocd_cr
    ! -- dummy
    class(GwtIstType) :: this !< GwtIstType object
    ! -- local
    !
    ! -- call standard BndType allocate scalars
    call this%BndType%allocate_scalars()
    !
    ! -- Allocate
    call mem_allocate(this%icimout, 'ICIMOUT', this%memoryPath)
    call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
    call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
    call mem_allocate(this%ioutsorbate, 'IOUTSORBATE', this%memoryPath)
    call mem_allocate(this%isrb, 'ISRB', this%memoryPath)
    call mem_allocate(this%idcy, 'IDCY', this%memoryPath)
    call mem_allocate(this%kiter, 'KITER', this%memoryPath)
    !
    ! -- Initialize
    this%icimout = 0
    this%ibudgetout = 0
    this%ibudcsv = 0
    this%ioutsorbate = 0
    this%isrb = 0
    this%idcy = 0
    this%kiter = 0
    !
    ! -- Create the ocd object, which is used to manage printing and saving
    !    of the immobile domain concentrations
    call ocd_cr(this%ocd)

Here is the call graph for this function:

get_ddconc()

real(dp) function gwtistmodule::get_ddconc ( real(dp), dimension(:), intent(in)  ddterm, real(dp), intent(in)  f, real(dp), intent(in)  cimold, real(dp), intent(in)  cnew  )

private

This function calculates the concentration of the immobile domain.

Parameters

[in]	ddterm	terms comprising the balance equation of the immobile domain
[in]	f	the f term used to calculate the immobile domain concentration
[in]	cimold	immobile domain concentration at end of last time step
[in]	cnew	concentration of the mobile domain at the end of the time step

Returns

calculated concentration of the immobile domain

Definition at line 1445 of file gwt-ist.f90.

    ! -- dummy
    real(DP), dimension(:), intent(in) :: ddterm !< terms comprising the balance equation of the immobile domain
    real(DP), intent(in) :: f !< the f term used to calculate the immobile domain concentration
    real(DP), intent(in) :: cimold !< immobile domain concentration at end of last time step
    real(DP), intent(in) :: cnew !< concentration of the mobile domain at the end of the time step
    ! -- result
    real(DP) :: cimnew !< calculated concentration of the immobile domain
    ! -- local
    !
    ! -- calculate ddconc
    cimnew = (ddterm(2) + ddterm(4)) * cimold + ddterm(9) * cnew - ddterm(7) &
             - ddterm(8)
    cimnew = cimnew / f

Here is the caller graph for this function:

get_ddterm()

subroutine gwtistmodule::get_ddterm ( real(dp), intent(in)  thetaim, real(dp), intent(in)  vcell, real(dp), intent(in)  delt, real(dp), intent(in)  swtpdt, real(dp), intent(in)  volfracim, real(dp), intent(in)  rhobim, real(dp), intent(in)  kdnew, real(dp), intent(in)  kdold, real(dp), intent(in)  lambda1im, real(dp), intent(in)  lambda2im, real(dp), intent(in)  gamma1im, real(dp), intent(in)  gamma2im, real(dp), intent(in)  zetaim, real(dp), dimension(:), intent(inout)  ddterm, real(dp), intent(inout)  f  )

private

This subroutine calculates the immobile domain (or dual domain) terms. The resulting ddterm and f terms are described in the GWT model report. The terms are concentration coefficients used in the balance equation for the immobile domain.

Parameters

[in]	thetaim	immobile domain porosity
[in]	vcell	volume of cell
[in]	delt	length of time step
[in]	swtpdt	cell saturation at end of time step
[in]	volfracim	volume fraction of immobile domain
[in]	rhobim	bulk density for the immobile domain (fim * rhob)
[in]	kdnew	effective distribution coefficient for new time
[in]	kdold	effective distribution coefficient for old time
[in]	lambda1im	first-order decay rate in aqueous phase
[in]	lambda2im	first-order decay rate in sorbed phase
[in]	gamma1im	zero-order decay rate in aqueous phase
[in]	gamma2im	zero-order decay rate in sorbed phase
[in]	zetaim	transfer coefficient between mobile and immobile domains
[in,out]	ddterm	nine terms comprising the balance equation of the immobile domain
[in,out]	f	the f term used to calculate the immobile domain concentration

Definition at line 1372 of file gwt-ist.f90.

    ! -- dummy
    real(DP), intent(in) :: thetaim !< immobile domain porosity
    real(DP), intent(in) :: vcell !< volume of cell
    real(DP), intent(in) :: delt !< length of time step
    real(DP), intent(in) :: swtpdt !< cell saturation at end of time step
    real(DP), intent(in) :: volfracim !< volume fraction of immobile domain
    real(DP), intent(in) :: rhobim !< bulk density for the immobile domain (fim * rhob)
    real(DP), intent(in) :: kdnew !< effective distribution coefficient for new time
    real(DP), intent(in) :: kdold !< effective distribution coefficient for old time
    real(DP), intent(in) :: lambda1im !< first-order decay rate in aqueous phase
    real(DP), intent(in) :: lambda2im !< first-order decay rate in sorbed phase
    real(DP), intent(in) :: gamma1im !< zero-order decay rate in aqueous phase
    real(DP), intent(in) :: gamma2im !< zero-order decay rate in sorbed phase
    real(DP), intent(in) :: zetaim !< transfer coefficient between mobile and immobile domains
    real(DP), dimension(:), intent(inout) :: ddterm !< nine terms comprising the balance equation of the immobile domain
    real(DP), intent(inout) :: f !< the f term used to calculate the immobile domain concentration
    ! -- local
    real(DP) :: tled
    !
    ! -- initialize
    tled = done / delt
    !
    ! -- Calculate terms.  These terms correspond to the concentration
    !    coefficients in equation 7-4 of the GWT model report.  However,
    !    an updated equation is presented as 9-9 in the supplemental technical
    !    information guide (mf6suptechinfo.pdf)
    ddterm(1) = thetaim * vcell * tled
    ddterm(2) = thetaim * vcell * tled
    ddterm(3) = volfracim * rhobim * vcell * kdnew * tled
    ddterm(4) = volfracim * rhobim * vcell * kdold * tled
    ddterm(5) = thetaim * lambda1im * vcell
    ddterm(6) = lambda2im * volfracim * rhobim * kdnew * vcell
    ddterm(7) = thetaim * gamma1im * vcell
    ddterm(8) = gamma2im * volfracim * rhobim * vcell
    ddterm(9) = vcell * swtpdt * zetaim
    !
    ! -- calculate denominator term, f
    f = ddterm(1) + ddterm(3) + ddterm(5) + ddterm(6) + ddterm(9)

Here is the caller graph for this function:

get_hcofrhs()

subroutine gwtistmodule::get_hcofrhs ( real(dp), dimension(:), intent(in)  ddterm, real(dp), intent(in)  f, real(dp), intent(in)  cimold, real(dp), intent(inout)  hcof, real(dp), intent(inout)  rhs  )

private

This subroutine calculates the hcof and rhs terms that must be added to the solution system of equations

Parameters

[in]	ddterm	terms comprising the balance equation of the immobile domain
[in]	f	the f term used to calculate the immobile domain concentration
[in]	cimold	immobile domain concentration at end of last time step
[in,out]	hcof	calculated contribution for the a-matrix diagonal position
[in,out]	rhs	calculated contribution for the solution right-hand side

Definition at line 1422 of file gwt-ist.f90.

    ! -- dummy
    real(DP), dimension(:), intent(in) :: ddterm !< terms comprising the balance equation of the immobile domain
    real(DP), intent(in) :: f !< the f term used to calculate the immobile domain concentration
    real(DP), intent(in) :: cimold !< immobile domain concentration at end of last time step
    real(DP), intent(inout) :: hcof !< calculated contribution for the a-matrix diagonal position
    real(DP), intent(inout) :: rhs !< calculated contribution for the solution right-hand side
    !
    ! -- calculate hcof
    hcof = ddterm(9)**2 / f - ddterm(9)
    !
    ! -- calculate rhs, and switch the sign because this term needs to
    !    be moved to the left hand side
    rhs = (ddterm(2) + ddterm(4)) * cimold - ddterm(7) - ddterm(8)
    rhs = rhs * ddterm(9) / f
    rhs = -rhs

Here is the caller graph for this function:

get_thetaim()

real(dp) function gwtistmodule::get_thetaim	(	class(gwtisttype)	this,
		integer(i4b), intent(in)	node
	)

Calculate and return thetaim, volume of immobile voids per aquifer volume

Parameters

	this	GwtIstType object
[in]	node	node number

Definition at line 1353 of file gwt-ist.f90.

    ! -- modules
    ! -- dummy
    class(GwtIstType) :: this !< GwtIstType object
    integer(I4B), intent(in) :: node !< node number
    ! -- return
    real(DP) :: thetaim
    !
    thetaim = this%volfrac(node) * this%porosity(node)

ist_ad()

subroutine gwtistmodule::ist_ad ( class(gwtisttype)  this )

private

Advance the IST Package and handle the adaptive time stepping feature by copying from new to old or old to new accordingly

Parameters

this	GwtIstType object

Definition at line 233 of file gwt-ist.f90.

    ! -- modules
    use simvariablesmodule, only: ifailedstepretry
    ! -- dummy variables
    class(GwtIstType) :: this !< GwtIstType object
    ! -- local variables
    integer(I4B) :: n
    !
    ! -- Call parent advance
    call this%BndType%bnd_ad()
    !
    ! -- set independent kiter counter to zero
    this%kiter = 0
    !
    ! -- copy cimnew into cimold or vice versa if this is a repeat of
    !    a failed time step
    if (ifailedstepretry == 0) then
      do n = 1, this%dis%nodes
        this%cimold(n) = this%cimnew(n)
      end do
    else
      do n = 1, this%dis%nodes
        this%cimnew(n) = this%cimold(n)
      end do
    end if

ist_allocate_arrays()

subroutine gwtistmodule::ist_allocate_arrays

(

class(gwtisttype), intent(inout)

this

)

Allocate and initialize package arrays.

Parameters

[in,out]

this

GwtIstType object

Definition at line 884 of file gwt-ist.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(GwtIstType), intent(inout) :: this !< GwtIstType object
    ! -- local
    integer(I4B) :: n
    !
    ! -- call standard BndType allocate scalars
    !    nbound and maxbound are 0 in order to keep memory footprint low
    call this%BndType%allocate_arrays()
    !
    ! -- allocate ist arrays of size nodes
    call mem_allocate(this%strg, this%dis%nodes, 'STRG', this%memoryPath)
    call mem_allocate(this%cim, this%dis%nodes, 'CIM', this%memoryPath)
    call mem_allocate(this%cimnew, this%dis%nodes, 'CIMNEW', this%memoryPath)
    call mem_allocate(this%cimold, this%dis%nodes, 'CIMOLD', this%memoryPath)
    call mem_allocate(this%porosity, this%dis%nodes, 'POROSITY', this%memoryPath)
    call mem_allocate(this%zetaim, this%dis%nodes, 'ZETAIM', this%memoryPath)
    call mem_allocate(this%volfrac, this%dis%nodes, 'VOLFRAC', this%memoryPath)
    if (this%isrb == 0) then
      call mem_allocate(this%bulk_density, 1, 'BULK_DENSITY', this%memoryPath)
      call mem_allocate(this%distcoef, 1, 'DISTCOEF', this%memoryPath)
      call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
      call mem_allocate(this%cimsrb, 1, 'SORBATE', this%memoryPath)
    else
      call mem_allocate(this%bulk_density, this%dis%nodes, 'BULK_DENSITY', &
                        this%memoryPath)
      call mem_allocate(this%distcoef, this%dis%nodes, 'DISTCOEF', &
                        this%memoryPath)
      call mem_allocate(this%cimsrb, this%dis%nodes, 'SORBATE', &
                        this%memoryPath)
      if (this%isrb == 1) then
        call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
      else
        call mem_allocate(this%sp2, this%dis%nodes, 'SP2', this%memoryPath)
      end if
    end if
    if (this%idcy == 0) then
      call mem_allocate(this%decay, 1, 'DECAY', this%memoryPath)
      call mem_allocate(this%decaylast, 1, 'DECAYLAST', this%memoryPath)
    else
      call mem_allocate(this%decay, this%dis%nodes, 'DECAY', this%memoryPath)
      call mem_allocate(this%decaylast, this%dis%nodes, 'DECAYLAST', &
                        this%memoryPath)
    end if
    if (this%isrb == 0 .and. this%idcy == 0) then
      call mem_allocate(this%decayslast, 1, 'DECAYSLAST', this%memoryPath)
    else
      call mem_allocate(this%decayslast, this%dis%nodes, 'DECAYSLAST', &
                        this%memoryPath)
    end if
    call mem_allocate(this%decay_sorbed, 1, 'DECAY_SORBED', this%memoryPath)
    !
    ! -- initialize
    do n = 1, this%dis%nodes
      this%strg(n) = dzero
      this%cim(n) = dzero
      this%cimnew(n) = dzero
      this%cimold(n) = dzero
      this%porosity(n) = dzero
      this%zetaim(n) = dzero
      this%volfrac(n) = dzero
    end do
    do n = 1, size(this%bulk_density)
      this%bulk_density(n) = dzero
      this%distcoef(n) = dzero
      this%cimsrb(n) = dzero
    end do
    do n = 1, size(this%sp2)
      this%sp2(n) = dzero
    end do
    do n = 1, size(this%decay)
      this%decay(n) = dzero
      this%decaylast(n) = dzero
    end do
    do n = 1, size(this%decayslast)
      this%decayslast(n) = dzero
    end do
    !
    ! -- Set pointers
    this%ocd%dis => this%dis

ist_ar()

subroutine gwtistmodule::ist_ar ( class(gwtisttype), intent(inout)  this )

private

Method to allocate and read static data for the package.

Parameters

[in,out]

this

GwtIstType object

Definition at line 162 of file gwt-ist.f90.

    ! -- modules
    use simmodule, only: store_error, count_errors
    use budgetmodule, only: budget_cr
    ! -- dummy
    class(GwtIstType), intent(inout) :: this !< GwtIstType object
    ! -- local
    integer(I4B) :: n
    !
    ! -- Allocate arrays
    call this%ist_allocate_arrays()
    !
    ! -- Now that arrays are allocated, check in the cimnew array to
    !    the output control manager for subsequent printing/saving
    call this%ocd%init_dbl('CIM', this%cimnew, this%dis, 'PRINT LAST ', &
                           'COLUMNS 10 WIDTH 11 DIGITS 4 GENERAL ', &
                           this%iout, dhnoflo)
    !
    ! -- read the data block
    call this%read_data()
    !
    ! -- set cimnew to the cim start values read from input
    do n = 1, this%dis%nodes
      this%cimnew(n) = this%cim(n)
    end do
    !
    ! -- add volfrac to the volfracim accumulator in mst package
    call this%mst%addto_volfracim(this%volfrac)
    !
    ! -- setup the immobile domain budget
    call budget_cr(this%budget, this%memoryPath)
    call this%budget%budget_df(nbditems, 'MASS', 'M', bdzone=this%packName)
    call this%budget%set_ibudcsv(this%ibudcsv)
    !
    ! -- Perform a check to ensure that sorption and decay are set
    !    consistently between the MST and IST packages.
    if (this%idcy /= this%mst%idcy) then
      call store_error('DECAY must be activated consistently between the &
        &MST and IST Packages.  Activate or deactivate DECAY for both &
        &Packages.')
    end if
    if (this%isrb /= this%mst%isrb) then
      call store_error('Sorption is active for the IST Package but it is not &
        &compatible with the sorption option selected for the MST Package.  &
        &If sorption is active for the IST Package, then the same type of &
        &sorption (LINEAR, FREUNDLICH, or LANGMUIR) must &
        &be specified in the options block of the MST Package.')
    end if
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if

Here is the call graph for this function:

ist_bd()

subroutine gwtistmodule::ist_bd ( class(gwtisttype)  this, type(budgettype), intent(inout)  model_budget  )

private

Add the flow between IST package and the model (ratin and ratout) to the model budget.

Parameters

	this	GwtIstType object
[in,out]	model_budget	model budget object

Definition at line 590 of file gwt-ist.f90.

    ! -- modules
    use tdismodule, only: delt
    use budgetmodule, only: budgettype, rate_accumulator
    ! -- dummy
    class(GwtIstType) :: this !< GwtIstType object
    type(BudgetType), intent(inout) :: model_budget !< model budget object
    ! -- local
    real(DP) :: ratin
    real(DP) :: ratout
    integer(I4B) :: isuppress_output
    isuppress_output = 0
    call rate_accumulator(this%strg(:), ratin, ratout)
    call model_budget%addentry(ratin, ratout, delt, this%text, &
                               isuppress_output, this%packName)

Here is the call graph for this function:

ist_calc_csrb()

subroutine gwtistmodule::ist_calc_csrb	(	class(gwtisttype)	this,
		real(dp), dimension(:), intent(in)	cim
	)

Parameters

	this	GwtMstType object
[in]	cim	immobile domain aqueous concentration at end of this time step

Definition at line 557 of file gwt-ist.f90.

    ! -- dummy
    class(GwtIstType) :: this !< GwtMstType object
    real(DP), intent(in), dimension(:) :: cim !< immobile domain aqueous concentration at end of this time step
    ! -- local
    integer(I4B) :: n
    real(DP) :: distcoef
    real(DP) :: csrb
 
    ! Calculate sorbed concentration
    do n = 1, size(cim)
      csrb = dzero
      if (this%ibound(n) > 0) then
        distcoef = this%distcoef(n)
        if (this%isrb == 1) then
          csrb = cim(n) * distcoef
        else if (this%isrb == 2) then
          csrb = get_freundlich_conc(cim(n), distcoef, this%sp2(n))
        else if (this%isrb == 3) then
          csrb = get_langmuir_conc(cim(n), distcoef, this%sp2(n))
        end if
      end if
      this%cimsrb(n) = csrb
    end do
 

Here is the call graph for this function:

ist_cq()

subroutine gwtistmodule::ist_cq	(	class(gwtisttype), intent(inout)	this,
		real(dp), dimension(:), intent(in)	x,
		real(dp), dimension(:), intent(inout), contiguous	flowja,
		integer(i4b), intent(in), optional	iadv
	)

Parameters

[in,out]	this	GwtIstType object
[in]	x	current dependent-variable value
[in,out]	flowja	flow between two connected control volumes
[in]	iadv	flag that indicates if this is an advance package

Definition at line 400 of file gwt-ist.f90.

    ! modules
    use tdismodule, only: delt
    use constantsmodule, only: dzero
    ! dummy
    class(GwtIstType), intent(inout) :: this !< GwtIstType object
    real(DP), dimension(:), intent(in) :: x !< current dependent-variable value
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
    integer(I4B), optional, intent(in) :: iadv !< flag that indicates if this is an advance package
    ! local
    integer(I4B) :: idiag
    integer(I4B) :: n
    real(DP) :: rate
    real(DP) :: swt, swtpdt
    real(DP) :: hhcof, rrhs
    real(DP) :: vcell
    real(DP) :: thetaim
    real(DP) :: zetaim
    real(DP) :: kdnew
    real(DP) :: kdold
    real(DP) :: volfracim
    real(DP) :: rhobim
    real(DP) :: lambda1im
    real(DP) :: lambda2im
    real(DP) :: gamma1im
    real(DP) :: gamma2im
    real(DP) :: cimnew
    real(DP) :: cimold
    real(DP) :: f
    real(DP) :: cimsrbold
    real(DP) :: cimsrbnew
    real(DP), dimension(9) :: ddterm
    ! formats
 
    ! initialize
    this%budterm(:, :) = dzero
 
    ! Calculate immobile domain transfer rate
    do n = 1, this%dis%nodes
 
      ! skip if transport inactive
      rate = dzero
      cimnew = dzero
      if (this%ibound(n) > 0) then
 
        ! calculate new and old water volumes
        vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
        swtpdt = this%fmi%gwfsat(n)
        swt = this%fmi%gwfsatold(n, delt)
        thetaim = this%get_thetaim(n)
 
        ! set exchange coefficient
        zetaim = this%zetaim(n)
 
        ! Calculate exchange with immobile domain
        rate = dzero
        hhcof = dzero
        rrhs = dzero
        kdnew = dzero
        kdold = dzero
        volfracim = dzero
        rhobim = dzero
        lambda1im = dzero
        lambda2im = dzero
        gamma1im = dzero
        gamma2im = dzero
        if (this%idcy == 1) lambda1im = this%decay(n)
        if (this%idcy == 2) then
          gamma1im = get_zero_order_decay(this%decay(n), this%decaylast(n), 0, &
                                          this%cimold(n), this%cimnew(n), delt)
        end if
 
        ! setup sorption variables
        if (this%isrb > 0) then
 
          ! initialize sorption variables
          volfracim = this%volfrac(n)
          rhobim = this%bulk_density(n)
 
          ! set isotherm dependent sorption variables
          select case (this%isrb)
          case (1)
            ! linear
            kdnew = this%distcoef(n)
            kdold = this%distcoef(n)
            cimsrbnew = this%cimnew(n) * kdnew
            cimsrbold = this%cimold(n) * kdold
          case (2)
            ! freundlich
            kdnew = get_freundlich_kd(this%cimnew(n), this%distcoef(n), &
                                      this%sp2(n))
            kdold = get_freundlich_kd(this%cimold(n), this%distcoef(n), &
                                      this%sp2(n))
            cimsrbnew = get_freundlich_conc(this%cimnew(n), this%distcoef(n), &
                                            this%sp2(n))
            cimsrbold = get_freundlich_conc(this%cimold(n), this%distcoef(n), &
                                            this%sp2(n))
          case (3)
            ! langmuir
            kdnew = get_langmuir_kd(this%cimnew(n), this%distcoef(n), &
                                    this%sp2(n))
            kdold = get_langmuir_kd(this%cimold(n), this%distcoef(n), &
                                    this%sp2(n))
            cimsrbnew = get_langmuir_conc(this%cimnew(n), this%distcoef(n), &
                                          this%sp2(n))
            cimsrbold = get_langmuir_conc(this%cimold(n), this%distcoef(n), &
                                          this%sp2(n))
          end select
 
          ! set decay of sorbed solute parameters
          if (this%idcy == 1) then
            lambda2im = this%decay_sorbed(n)
          else if (this%idcy == 2) then
            gamma2im = get_zero_order_decay(this%decay_sorbed(n), &
                                            this%decayslast(n), &
                                            0, cimsrbold, &
                                            cimsrbnew, delt)
          end if
        end if
 
        ! calculate the terms and then get the hcof and rhs contributions
        call get_ddterm(thetaim, vcell, delt, swtpdt, &
                        volfracim, rhobim, kdnew, kdold, lambda1im, lambda2im, &
                        gamma1im, gamma2im, zetaim, ddterm, f)
        cimold = this%cimold(n)
        call get_hcofrhs(ddterm, f, cimold, hhcof, rrhs)
 
        ! calculate rate from hcof and rhs
        rate = hhcof * x(n) - rrhs
 
        ! calculate immobile domain concentration
        cimnew = get_ddconc(ddterm, f, cimold, x(n))
 
        ! accumulate the budget terms
        call accumulate_budterm(this%budterm, ddterm, cimnew, cimold, x(n), &
                                this%idcy)
      end if
 
      ! store rate and add to flowja
      this%strg(n) = rate
      idiag = this%dis%con%ia(n)
      flowja(idiag) = flowja(idiag) + rate
 
      ! store immobile domain concentration
      this%cimnew(n) = cimnew
 
    end do
 
    ! calculate csrb
    if (this%isrb /= 0) then
      call this%ist_calc_csrb(this%cimnew)
    end if
 

Here is the call graph for this function:

ist_create()

subroutine, public gwtistmodule::ist_create	(	class(bndtype), pointer	packobj,
		integer(i4b), intent(in)	id,
		integer(i4b), intent(in)	ibcnum,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout,
		character(len=*), intent(in)	namemodel,
		character(len=*), intent(in)	pakname,
		type(tspfmitype), pointer	fmi,
		type(gwtmsttype), pointer	mst
	)

Create a new IST object

Parameters

	packobj	BndType pointer that will point to new IST Package
[in]	id	name of the model
[in]	ibcnum	consecutive package number
[in]	inunit	unit number of package input file
[in]	iout	unit number of model listing file
[in]	namemodel	name of the model
[in]	pakname	name of the package

Definition at line 115 of file gwt-ist.f90.

    ! -- dummy
    class(BndType), pointer :: packobj !< BndType pointer that will point to new IST Package
    integer(I4B), intent(in) :: id !< name of the model
    integer(I4B), intent(in) :: ibcnum !< consecutive package number
    integer(I4B), intent(in) :: inunit !< unit number of package input file
    integer(I4B), intent(in) :: iout !< unit number of model listing file
    character(len=*), intent(in) :: namemodel !< name of the model
    character(len=*), intent(in) :: pakname !< name of the package
    ! -- local
    type(GwtIstType), pointer :: istobj
    type(TspFmiType), pointer :: fmi
    type(GwtMstType), pointer :: mst
    !
    ! -- allocate the object and assign values to object variables
    allocate (istobj)
    packobj => istobj
    !
    ! -- create name and memory path
    call packobj%set_names(ibcnum, namemodel, pakname, ftype)
    packobj%text = text
    !
    ! -- allocate scalars
    call packobj%allocate_scalars()
    !
    ! -- initialize package
    call packobj%pack_initialize()
    !
    ! -- store values
    packobj%inunit = inunit
    packobj%iout = iout
    packobj%id = id
    packobj%ibcnum = ibcnum
    packobj%ncolbnd = 1
    packobj%iscloc = 1
    !
    ! -- Point IST specific variables
    istobj%fmi => fmi
    istobj%mst => mst

Here is the caller graph for this function:

ist_da()

subroutine gwtistmodule::ist_da

(

class(gwtisttype)

this

)

Deallocate package scalars and arrays.

Parameters

this	GwtIstType object

Definition at line 794 of file gwt-ist.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(GwtIstType) :: this !< GwtIstType object
    !
    ! -- Deallocate arrays if package was active
    if (this%inunit > 0) then
      call mem_deallocate(this%icimout)
      call mem_deallocate(this%ibudgetout)
      call mem_deallocate(this%ibudcsv)
      call mem_deallocate(this%ioutsorbate)
      call mem_deallocate(this%idcy)
      call mem_deallocate(this%isrb)
      call mem_deallocate(this%kiter)
      call mem_deallocate(this%cim)
      call mem_deallocate(this%cimnew)
      call mem_deallocate(this%cimold)
      call mem_deallocate(this%cimsrb)
      call mem_deallocate(this%zetaim)
      call mem_deallocate(this%porosity)
      call mem_deallocate(this%volfrac)
      call mem_deallocate(this%bulk_density)
      call mem_deallocate(this%distcoef)
      call mem_deallocate(this%sp2)
      call mem_deallocate(this%decay)
      call mem_deallocate(this%decaylast)
      call mem_deallocate(this%decayslast)
      call mem_deallocate(this%decay_sorbed)
      call mem_deallocate(this%strg)
      this%fmi => null()
      this%mst => null()
    end if
    !
    ! -- Scalars
    !
    ! -- Objects
    call this%budget%budget_da()
    deallocate (this%budget)
    call this%ocd%ocd_da()
    deallocate (this%ocd)
    !
    ! -- deallocate parent
    call this%BndType%bnd_da()

ist_fc()

subroutine gwtistmodule::ist_fc	(	class(gwtisttype)	this,
		real(dp), dimension(:), intent(inout)	rhs,
		integer(i4b), dimension(:), intent(in)	ia,
		integer(i4b), dimension(:), intent(in)	idxglo,
		class(matrixbasetype), pointer	matrix_sln
	)

Parameters

	this	GwtIstType object
[in,out]	rhs	right-hand side vector for model
[in]	ia	solution CRS row pointers
[in]	idxglo	mapping vector for model (local) to solution (global)
	matrix_sln	solution coefficient matrix

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

    ! modules
    use tdismodule, only: delt
    ! dummy
    class(GwtIstType) :: this !< GwtIstType object
    real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model
    integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers
    integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global)
    class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
    ! local
    integer(I4B) :: n, idiag
    real(DP) :: tled
    real(DP) :: hhcof, rrhs
    real(DP) :: swt, swtpdt
    real(DP) :: vcell
    real(DP) :: thetaim
    real(DP) :: zetaim
    real(DP) :: kdnew
    real(DP) :: kdold
    real(DP) :: volfracim
    real(DP) :: rhobim
    real(DP) :: lambda1im
    real(DP) :: lambda2im
    real(DP) :: gamma1im
    real(DP) :: gamma2im
    real(DP) :: cimold
    real(DP) :: f
    real(DP) :: cimsrbold
    real(DP) :: cimsrbnew
    real(DP), dimension(9) :: ddterm
 
    ! set variables
    tled = done / delt
    this%kiter = this%kiter + 1
 
    ! loop through each node and calculate immobile domain 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)
      swt = this%fmi%gwfsatold(n, delt)
      thetaim = this%get_thetaim(n)
      idiag = ia(n)
 
      ! set exchange coefficient
      zetaim = this%zetaim(n)
 
      ! Add dual domain mass transfer contributions to rhs and hcof
      ! dcimsrbdc = DZERO
      kdnew = dzero
      kdold = dzero
      volfracim = dzero
      rhobim = dzero
      lambda1im = dzero
      lambda2im = dzero
      gamma1im = dzero
      gamma2im = dzero
 
      ! set variables for decay of aqueous solute
      if (this%idcy == 1) lambda1im = this%decay(n)
      if (this%idcy == 2) then
        gamma1im = get_zero_order_decay(this%decay(n), this%decaylast(n), &
                                        this%kiter, this%cimold(n), &
                                        this%cimnew(n), delt)
        this%decaylast(n) = gamma1im
      end if
 
      ! setup sorption variables
      if (this%isrb > 0) then
 
        ! initialize sorption variables
        volfracim = this%volfrac(n)
        rhobim = this%bulk_density(n)
 
        ! set isotherm dependent sorption variables
        select case (this%isrb)
        case (1)
          ! linear
          kdnew = this%distcoef(n)
          kdold = this%distcoef(n)
          cimsrbnew = this%cimnew(n) * kdnew
          cimsrbold = this%cimold(n) * kdold
        case (2)
          ! freundlich
          kdnew = get_freundlich_kd(this%cimnew(n), this%distcoef(n), &
                                    this%sp2(n))
          kdold = get_freundlich_kd(this%cimold(n), this%distcoef(n), &
                                    this%sp2(n))
          cimsrbnew = get_freundlich_conc(this%cimnew(n), this%distcoef(n), &
                                          this%sp2(n))
          cimsrbold = get_freundlich_conc(this%cimold(n), this%distcoef(n), &
                                          this%sp2(n))
        case (3)
          ! langmuir
          kdnew = get_langmuir_kd(this%cimnew(n), this%distcoef(n), &
                                  this%sp2(n))
          kdold = get_langmuir_kd(this%cimold(n), this%distcoef(n), &
                                  this%sp2(n))
          cimsrbnew = get_langmuir_conc(this%cimnew(n), this%distcoef(n), &
                                        this%sp2(n))
          cimsrbold = get_langmuir_conc(this%cimold(n), this%distcoef(n), &
                                        this%sp2(n))
        end select
 
        ! set decay of sorbed solute parameters
        if (this%idcy == 1) then
          lambda2im = this%decay_sorbed(n)
        else if (this%idcy == 2) then
          gamma2im = get_zero_order_decay(this%decay_sorbed(n), &
                                          this%decayslast(n), &
                                          this%kiter, cimsrbold, &
                                          cimsrbnew, delt)
          this%decayslast(n) = gamma2im
        end if
      end if
 
      ! calculate dual domain terms and get the hcof and rhs contributions
      call get_ddterm(thetaim, vcell, delt, swtpdt, &
                      volfracim, rhobim, kdnew, kdold, lambda1im, lambda2im, &
                      gamma1im, gamma2im, zetaim, ddterm, f)
      cimold = this%cimold(n)
      call get_hcofrhs(ddterm, f, cimold, hhcof, rrhs)
 
      ! update solution accumulators
      call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
      rhs(n) = rhs(n) + rrhs
 
    end do
 

Here is the call graph for this function:

ist_ot_bdsummary()

subroutine gwtistmodule::ist_ot_bdsummary ( class(gwtisttype)  this, integer(i4b), intent(in)  kstp, integer(i4b), intent(in)  kper, integer(i4b), intent(in)  iout, integer(i4b), intent(in)  ibudfl  )

private

Output advanced boundary package budget summary. This method only needs to be overridden for advanced packages that save budget summaries to the model listing file.

Parameters

	this	GwtIstType object
[in]	kstp	time step number
[in]	kper	period number
[in]	iout	flag and unit number for the model listing file
[in]	ibudfl	flag indicating budget should be written

Definition at line 763 of file gwt-ist.f90.

    ! -- modules
    use tdismodule, only: delt, totim
    ! -- dummy variables
    class(GwtIstType) :: this !< GwtIstType object
    integer(I4B), intent(in) :: kstp !< time step number
    integer(I4B), intent(in) :: kper !< period number
    integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
    integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
    ! -- local
    integer(I4B) :: isuppress_output = 0
    !
    ! -- Fill budget terms
    call this%budget%reset()
    call this%budget%addentry(this%budterm, delt, budtxt, isuppress_output)
    !
    ! -- Write budget to list file
    call this%budget%finalize_step(delt)
    if (ibudfl /= 0) then
      call this%budget%budget_ot(kstp, kper, iout)
    end if
    !
    ! -- Write budget csv
    call this%budget%writecsv(totim)

ist_ot_dv()

subroutine gwtistmodule::ist_ot_dv	(	class(gwtisttype)	this,
		integer(i4b), intent(in)	idvsave,
		integer(i4b), intent(in)	idvprint
	)

Parameters

	this	BndType object
[in]	idvsave	flag and unit number for dependent-variable output
[in]	idvprint	flag indicating if dependent-variable should be written to the model listing file

Definition at line 676 of file gwt-ist.f90.

    ! dummy variables
    class(GwtIstType) :: this !< BndType object
    integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output
    integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file
 
    call this%output_immobile_concentration(idvsave, idvprint)
    call this%output_immobile_sorbate_concentration(idvsave, idvprint)
 

ist_ot_model_flows()

subroutine gwtistmodule::ist_ot_model_flows	(	class(gwtisttype)	this,
		integer(i4b), intent(in)	icbcfl,
		integer(i4b), intent(in)	ibudfl,
		integer(i4b), intent(in)	icbcun,
		integer(i4b), dimension(:), intent(in), optional	imap
	)

Output flow terms between the IST package and model to a binary file and/or print flows to the model listing file.

Parameters

	this	GwtIstType object
[in]	icbcfl	flag for cell-by-cell output
[in]	ibudfl	flag indication if cell-by-cell data should be saved
[in]	icbcun	unit number for cell-by-cell output
[in]	imap	mapping vector

Definition at line 613 of file gwt-ist.f90.

    ! -- modules
    use constantsmodule, only: dzero
    ! -- dummy
    class(GwtIstType) :: this !< GwtIstType object
    integer(I4B), intent(in) :: icbcfl !< flag for cell-by-cell output
    integer(I4B), intent(in) :: ibudfl !< flag indication if cell-by-cell data should be saved
    integer(I4B), intent(in) :: icbcun !< unit number for cell-by-cell output
    integer(I4B), dimension(:), optional, intent(in) :: imap !< mapping vector
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: ibinun
    integer(I4B) :: nbound
    integer(I4B) :: naux
    real(DP) :: rate
    real(DP), dimension(0) :: auxrow
    !
    ! -- 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 cell-by-cell flows will be saved as a list, write header.
    if (ibinun /= 0) then
      nbound = this%dis%nodes
      naux = 0
      call this%dis%record_srcdst_list_header(this%text, this%name_model, &
                                              this%name_model, this%name_model, &
                                              this%packName, naux, this%auxname, &
                                              ibinun, nbound, this%iout)
    end if
    !
    ! -- Calculate immobile domain rhs and hcof
    do n = 1, this%dis%nodes
      !
      ! -- skip if transport inactive
      rate = dzero
      if (this%ibound(n) > 0) then
        !
        ! -- set rate from this%strg
        rate = this%strg(n)
      end if
      !
      ! -- If saving cell-by-cell flows in list, write flow
      if (ibinun /= 0) then
        call this%dis%record_mf6_list_entry(ibinun, n, n, rate, &
                                            naux, auxrow, &
                                            olconv=.true., &
                                            olconv2=.true.)
      end if
      !
    end do

ist_read_dimensions()

subroutine gwtistmodule::ist_read_dimensions

(

class(gwtisttype), intent(inout)

this

)

Read dimensions for package.

Parameters

[in,out]

this

GwtIstType object

Definition at line 1113 of file gwt-ist.f90.

    ! -- dummy
    class(GwtIstType), intent(inout) :: this !< GwtIstType object
    ! -- local
    ! -- format

ist_rp()

subroutine gwtistmodule::ist_rp

(

class(gwtisttype), intent(inout)

this

)

Method to read and prepare package data

Parameters

[in,out]

this

GwtIstType object

Definition at line 220 of file gwt-ist.f90.

    ! -- dummy
    class(GwtIstType), intent(inout) :: this !< GwtIstType object
    ! -- local
    ! -- format

output_immobile_concentration()

subroutine gwtistmodule::output_immobile_concentration ( class(gwtisttype)  this, integer(i4b), intent(in)  idvsave, integer(i4b), intent(in)  idvprint  )

private

Parameters

	this	BndType object
[in]	idvsave	flag and unit number for dependent-variable output
[in]	idvprint	flag indicating if dependent-variable should be written to the model listing file

Definition at line 689 of file gwt-ist.f90.

    ! modules
    use tdismodule, only: kstp, endofperiod
    ! dummy variables
    class(GwtIstType) :: this !< BndType object
    integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output
    integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file
    ! local
    integer(I4B) :: ipflg
    integer(I4B) :: ibinun
    !
    ! Save cim to a binary file. ibinun is a flag where 1 indicates that
    ! cim should be written to a binary file if a binary file is open for it.
    ipflg = 0
    ibinun = 1
    if (idvsave == 0) ibinun = 0
    if (ibinun /= 0) then
      call this%ocd%ocd_ot(ipflg, kstp, endofperiod, this%iout, &
                           iprint_opt=0, isav_opt=ibinun)
    end if
    !
    ! Print immobile domain concentrations to listing file
    if (idvprint /= 0) then
      call this%ocd%ocd_ot(ipflg, kstp, endofperiod, this%iout, &
                           iprint_opt=idvprint, isav_opt=0)
    end if
 

output_immobile_sorbate_concentration()

subroutine gwtistmodule::output_immobile_sorbate_concentration	(	class(gwtisttype)	this,
		integer(i4b), intent(in)	idvsave,
		integer(i4b), intent(in)	idvprint
	)

Parameters

	this	BndType object
[in]	idvsave	flag and unit number for dependent-variable output
[in]	idvprint	flag indicating if dependent-variable should be written to the model listing file

Definition at line 720 of file gwt-ist.f90.

    ! modules
    ! dummy
    class(GwtIstType) :: this !< BndType object
    integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output
    integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file
    ! local
    character(len=1) :: cdatafmp = ' ', editdesc = ' '
    integer(I4B) :: ibinun
    integer(I4B) :: iprint, nvaluesp, nwidthp
    real(DP) :: dinact
 
    ! Save cimsrb to a binary file. ibinun is a flag where 1 indicates that
    ! cim should be written to a binary file if a binary file is open for it.
    ! 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%cimsrb, this%iout, iprint, ibinun, &
                                   '         SORBATE', cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      end if
    end if
 

read_data()

subroutine gwtistmodule::read_data ( class(gwtisttype)  this )

private

Read data for package.

Parameters

this	GwtIstType object

Definition at line 1125 of file gwt-ist.f90.

    ! -- modules
    use constantsmodule, only: linelength
    use simmodule, only: store_error, count_errors
    use memorymanagermodule, only: mem_reallocate, mem_reassignptr
    ! -- dummy
    class(GwtIstType) :: this !< GwtIstType object
    ! -- local
    character(len=LINELENGTH) :: errmsg, keyword
    character(len=:), allocatable :: line
    integer(I4B) :: istart, istop, lloc, ierr
    logical :: isfound, endOfBlock
    logical, dimension(9) :: lname
    character(len=24), dimension(9) :: aname
    ! -- formats
    ! -- data
    data aname(1)/'            BULK DENSITY'/
    data aname(2)/'DISTRIBUTION COEFFICIENT'/
    data aname(3)/'              DECAY RATE'/
    data aname(4)/'       DECAY SORBED RATE'/
    data aname(5)/'   INITIAL IMMOBILE CONC'/
    data aname(6)/'  FIRST ORDER TRANS RATE'/
    data aname(7)/'IMMOBILE DOMAIN POROSITY'/
    data aname(8)/'IMMOBILE VOLUME FRACTION'/
    data aname(9)/'   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 ('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(1))
          lname(1) = .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(2))
          lname(2) = .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(3))
          lname(3) = .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(4))
          lname(4) = .true.
        case ('CIM')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%cim, &
                                        aname(5))
          lname(5) = .true.
        case ('ZETAIM')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%zetaim, &
                                        aname(6))
          lname(6) = .true.
        case ('POROSITY')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%porosity, &
                                        aname(7))
          lname(7) = .true.
        case ('VOLFRAC')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%volfrac, &
                                        aname(8))
          lname(8) = .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(9))
          lname(9) = .true.
        case ('THETAIM')
          write (errmsg, '(a)') &
            'THETAIM is no longer supported. See Chapter 9 in &
            &mf6suptechinfo.pdf for revised parameterization of mobile and &
            &immobile domain simulations.'
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        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 sorption variables
    if (this%isrb > 0) then
      if (.not. lname(1)) 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(2)) 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 .and. .not. lname(9)) then
        write (errmsg, '(a)') 'Nonlinear sorption is active but SP2 &
          &not specified.  SP2 must be specified in GRIDDATA block when &
          &FREUNDLICH or LANGMUIR sorption is active.'
        call store_error(errmsg)
      end if
    else
      if (lname(1)) then
        write (this%iout, '(1x,a)') 'Warning.  Sorption is not active but &
          &BULK_DENSITY was specified.  BULK_DENSITY will have no affect on &
          &simulation results.'
      end if
      if (lname(2)) then
        write (this%iout, '(1x,a)') 'Warning.  Sorption is not active but &
          &distribution coefficient was specified.  DISTCOEF will have &
          &no affect on simulation results.'
      end if
      if (lname(9)) then
        write (this%iout, '(1x,a)') 'Warning.  Sorption is not active but &
          &SP2 was specified.  SP2 will have no affect on simulation results.'
      end if
    end if
    !
    ! -- Check for required decay/production rate coefficients
    if (this%idcy > 0) then
      if (.not. lname(3)) then
        write (errmsg, '(a)') 'First- or zero-order decay is &
          &active but the first rate coefficient was not specified.  &
          &Decay must be specified in GRIDDATA block.'
        call store_error(errmsg)
      end if
      if (.not. lname(4)) then
        !
        ! -- If DECAY_SORBED not specified and sorption is active, then set
        !    decay_sorbed equal to decay
        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(3)) then
        write (this%iout, '(1x,a)') 'Warning.  First- or zero-order decay &
          &is not active but DECAY was specified.  DECAY will &
          &have no affect on simulation results.'
      end if
      if (lname(4)) then
        write (this%iout, '(1x,a)') 'Warning.  First- or zero-order decay &
          &is not active but DECAY_SORBED was specified.  &
          &DECAY_SORBED will have no affect on simulation &
          &results.'
      end if
    end if
    !
    ! -- Check for required dual domain arrays or warn if they are specified
    !    but won't be used.
    if (.not. lname(5)) then
      write (this%iout, '(1x,a)') 'Warning.  Dual domain is active but &
        &initial immobile domain concentration was not specified.  &
        &Setting CIM to zero.'
    end if
    if (.not. lname(6)) then
      write (errmsg, '(a)') 'Dual domain is active but dual &
        &domain mass transfer rate (ZETAIM) was not specified.  ZETAIM &
        &must be specified in GRIDDATA block.'
      call store_error(errmsg)
    end if
    if (.not. lname(7)) then
      write (errmsg, '(a)') 'Dual domain is active but &
        &immobile domain POROSITY was not specified.  POROSITY &
        &must be specified in GRIDDATA block.'
      call store_error(errmsg)
    end if
    if (.not. lname(8)) then
      write (errmsg, '(a)') 'Dual domain is active but &
        &immobile domain VOLFRAC was not specified.  VOLFRAC &
        &must be specified in GRIDDATA block. This is a new &
        &requirement for MODFLOW versions later than version &
        &6.4.1.'
      call store_error(errmsg)
    end if
    !
    ! -- terminate if errors
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if

Here is the call graph for this function:

read_options()

subroutine gwtistmodule::read_options

(

class(gwtisttype), intent(inout)

this

)

Read options for boundary packages.

Parameters

[in,out]

this

GwtIstType object

Definition at line 973 of file gwt-ist.f90.

    ! -- modules
    use constantsmodule, only: linelength, mnormal
    use simmodule, only: store_error
    use openspecmodule, only: access, form
    use inputoutputmodule, only: getunit, assign_iounit, openfile
    ! -- dummy
    class(GwtIstType), intent(inout) :: this !< GwtIstType object
    ! -- local
    character(len=LINELENGTH) :: errmsg, keyword
    character(len=LINELENGTH) :: sorption
    character(len=LINELENGTH) :: fname
    character(len=:), allocatable :: keyword2
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    logical :: found
    ! -- 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 SELECTED. ')"
    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 :: fmtistbin = &
      "(4x, 'IST ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, &
      &/4x, 'OPENED ON UNIT: ', I0)"
    !
    ! -- 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 IMMOBILE 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 ('CIM')
          call this%parser%GetRemainingLine(keyword2)
          call this%ocd%set_option(keyword2, this%inunit, this%iout)
        case ('BUDGET')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            call assign_iounit(this%ibudgetout, this%inunit, "BUDGET fileout")
            call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
                          form, access, 'REPLACE', mode_opt=mnormal)
            write (this%iout, fmtistbin) 'BUDGET', trim(adjustl(fname)), &
              this%ibudgetout
            found = .true.
          else
            call store_error('Optional BUDGET keyword must &
                             &be followed by FILEOUT')
          end if
        case ('BUDGETCSV')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            call assign_iounit(this%ibudcsv, this%inunit, "BUDGETCSV fileout")
            call openfile(this%ibudcsv, this%iout, fname, 'CSV', &
                          filstat_opt='REPLACE')
            write (this%iout, fmtistbin) 'BUDGET CSV', trim(adjustl(fname)), &
              this%ibudcsv
          else
            call store_error('Optional BUDGETCSV keyword must be followed by &
              &FILEOUT')
          end if
        case ('SORBTION')
          call store_error('SORBTION is not a valid option.  Use &
                           &SORPTION instead.')
          call this%parser%StoreErrorUnit()
        case ('SORPTION')
          call this%parser%GetStringCaps(sorption)
          select case (sorption)
          case ('LINEAR', '')
            this%isrb = 1
            write (this%iout, fmtlinear)
          case ('FREUNDLICH')
            this%isrb = 2
            write (this%iout, fmtfreundlich)
          case ('LANGMUIR')
            this%isrb = 3
            write (this%iout, fmtlangmuir)
          case default
            call store_error('Unknown sorption type was specified ' &
                             & //'('//trim(adjustl(sorption))//').'// &
                             &' Sorption must be specified as LINEAR, &
                             &FREUNDLICH, or LANGMUIR.')
            call this%parser%StoreErrorUnit()
          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 ('SORBATE')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            this%ioutsorbate = getunit()
            call openfile(this%ioutsorbate, this%iout, fname, 'DATA(BINARY)', &
                          form, access, 'REPLACE', mode_opt=mnormal)
            write (this%iout, fmtistbin) &
              'SORBATE', fname, this%ioutsorbate
          else
            errmsg = 'Optional SORBATE keyword must be '// &
                     'followed by FILEOUT.'
            call store_error(errmsg)
          end if
        case default
          write (errmsg, '(a,a)') 'Unknown IST option: ', &
            trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END OF IMMOBILE STORAGE AND TRANSFER &
                                &OPTIONS'
    end if

Here is the call graph for this function:

Variable Documentation

budtxt

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

private

Definition at line 38 of file gwt-ist.f90.

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

ftype

character(len=lenftype) gwtistmodule::ftype = 'IST'

private

Definition at line 35 of file gwt-ist.f90.

  character(len=LENFTYPE) :: ftype = 'IST'

nbditems

integer(i4b), parameter gwtistmodule::nbditems = 5

private

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

  integer(I4B), parameter :: NBDITEMS = 5

text

character(len=lenpackagename) gwtistmodule::text = ' IMMOBILE DOMAIN'

private

Definition at line 36 of file gwt-ist.f90.

  character(len=LENPACKAGENAME) :: text = ' IMMOBILE DOMAIN'