MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
Data Types | Enumerations | Functions/Subroutines | Variables
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 , sorption_off = 0 ,
  sorption_linear = 1 , sorption_freund = 2 , sorption_lang = 3
}
 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_srb_term (isrb, volfracm, rhobm, vcell, tled, cnew, cold, swnew, swold, const1, const2, rate, hcofval, rhsval)
 @ brief Calculate sorption terms More...
 
subroutine mst_fc_dcy_srb (this, nodes, cold, nja, matrix_sln, idxglo, rhs, cnew, kiter)
 @ brief Fill sorption-decay coefficient method for package More...
 
subroutine mst_cq (this, nodes, cnew, cold, flowja)
 @ brief Calculate flows for package More...
 
subroutine mst_cq_sto (this, nodes, cnew, cold, flowja)
 @ brief Calculate storage terms for package More...
 
subroutine mst_cq_dcy (this, nodes, cnew, cold, flowja)
 @ brief Calculate decay terms for package More...
 
subroutine mst_cq_srb (this, nodes, cnew, cold, flowja)
 @ brief Calculate sorption terms for package More...
 
subroutine mst_cq_dcy_srb (this, nodes, cnew, cold, flowja)
 @ brief Calculate decay-sorption terms for package More...
 
subroutine mst_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_freundlich_conc (conc, kf, a)
 @ brief Calculate sorption concentration using Freundlich More...
 
real(dp) function get_langmuir_conc (conc, kl, sbar)
 @ brief Calculate sorption concentration using Langmuir More...
 
real(dp) function get_freundlich_derivative (conc, kf, a)
 @ brief Calculate sorption derivative using Freundlich More...
 
real(dp) function get_langmuir_derivative (conc, kl, sbar)
 @ brief Calculate sorption derivative using Langmuir More...
 
real(dp) function get_freundlich_kd (conc, kf, a)
 @ brief Get effective Freundlich distribution coefficient More...
 
real(dp) function get_langmuir_kd (conc, kl, sbar)
 @ brief Get effective Langmuir distribution coefficient 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.

sorption_off 

Sorption is inactive (default)

sorption_linear 

Linear sorption between aqueous and solid phases.

sorption_freund 

Freundlich sorption between aqueous and solid phases.

sorption_lang 

Langmuir sorption between aqueous and solid phases.

Definition at line 34 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
thisGwtMstType object
[in]volfracimimmobile domain volume fraction that contributes to total immobile volume fraction

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

1480  ! -- dummy
1481  class(GwtMstType) :: this !< GwtMstType object
1482  real(DP), dimension(:), intent(in) :: volfracim !< immobile domain volume fraction that contributes to total immobile volume fraction
1483  ! -- local
1484  integer(I4B) :: n
1485  !
1486  ! -- Add to volfracim
1487  do n = 1, this%dis%nodes
1488  this%volfracim(n) = this%volfracim(n) + volfracim(n)
1489  end do
1490  !
1491  ! -- An immobile domain is adding a volume fraction, so update thetam
1492  ! accordingly.
1493  do n = 1, this%dis%nodes
1494  this%thetam(n) = this%get_volfracm(n) * this%porosity(n)
1495  end do

◆ allocate_arrays()

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

Method to allocate arrays for the package.

Parameters
thisGwtMstType object
[in]nodesnumber of nodes

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

1087  ! -- modules
1088  use memorymanagermodule, only: mem_allocate
1089  use constantsmodule, only: dzero
1090  ! -- dummy
1091  class(GwtMstType) :: this !< GwtMstType object
1092  integer(I4B), intent(in) :: nodes !< number of nodes
1093  ! -- local
1094  integer(I4B) :: n
1095  !
1096  ! -- Allocate
1097  ! -- sto
1098  call mem_allocate(this%porosity, nodes, 'POROSITY', this%memoryPath)
1099  call mem_allocate(this%thetam, nodes, 'THETAM', this%memoryPath)
1100  call mem_allocate(this%volfracim, nodes, 'VOLFRACIM', this%memoryPath)
1101  call mem_allocate(this%ratesto, nodes, 'RATESTO', this%memoryPath)
1102  !
1103  ! -- dcy
1104  if (this%idcy == decay_off) then
1105  call mem_allocate(this%ratedcy, 1, 'RATEDCY', this%memoryPath)
1106  call mem_allocate(this%decay, 1, 'DECAY', this%memoryPath)
1107  call mem_allocate(this%decaylast, 1, 'DECAYLAST', this%memoryPath)
1108  else
1109  call mem_allocate(this%ratedcy, this%dis%nodes, 'RATEDCY', this%memoryPath)
1110  call mem_allocate(this%decay, nodes, 'DECAY', this%memoryPath)
1111  call mem_allocate(this%decaylast, nodes, 'DECAYLAST', this%memoryPath)
1112  end if
1113  if (this%idcy /= decay_off .and. this%isrb /= sorption_off) then
1114  call mem_allocate(this%ratedcys, this%dis%nodes, 'RATEDCYS', &
1115  this%memoryPath)
1116  call mem_allocate(this%decayslast, this%dis%nodes, 'DECAYSLAST', &
1117  this%memoryPath)
1118  else
1119  call mem_allocate(this%ratedcys, 1, 'RATEDCYS', this%memoryPath)
1120  call mem_allocate(this%decayslast, 1, 'DECAYSLAST', this%memoryPath)
1121  end if
1122  call mem_allocate(this%decay_sorbed, 1, 'DECAY_SORBED', &
1123  this%memoryPath)
1124  !
1125  ! -- srb
1126  if (this%isrb == sorption_off) then
1127  call mem_allocate(this%bulk_density, 1, 'BULK_DENSITY', this%memoryPath)
1128  call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
1129  call mem_allocate(this%distcoef, 1, 'DISTCOEF', this%memoryPath)
1130  call mem_allocate(this%ratesrb, 1, 'RATESRB', this%memoryPath)
1131  call mem_allocate(this%csrb, 1, 'CSRB', this%memoryPath)
1132  else
1133  call mem_allocate(this%bulk_density, nodes, 'BULK_DENSITY', this%memoryPath)
1134  call mem_allocate(this%distcoef, nodes, 'DISTCOEF', this%memoryPath)
1135  call mem_allocate(this%ratesrb, nodes, 'RATESRB', this%memoryPath)
1136  call mem_allocate(this%csrb, nodes, 'CSRB', this%memoryPath)
1137  if (this%isrb == sorption_linear) then
1138  call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
1139  else
1140  call mem_allocate(this%sp2, nodes, 'SP2', this%memoryPath)
1141  end if
1142  end if
1143  !
1144  ! -- Initialize
1145  do n = 1, nodes
1146  this%porosity(n) = dzero
1147  this%thetam(n) = dzero
1148  this%volfracim(n) = dzero
1149  this%ratesto(n) = dzero
1150  end do
1151  do n = 1, size(this%decay)
1152  this%decay(n) = dzero
1153  this%ratedcy(n) = dzero
1154  this%decaylast(n) = dzero
1155  end do
1156  do n = 1, size(this%bulk_density)
1157  this%bulk_density(n) = dzero
1158  this%distcoef(n) = dzero
1159  this%ratesrb(n) = dzero
1160  this%csrb(n) = dzero
1161  end do
1162  do n = 1, size(this%sp2)
1163  this%sp2(n) = dzero
1164  end do
1165  do n = 1, size(this%ratedcys)
1166  this%ratedcys(n) = dzero
1167  this%decayslast(n) = dzero
1168  end do
memorymanagermodule::mem_allocate
Definition: MemoryManager.f90:54
constantsmodule
This module contains simulation constants.
Definition: Constants.f90:9
constantsmodule::dzero
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65

◆ allocate_scalars()

subroutine gwtmstmodule::allocate_scalars ( class(gwtmsttype)  this)

Method to allocate scalar variables for the package.

Parameters
thisGwtMstType object

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

1063  ! -- modules
1064  use memorymanagermodule, only: mem_allocate, mem_setptr
1065  ! -- dummy
1066  class(GwtMstType) :: this !< GwtMstType object
1067  !
1068  ! -- Allocate scalars in NumericalPackageType
1069  call this%NumericalPackageType%allocate_scalars()
1070  !
1071  ! -- Allocate
1072  call mem_allocate(this%isrb, 'ISRB', this%memoryPath)
1073  call mem_allocate(this%ioutsorbate, 'IOUTSORBATE', this%memoryPath)
1074  call mem_allocate(this%idcy, 'IDCY', this%memoryPath)
1075  !
1076  ! -- Initialize
1077  this%isrb = izero
1078  this%ioutsorbate = 0
1079  this%idcy = izero
memorymanagermodule::mem_setptr
Definition: MemoryManager.f90:89

◆ get_freundlich_conc()

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

Function to calculate sorption concentration using Freundlich

Parameters
[in]concsolute concentration
[in]kffreundlich constant
[in]afreundlich exponent

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

1517  ! -- dummy
1518  real(DP), intent(in) :: conc !< solute concentration
1519  real(DP), intent(in) :: kf !< freundlich constant
1520  real(DP), intent(in) :: a !< freundlich exponent
1521  ! -- return
1522  real(DP) :: cbar
1523  !
1524  if (conc > dzero) then
1525  cbar = kf * conc**a
1526  else
1527  cbar = dzero
1528  end if
Here is the caller graph for this function:

◆ get_freundlich_derivative()

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

Function to calculate sorption derivative using Freundlich

Parameters
[in]concsolute concentration
[in]kffreundlich constant
[in]afreundlich exponent

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

1555  ! -- dummy
1556  real(DP), intent(in) :: conc !< solute concentration
1557  real(DP), intent(in) :: kf !< freundlich constant
1558  real(DP), intent(in) :: a !< freundlich exponent
1559  ! -- return
1560  real(DP) :: derv
1561  !
1562  if (conc > dzero) then
1563  derv = kf * a * conc**(a - done)
1564  else
1565  derv = dzero
1566  end if
Here is the caller graph for this function:

◆ get_freundlich_kd()

real(dp) function gwtmstmodule::get_freundlich_kd ( real(dp), intent(in)  conc,
real(dp), intent(in)  kf,
real(dp), intent(in)  a 
)
Parameters
[in]concsolute concentration
[in]kffreundlich constant
[in]afreundlich exponent
Returns
effective distribution coefficient

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

1591  ! -- dummy
1592  real(DP), intent(in) :: conc !< solute concentration
1593  real(DP), intent(in) :: kf !< freundlich constant
1594  real(DP), intent(in) :: a !< freundlich exponent
1595  ! -- return
1596  real(DP) :: kd !< effective distribution coefficient
1597  !
1598  if (conc > dzero) then
1599  kd = kf * conc**(a - done)
1600  else
1601  kd = dzero
1602  end if
Here is the caller graph for this function:

◆ get_langmuir_conc()

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

Function to calculate sorption concentration using Langmuir

Parameters
[in]concsolute concentration
[in]kllangmuir constant
[in]sbarlangmuir sorption sites

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

1536  ! -- dummy
1537  real(DP), intent(in) :: conc !< solute concentration
1538  real(DP), intent(in) :: kl !< langmuir constant
1539  real(DP), intent(in) :: sbar !< langmuir sorption sites
1540  ! -- return
1541  real(DP) :: cbar
1542  !
1543  if (conc > dzero) then
1544  cbar = (kl * sbar * conc) / (done + kl * conc)
1545  else
1546  cbar = dzero
1547  end if
Here is the caller graph for this function:

◆ get_langmuir_derivative()

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

Function to calculate sorption derivative using Langmuir

Parameters
[in]concsolute concentration
[in]kllangmuir constant
[in]sbarlangmuir sorption sites

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

1574  ! -- dummy
1575  real(DP), intent(in) :: conc !< solute concentration
1576  real(DP), intent(in) :: kl !< langmuir constant
1577  real(DP), intent(in) :: sbar !< langmuir sorption sites
1578  ! -- return
1579  real(DP) :: derv
1580  !
1581  if (conc > dzero) then
1582  derv = (kl * sbar) / (done + kl * conc)**dtwo
1583  else
1584  derv = dzero
1585  end if
Here is the caller graph for this function:

◆ get_langmuir_kd()

real(dp) function gwtmstmodule::get_langmuir_kd ( real(dp), intent(in)  conc,
real(dp), intent(in)  kl,
real(dp), intent(in)  sbar 
)
Parameters
[in]concsolute concentration
[in]kllangmuir constant
[in]sbarlangmuir sorption sites
Returns
effective distribution coefficient

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

1608  ! -- dummy
1609  real(DP), intent(in) :: conc !< solute concentration
1610  real(DP), intent(in) :: kl !< langmuir constant
1611  real(DP), intent(in) :: sbar !< langmuir sorption sites
1612  ! -- return
1613  real(DP) :: kd !< effective distribution coefficient
1614  !
1615  if (conc > dzero) then
1616  kd = (kl * sbar) / (done + kl * conc)
1617  else
1618  kd = dzero
1619  end if
Here is the caller graph for this function:

◆ get_volfracm()

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

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

Parameters
thisGwtMstType object
[in]nodenode number

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

1503  ! -- dummy
1504  class(GwtMstType) :: this !< GwtMstType object
1505  integer(I4B), intent(in) :: node !< node number
1506  ! -- return
1507  real(DP) :: volfracm
1508  !
1509  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_usruser-entered decay rate
[in]decay_rate_lastdecay rate used for last iteration
[in]kiterPicard iteration counter
[in]coldconcentration at end of last time step
[in]cnewconcentration at end of this time step
[in]deltlength of time step
Returns
returned value for decay rate

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

1632  ! -- dummy
1633  real(DP), intent(in) :: decay_rate_usr !< user-entered decay rate
1634  real(DP), intent(in) :: decay_rate_last !< decay rate used for last iteration
1635  integer(I4B), intent(in) :: kiter !< Picard iteration counter
1636  real(DP), intent(in) :: cold !< concentration at end of last time step
1637  real(DP), intent(in) :: cnew !< concentration at end of this time step
1638  real(DP), intent(in) :: delt !< length of time step
1639  ! -- return
1640  real(DP) :: decay_rate !< returned value for decay rate
1641  !
1642  ! -- Return user rate if production, otherwise constrain, if necessary
1643  if (decay_rate_usr < dzero) then
1644  !
1645  ! -- Production, no need to limit rate
1646  decay_rate = decay_rate_usr
1647  else
1648  !
1649  ! -- Need to ensure decay does not result in negative
1650  ! concentration, so reduce the rate if it would result in
1651  ! removing more mass than is in the cell.
1652  if (kiter == 1) then
1653  decay_rate = min(decay_rate_usr, cold / delt)
1654  else
1655  decay_rate = decay_rate_last
1656  if (cnew < dzero) then
1657  decay_rate = decay_rate_last + cnew / delt
1658  else if (cnew > cold) then
1659  decay_rate = decay_rate_last + cnew / delt
1660  end if
1661  decay_rate = min(decay_rate_usr, decay_rate)
1662  end if
1663  decay_rate = max(decay_rate, dzero)
1664  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 1447 of file gwt-mst.f90.

1448  use gwtmstinputmodule, only: gwtmstparamfoundtype
1449  class(GwTMstType) :: this
1450  type(GwtMstParamFoundType), intent(in) :: found
1451 
1452  write (this%iout, '(1x,a)') 'PROCESSING GRIDDATA'
1453  if (found%porosity) then
1454  write (this%iout, '(4x,a)') 'MOBILE DOMAIN POROSITY set from input file'
1455  end if
1456  if (found%decay) then
1457  write (this%iout, '(4x,a)') 'DECAY RATE set from input file'
1458  end if
1459  if (found%decay_sorbed) then
1460  write (this%iout, '(4x,a)') 'DECAY SORBED RATE set from input file'
1461  end if
1462  if (found%bulk_density) then
1463  write (this%iout, '(4x,a)') 'BULK DENSITY set from input file'
1464  end if
1465  if (found%distcoef) then
1466  write (this%iout, '(4x,a)') 'DISTRIBUTION COEFFICIENT set from input file'
1467  end if
1468  if (found%sp2) then
1469  write (this%iout, '(4x,a)') 'SECOND SORPTION PARAM set from input file'
1470  end if
1471  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 1228 of file gwt-mst.f90.

1229  use gwtmstinputmodule, only: gwtmstparamfoundtype
1230  class(GwTMstType) :: this
1231  type(GwtMstParamFoundType), intent(in) :: found
1232  character(len=*), intent(in) :: sorbate_fname
1233  ! -- formats
1234  character(len=*), parameter :: fmtisvflow = &
1235  "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
1236  &WHENEVER ICBCFL IS NOT ZERO.')"
1237  character(len=*), parameter :: fmtlinear = &
1238  &"(4x,'LINEAR SORPTION IS ACTIVE. ')"
1239  character(len=*), parameter :: fmtfreundlich = &
1240  &"(4x,'FREUNDLICH SORPTION IS ACTIVE. ')"
1241  character(len=*), parameter :: fmtlangmuir = &
1242  &"(4x,'LANGMUIR SORPTION IS ACTIVE. ')"
1243  character(len=*), parameter :: fmtidcy1 = &
1244  &"(4x,'FIRST-ORDER DECAY IS ACTIVE. ')"
1245  character(len=*), parameter :: fmtidcy2 = &
1246  &"(4x,'ZERO-ORDER DECAY IS ACTIVE. ')"
1247  character(len=*), parameter :: fmtfileout = &
1248  "(4x,'MST ',1x,a,1x,' WILL BE SAVED TO FILE: ',a,/4x,&
1249  &'OPENED ON UNIT: ',I7)"
1250 
1251  write (this%iout, '(1x,a)') 'PROCESSING MOBILE STORAGE AND TRANSFER OPTIONS'
1252  if (found%save_flows) then
1253  write (this%iout, fmtisvflow)
1254  end if
1255  if (found%order1_decay) then
1256  write (this%iout, fmtidcy1)
1257  end if
1258  if (found%order0_decay) then
1259  write (this%iout, fmtidcy2)
1260  end if
1261  if (found%sorption) then
1262  select case (this%isrb)
1263  case (sorption_linear)
1264  write (this%iout, fmtlinear)
1265  case (sorption_freund)
1266  write (this%iout, fmtfreundlich)
1267  case (sorption_lang)
1268  write (this%iout, fmtlangmuir)
1269  end select
1270  end if
1271  if (found%sorbatefile) then
1272  write (this%iout, fmtfileout) &
1273  'SORBATE', sorbate_fname, this%ioutsorbate
1274  end if
1275  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]thisGwtMstType object
[in]dispointer to dis package
iboundpointer to GWT ibound array

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

144  ! -- modules
145  ! -- dummy
146  class(GwtMstType), intent(inout) :: this !< GwtMstType object
147  class(DisBaseType), pointer, intent(in) :: dis !< pointer to dis package
148  integer(I4B), dimension(:), pointer, contiguous :: ibound !< pointer to GWT ibound array
149  ! -- local
150  ! -- formats
151  character(len=*), parameter :: fmtmst = &
152  "(1x,/1x,'MST -- MOBILE STORAGE AND TRANSFER PACKAGE, VERSION 1, &
153  &7/29/2020 INPUT READ FROM MEMPATH: ', A, //)"
154  !
155  ! --print a message identifying the mobile storage and transfer package.
156  write (this%iout, fmtmst) this%input_mempath
157  !
158  ! -- Source options
159  call this%source_options()
160  !
161  ! -- store pointers to arguments that were passed in
162  this%dis => dis
163  this%ibound => ibound
164  !
165  ! -- Allocate arrays
166  call this%allocate_arrays(dis%nodes)
167  !
168  ! -- source the data block
169  call this%source_data()

◆ mst_bd()

subroutine gwtmstmodule::mst_bd ( class(gwtmsttype)  this,
integer(i4b), intent(in)  isuppress_output,
type(budgettype), intent(inout)  model_budget 
)
Parameters
thisGwtMstType object
[in]isuppress_outputflag to suppress output
[in,out]model_budgetmodel budget object

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

889  ! -- modules
890  use tdismodule, only: delt
891  use budgetmodule, only: budgettype, rate_accumulator
892  ! -- dummy
893  class(GwtMstType) :: this !< GwtMstType object
894  integer(I4B), intent(in) :: isuppress_output !< flag to suppress output
895  type(BudgetType), intent(inout) :: model_budget !< model budget object
896  ! -- local
897  real(DP) :: rin
898  real(DP) :: rout
899  !
900  ! -- sto
901  call rate_accumulator(this%ratesto, rin, rout)
902  call model_budget%addentry(rin, rout, delt, budtxt(1), &
903  isuppress_output, rowlabel=this%packName)
904  !
905  ! -- dcy
906  if (this%idcy /= decay_off) then
907  call rate_accumulator(this%ratedcy, rin, rout)
908  call model_budget%addentry(rin, rout, delt, budtxt(2), &
909  isuppress_output, rowlabel=this%packName)
910  end if
911  !
912  ! -- srb
913  if (this%isrb /= sorption_off) then
914  call rate_accumulator(this%ratesrb, rin, rout)
915  call model_budget%addentry(rin, rout, delt, budtxt(3), &
916  isuppress_output, rowlabel=this%packName)
917  end if
918  !
919  ! -- srb dcy
920  if (this%isrb /= sorption_off .and. this%idcy /= decay_off) then
921  call rate_accumulator(this%ratedcys, rin, rout)
922  call model_budget%addentry(rin, rout, delt, budtxt(4), &
923  isuppress_output, rowlabel=this%packName)
924  end if
budgetmodule
This module contains the BudgetModule.
Definition: Budget.f90:20
budgetmodule::rate_accumulator
subroutine, public rate_accumulator(flow, rin, rout)
@ brief Rate accumulator subroutine
Definition: Budget.f90:632
tdismodule::delt
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
budgetmodule::budgettype
Derived type for the Budget object.
Definition: Budget.f90:39
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
thisGwtMstType object
[in]cnewconcentration at end of this time step

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

859  ! -- dummy
860  class(GwtMstType) :: this !< GwtMstType object
861  real(DP), intent(in), dimension(:) :: cnew !< concentration at end of this time step
862  ! -- local
863  integer(I4B) :: n
864  real(DP) :: distcoef
865  real(DP) :: csrb
866 
867  ! Calculate sorbed concentration
868  do n = 1, size(cnew)
869  csrb = dzero
870  if (this%ibound(n) > 0) then
871  distcoef = this%distcoef(n)
872  select case (this%isrb)
873  case (sorption_linear)
874  csrb = cnew(n) * distcoef
875  case (sorption_freund)
876  csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
877  case (sorption_lang)
878  csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
879  end select
880  end if
881  this%csrb(n) = csrb
882  end do
883 
Here is the call graph for this function:

◆ mst_cq()

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

Method to calculate flows for the package.

Parameters
thisGwtMstType object
[in]nodesnumber of nodes
[in]cnewconcentration at end of this time step
[in]coldconcentration at end of last time step
[in,out]flowjaflow between two connected control volumes

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

558  ! -- dummy
559  class(GwtMstType) :: this !< GwtMstType object
560  integer(I4B), intent(in) :: nodes !< number of nodes
561  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
562  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
563  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
564  !
565  ! - storage
566  call this%mst_cq_sto(nodes, cnew, cold, flowja)
567  !
568  ! -- decay
569  if (this%idcy /= decay_off) then
570  call this%mst_cq_dcy(nodes, cnew, cold, flowja)
571  end if
572  !
573  ! -- sorption
574  if (this%isrb /= sorption_off) then
575  call this%mst_cq_srb(nodes, cnew, cold, flowja)
576  end if
577  !
578  ! -- decay sorbed
579  if (this%isrb /= sorption_off .and. this%idcy /= decay_off) then
580  call this%mst_cq_dcy_srb(nodes, cnew, cold, flowja)
581  end if
582  !
583  ! -- calculate csrb
584  if (this%isrb /= sorption_off) then
585  call this%mst_calc_csrb(cnew)
586  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
thisGwtMstType object
[in]nodesnumber of nodes
[in]cnewconcentration at end of this time step
[in]coldconcentration at end of last time step
[in,out]flowjaflow between two connected control volumes

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

642  ! -- modules
643  use tdismodule, only: delt
644  ! -- dummy
645  class(GwtMstType) :: this !< GwtMstType object
646  integer(I4B), intent(in) :: nodes !< number of nodes
647  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
648  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
649  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
650  ! -- local
651  integer(I4B) :: n
652  integer(I4B) :: idiag
653  real(DP) :: rate
654  real(DP) :: swtpdt
655  real(DP) :: hhcof, rrhs
656  real(DP) :: vcell
657  real(DP) :: decay_rate
658  !
659  ! -- initialize
660  !
661  ! -- Calculate decay change
662  do n = 1, nodes
663  !
664  ! -- skip if transport inactive
665  this%ratedcy(n) = dzero
666  if (this%ibound(n) <= 0) cycle
667  !
668  ! -- calculate new and old water volumes
669  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
670  swtpdt = this%fmi%gwfsat(n)
671  !
672  ! -- calculate decay gains and losses
673  rate = dzero
674  hhcof = dzero
675  rrhs = dzero
676  if (this%idcy == decay_first_order) then
677  if (cnew(n) > dzero) then
678  hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
679  end if
680  elseif (this%idcy == decay_zero_order) then
681  decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), &
682  0, cold(n), cnew(n), delt)
683  rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
684  end if
685  rate = hhcof * cnew(n) - rrhs
686  this%ratedcy(n) = rate
687  idiag = this%dis%con%ia(n)
688  flowja(idiag) = flowja(idiag) + rate
689  !
690  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
thisGwtMstType object
[in]nodesnumber of nodes
[in]cnewconcentration at end of this time step
[in]coldconcentration at end of last time step
[in,out]flowjaflow between two connected control volumes

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

756  ! -- modules
757  use tdismodule, only: delt
758  ! -- dummy
759  class(GwtMstType) :: this !< GwtMstType object
760  integer(I4B), intent(in) :: nodes !< number of nodes
761  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
762  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
763  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
764  ! -- local
765  integer(I4B) :: n
766  integer(I4B) :: idiag
767  real(DP) :: rate
768  real(DP) :: hhcof, rrhs
769  real(DP) :: vcell
770  real(DP) :: swnew
771  real(DP) :: distcoef
772  real(DP) :: volfracm
773  real(DP) :: rhobm
774  real(DP) :: term
775  real(DP) :: csrb
776  real(DP) :: csrbnew
777  real(DP) :: csrbold
778  real(DP) :: decay_rate
779  !
780  ! -- Calculate sorbed decay change
781  ! This routine will only be called if sorption and decay are active
782  do n = 1, nodes
783  !
784  ! -- initialize rates
785  this%ratedcys(n) = dzero
786  !
787  ! -- skip if transport inactive
788  if (this%ibound(n) <= 0) cycle
789  !
790  ! -- set variables
791  hhcof = dzero
792  rrhs = dzero
793  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
794  swnew = this%fmi%gwfsat(n)
795  distcoef = this%distcoef(n)
796  volfracm = this%get_volfracm(n)
797  rhobm = this%bulk_density(n)
798  term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
799  !
800  ! -- add sorbed mass decay rate terms to accumulators
801  select case (this%idcy)
802  case (decay_first_order)
803 
804  !
805  select case (this%isrb)
806  case (sorption_linear)
807  !
808  ! -- first order decay rate is a function of concentration, so add
809  ! to left hand side
810  if (cnew(n) > dzero) then
811  hhcof = -term * distcoef
812  end if
813  case (sorption_freund)
814  !
815  ! -- nonlinear Freundlich sorption, so add to RHS
816  csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
817  rrhs = term * csrb
818  case (sorption_lang)
819  !
820  ! -- nonlinear Lanmuir sorption, so add to RHS
821  csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
822  rrhs = term * csrb
823  end select
824  case (decay_zero_order)
825  !
826  ! -- Call function to get zero-order decay rate, which may be changed
827  ! from the user-specified rate to prevent negative concentrations
828  if (distcoef > dzero) then
829  select case (this%isrb)
830  case (sorption_linear)
831  csrbold = cold(n) * distcoef
832  csrbnew = cnew(n) * distcoef
833  case (sorption_freund)
834  csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n))
835  csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
836  case (sorption_lang)
837  csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n))
838  csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
839  end select
840  decay_rate = get_zero_order_decay(this%decay_sorbed(n), &
841  this%decayslast(n), &
842  0, csrbold, csrbnew, delt)
843  rrhs = decay_rate * volfracm * rhobm * swnew * vcell
844  end if
845  end select
846  !
847  ! -- calculate rate
848  rate = hhcof * cnew(n) - rrhs
849  this%ratedcys(n) = rate
850  idiag = this%dis%con%ia(n)
851  flowja(idiag) = flowja(idiag) + rate
852  !
853  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
thisGwtMstType object
[in]nodesnumber of nodes
[in]cnewconcentration at end of this time step
[in]coldconcentration at end of last time step
[in,out]flowjaflow between two connected control volumes

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

698  ! -- modules
699  use tdismodule, only: delt
700  ! -- dummy
701  class(GwtMstType) :: this !< GwtMstType object
702  integer(I4B), intent(in) :: nodes !< number of nodes
703  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
704  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
705  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
706  ! -- local
707  integer(I4B) :: n
708  integer(I4B) :: idiag
709  real(DP) :: rate
710  real(DP) :: tled
711  real(DP) :: swt, swtpdt
712  real(DP) :: vcell
713  real(DP) :: volfracm
714  real(DP) :: rhobm
715  real(DP) :: const1
716  real(DP) :: const2
717  !
718  ! -- initialize
719  tled = done / delt
720  !
721  ! -- Calculate sorption change
722  do n = 1, nodes
723  !
724  ! -- initialize rates
725  this%ratesrb(n) = dzero
726  !
727  ! -- skip if transport inactive
728  if (this%ibound(n) <= 0) cycle
729  !
730  ! -- assign variables
731  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
732  swtpdt = this%fmi%gwfsat(n)
733  swt = this%fmi%gwfsatold(n, delt)
734  volfracm = this%get_volfracm(n)
735  rhobm = this%bulk_density(n)
736  const1 = this%distcoef(n)
737  const2 = 0.
738  if (this%isrb == sorption_freund .or. this%isrb == sorption_lang) then
739  const2 = this%sp2(n)
740  end if
741  call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
742  cold(n), swtpdt, swt, const1, const2, &
743  rate=rate)
744  this%ratesrb(n) = rate
745  idiag = this%dis%con%ia(n)
746  flowja(idiag) = flowja(idiag) + rate
747  !
748  end do
Here is the call graph for this function:

◆ mst_cq_sto()

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

Method to calculate storage terms for the package.

Parameters
thisGwtMstType object
[in]nodesnumber of nodes
[in]cnewconcentration at end of this time step
[in]coldconcentration at end of last time step
[in,out]flowjaflow between two connected control volumes

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

594  ! -- modules
595  use tdismodule, only: delt
596  ! -- dummy
597  class(GwtMstType) :: this !< GwtMstType object
598  integer(I4B), intent(in) :: nodes !< number of nodes
599  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
600  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
601  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
602  ! -- local
603  integer(I4B) :: n
604  integer(I4B) :: idiag
605  real(DP) :: rate
606  real(DP) :: tled
607  real(DP) :: vnew, vold
608  real(DP) :: hhcof, rrhs
609  !
610  ! -- initialize
611  tled = done / delt
612  !
613  ! -- Calculate storage change
614  do n = 1, nodes
615  this%ratesto(n) = dzero
616  !
617  ! -- skip if transport inactive
618  if (this%ibound(n) <= 0) cycle
619  !
620  ! -- calculate new and old water volumes
621  vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
622  this%fmi%gwfsat(n) * this%thetam(n)
623  vold = vnew
624  if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) * delt
625  if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) * delt
626  !
627  ! -- calculate rate
628  hhcof = -vnew * tled
629  rrhs = -vold * tled * cold(n)
630  rate = hhcof * cnew(n) - rrhs
631  this%ratesto(n) = rate
632  idiag = this%dis%con%ia(n)
633  flowja(idiag) = flowja(idiag) + rate
634  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
mstobjunallocated new mst object to create
[in]name_modelname of the model
[in]input_mempathmemory path of input
[in]inunitunit number of WEL package input file
[in]ioutunit number of model listing file
[in]fmifmi package for this GWT model

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

116  ! -- dummy
117  type(GwtMstType), pointer :: mstobj !< unallocated new mst object to create
118  character(len=*), intent(in) :: name_model !< name of the model
119  character(len=*), intent(in) :: input_mempath !< memory path of input
120  integer(I4B), intent(in) :: inunit !< unit number of WEL package input file
121  integer(I4B), intent(in) :: iout !< unit number of model listing file
122  type(TspFmiType), intent(in), target :: fmi !< fmi package for this GWT model
123  !
124  ! -- Create the object
125  allocate (mstobj)
126  !
127  ! -- create name and memory path
128  call mstobj%set_names(1, name_model, 'MST', 'MST', input_mempath)
129  !
130  ! -- Allocate scalars
131  call mstobj%allocate_scalars()
132  !
133  ! -- Set variables
134  mstobj%inunit = inunit
135  mstobj%iout = iout
136  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
thisGwtMstType object

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

1023  ! -- modules
1024  use memorymanagermodule, only: mem_deallocate
1025  ! -- dummy
1026  class(GwtMstType) :: this !< GwtMstType object
1027  !
1028  ! -- Deallocate arrays if package was active
1029  if (this%inunit > 0) then
1030  call mem_deallocate(this%porosity)
1031  call mem_deallocate(this%thetam)
1032  call mem_deallocate(this%volfracim)
1033  call mem_deallocate(this%ratesto)
1034  call mem_deallocate(this%idcy)
1035  call mem_deallocate(this%decay)
1036  call mem_deallocate(this%decay_sorbed)
1037  call mem_deallocate(this%ratedcy)
1038  call mem_deallocate(this%decaylast)
1039  call mem_deallocate(this%decayslast)
1040  call mem_deallocate(this%isrb)
1041  call mem_deallocate(this%ioutsorbate)
1042  call mem_deallocate(this%bulk_density)
1043  call mem_deallocate(this%distcoef)
1044  call mem_deallocate(this%sp2)
1045  call mem_deallocate(this%ratesrb)
1046  call mem_deallocate(this%csrb)
1047  call mem_deallocate(this%ratedcys)
1048  this%ibound => null()
1049  this%fmi => null()
1050  end if
1051  !
1052  ! -- Scalars
1053  !
1054  ! -- deallocate parent
1055  call this%NumericalPackageType%da()
memorymanagermodule::mem_deallocate
Definition: MemoryManager.f90:122

◆ 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
thisGwtMstType object
[in]nodesnumber of nodes
[in]coldconcentration at end of last time step
[in]njanumber of GWT connections
matrix_slnsolution matrix
[in]idxglomapping vector for model (local) to solution (global)
[in,out]rhsright-hand side vector for model
[in]cnewconcentration at end of this time step
[in]kitersolution outer iteration number

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

178  ! -- modules
179  ! -- dummy
180  class(GwtMstType) :: this !< GwtMstType object
181  integer, intent(in) :: nodes !< number of nodes
182  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
183  integer(I4B), intent(in) :: nja !< number of GWT connections
184  class(MatrixBaseType), pointer :: matrix_sln !< solution matrix
185  integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
186  real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
187  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
188  integer(I4B), intent(in) :: kiter !< solution outer iteration number
189  !
190  ! -- storage contribution
191  call this%mst_fc_sto(nodes, cold, nja, matrix_sln, idxglo, rhs)
192  !
193  ! -- decay contribution
194  if (this%idcy /= decay_off) then
195  call this%mst_fc_dcy(nodes, cold, cnew, nja, matrix_sln, idxglo, &
196  rhs, kiter)
197  end if
198  !
199  ! -- sorption contribution
200  if (this%isrb /= sorption_off) then
201  call this%mst_fc_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
202  end if
203  !
204  ! -- decay sorbed contribution
205  if (this%isrb /= sorption_off .and. this%idcy /= decay_off) then
206  call this%mst_fc_dcy_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, &
207  cnew, kiter)
208  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
thisGwtMstType object
[in]nodesnumber of nodes
[in]coldconcentration at end of last time step
[in]cnewconcentration at end of this time step
[in]njanumber of GWT connections
matrix_slnsolution coefficient matrix
[in]idxglomapping vector for model (local) to solution (global)
[in,out]rhsright-hand side vector for model
[in]kitersolution outer iteration number

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

263  ! -- modules
264  use tdismodule, only: delt
265  ! -- dummy
266  class(GwtMstType) :: this !< GwtMstType object
267  integer, intent(in) :: nodes !< number of nodes
268  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
269  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
270  integer(I4B), intent(in) :: nja !< number of GWT connections
271  class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
272  integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
273  real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
274  integer(I4B), intent(in) :: kiter !< solution outer iteration number
275  ! -- local
276  integer(I4B) :: n, idiag
277  real(DP) :: hhcof, rrhs
278  real(DP) :: swtpdt
279  real(DP) :: vcell
280  real(DP) :: decay_rate
281  !
282  ! -- loop through and calculate decay contribution to hcof and rhs
283  do n = 1, this%dis%nodes
284  !
285  ! -- skip if transport inactive
286  if (this%ibound(n) <= 0) cycle
287  !
288  ! -- calculate new and old water volumes
289  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
290  swtpdt = this%fmi%gwfsat(n)
291  !
292  ! -- add decay rate terms to accumulators
293  idiag = this%dis%con%ia(n)
294  select case (this%idcy)
295  case (decay_first_order)
296  !
297  ! -- first order decay rate is a function of concentration, so add
298  ! to left hand side
299  if (cnew(n) > dzero) then
300  hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
301  call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
302  end if
303  case (decay_zero_order)
304  !
305  ! -- Call function to get zero-order decay rate, which may be changed
306  ! from the user-specified rate to prevent negative concentrations
307  decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), &
308  kiter, cold(n), cnew(n), delt)
309  this%decaylast(n) = decay_rate
310  rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
311  rhs(n) = rhs(n) + rrhs
312  end select
313  !
314  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
thisGwtMstType object
[in]nodesnumber of nodes
[in]coldconcentration at end of last time step
[in]njanumber of GWT connections
matrix_slnsolution coefficient matrix
[in]idxglomapping vector for model (local) to solution (global)
[in,out]rhsright-hand side vector for model
[in]cnewconcentration at end of this time step
[in]kitersolution outer iteration number

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

456  ! -- modules
457  use tdismodule, only: delt
458  ! -- dummy
459  class(GwtMstType) :: this !< GwtMstType object
460  integer, intent(in) :: nodes !< number of nodes
461  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
462  integer(I4B), intent(in) :: nja !< number of GWT connections
463  class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
464  integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
465  real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
466  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
467  integer(I4B), intent(in) :: kiter !< solution outer iteration number
468  ! -- local
469  integer(I4B) :: n, idiag
470  real(DP) :: hhcof, rrhs
471  real(DP) :: vcell
472  real(DP) :: swnew
473  real(DP) :: distcoef
474  real(DP) :: volfracm
475  real(DP) :: rhobm
476  real(DP) :: term
477  real(DP) :: csrb
478  real(DP) :: decay_rate
479  real(DP) :: csrbold
480  real(DP) :: csrbnew
481  !
482  ! -- loop through and calculate sorption contribution to hcof and rhs
483  do n = 1, this%dis%nodes
484  !
485  ! -- skip if transport inactive
486  if (this%ibound(n) <= 0) cycle
487  !
488  ! -- set variables
489  hhcof = dzero
490  rrhs = dzero
491  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
492  swnew = this%fmi%gwfsat(n)
493  distcoef = this%distcoef(n)
494  idiag = this%dis%con%ia(n)
495  volfracm = this%get_volfracm(n)
496  rhobm = this%bulk_density(n)
497  term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
498  !
499  ! -- add sorbed mass decay rate terms to accumulators
500  select case (this%idcy)
501  case (decay_first_order)
502  select case (this%isrb)
503  case (sorption_linear)
504  !
505  ! -- first order decay rate is a function of concentration, so add
506  ! to left hand side
507  if (cnew(n) > dzero) then
508  hhcof = -term * distcoef
509  end if
510  case (sorption_freund)
511  !
512  ! -- nonlinear Freundlich sorption, so add to RHS
513  csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
514  rrhs = term * csrb
515  case (sorption_lang)
516  !
517  ! -- nonlinear Lanmuir sorption, so add to RHS
518  csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
519  rrhs = term * csrb
520  end select
521  case (decay_zero_order)
522  !
523  ! -- call function to get zero-order decay rate, which may be changed
524  ! from the user-specified rate to prevent negative concentrations
525  if (distcoef > dzero) then
526  select case (this%isrb)
527  case (sorption_linear)
528  csrbold = cold(n) * distcoef
529  csrbnew = cnew(n) * distcoef
530  case (sorption_freund)
531  csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n))
532  csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
533  case (sorption_lang)
534  csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n))
535  csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
536  end select
537  !
538  decay_rate = get_zero_order_decay(this%decay_sorbed(n), &
539  this%decayslast(n), &
540  kiter, csrbold, csrbnew, delt)
541  this%decayslast(n) = decay_rate
542  rrhs = decay_rate * volfracm * rhobm * swnew * vcell
543  end if
544  end select
545  !
546  ! -- Add hhcof to diagonal and rrhs to right-hand side
547  call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
548  rhs(n) = rhs(n) + rrhs
549  !
550  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
thisGwtMstType object
[in]nodesnumber of nodes
[in]coldconcentration at end of last time step
[in]njanumber of GWT connections
matrix_slnsolution coefficient matrix
[in]idxglomapping vector for model (local) to solution (global)
[in,out]rhsright-hand side vector for model
[in]cnewconcentration at end of this time step

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

323  ! -- modules
324  use tdismodule, only: delt
325  ! -- dummy
326  class(GwtMstType) :: this !< GwtMstType object
327  integer, intent(in) :: nodes !< number of nodes
328  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
329  integer(I4B), intent(in) :: nja !< number of GWT connections
330  class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
331  integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
332  real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
333  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
334  ! -- local
335  integer(I4B) :: n, idiag
336  real(DP) :: tled
337  real(DP) :: hhcof, rrhs
338  real(DP) :: swt, swtpdt
339  real(DP) :: vcell
340  real(DP) :: const1
341  real(DP) :: const2
342  real(DP) :: volfracm
343  real(DP) :: rhobm
344  !
345  ! -- set variables
346  tled = done / delt
347  !
348  ! -- loop through and calculate sorption contribution to hcof and rhs
349  do n = 1, this%dis%nodes
350  !
351  ! -- skip if transport inactive
352  if (this%ibound(n) <= 0) cycle
353  !
354  ! -- assign variables
355  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
356  swtpdt = this%fmi%gwfsat(n)
357  swt = this%fmi%gwfsatold(n, delt)
358  idiag = this%dis%con%ia(n)
359  const1 = this%distcoef(n)
360  const2 = 0.
361  if (this%isrb == sorption_freund .or. this%isrb == sorption_lang) then
362  const2 = this%sp2(n)
363  end if
364  volfracm = this%get_volfracm(n)
365  rhobm = this%bulk_density(n)
366  call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
367  cold(n), swtpdt, swt, const1, const2, &
368  hcofval=hhcof, rhsval=rrhs)
369  !
370  ! -- Add hhcof to diagonal and rrhs to right-hand side
371  call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
372  rhs(n) = rhs(n) + rrhs
373  !
374  end do
Here is the call graph for this function:

◆ mst_fc_sto()

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

Method to calculate and fill storage coefficients for the package.

Parameters
thisGwtMstType object
[in]nodesnumber of nodes
[in]coldconcentration at end of last time step
[in]njanumber of GWT connections
matrix_slnsolution coefficient matrix
[in]idxglomapping vector for model (local) to solution (global)
[in,out]rhsright-hand side vector for model

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

216  ! -- modules
217  use tdismodule, only: delt
218  ! -- dummy
219  class(GwtMstType) :: this !< GwtMstType object
220  integer, intent(in) :: nodes !< number of nodes
221  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
222  integer(I4B), intent(in) :: nja !< number of GWT connections
223  class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix
224  integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
225  real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
226  ! -- local
227  integer(I4B) :: n, idiag
228  real(DP) :: tled
229  real(DP) :: hhcof, rrhs
230  real(DP) :: vnew, vold
231  !
232  ! -- set variables
233  tled = done / delt
234  !
235  ! -- loop through and calculate storage contribution to hcof and rhs
236  do n = 1, this%dis%nodes
237  !
238  ! -- skip if transport inactive
239  if (this%ibound(n) <= 0) cycle
240  !
241  ! -- calculate new and old water volumes
242  vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
243  this%fmi%gwfsat(n) * this%thetam(n)
244  vold = vnew
245  if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) * delt
246  if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) * delt
247  !
248  ! -- add terms to diagonal and rhs accumulators
249  hhcof = -vnew * tled
250  rrhs = -vold * tled * cold(n)
251  idiag = this%dis%con%ia(n)
252  call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
253  rhs(n) = rhs(n) + rrhs
254  end do

◆ mst_ot_dv()

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

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

985  ! -- dummy
986  class(GwtMstType) :: this
987  integer(I4B), intent(in) :: idvsave
988  ! -- local
989  character(len=1) :: cdatafmp = ' ', editdesc = ' '
990  integer(I4B) :: ibinun
991  integer(I4B) :: iprint
992  integer(I4B) :: nvaluesp
993  integer(I4B) :: nwidthp
994  real(DP) :: dinact
995 
996  ! Set unit number for sorbate output
997  if (this%ioutsorbate /= 0) then
998  ibinun = 1
999  else
1000  ibinun = 0
1001  end if
1002  if (idvsave == 0) ibinun = 0
1003 
1004  ! save sorbate concentration array
1005  if (ibinun /= 0) then
1006  iprint = 0
1007  dinact = dhnoflo
1008  if (this%ioutsorbate /= 0) then
1009  ibinun = this%ioutsorbate
1010  call this%dis%record_array(this%csrb, this%iout, iprint, ibinun, &
1011  ' SORBATE', cdatafmp, nvaluesp, &
1012  nwidthp, editdesc, dinact)
1013  end if
1014  end if
1015 

◆ 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
thisGwtMstType object
[in]icbcflflag and unit number for cell-by-cell output
[in]icbcunflag indication if cell-by-cell data should be saved

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

932  ! -- dummy
933  class(GwtMstType) :: this !< GwtMstType object
934  integer(I4B), intent(in) :: icbcfl !< flag and unit number for cell-by-cell output
935  integer(I4B), intent(in) :: icbcun !< flag indication if cell-by-cell data should be saved
936  ! -- local
937  integer(I4B) :: ibinun
938  integer(I4B) :: iprint, nvaluesp, nwidthp
939  character(len=1) :: cdatafmp = ' ', editdesc = ' '
940  real(DP) :: dinact
941  !
942  ! -- Set unit number for binary output
943  if (this%ipakcb < 0) then
944  ibinun = icbcun
945  elseif (this%ipakcb == 0) then
946  ibinun = 0
947  else
948  ibinun = this%ipakcb
949  end if
950  if (icbcfl == 0) ibinun = 0
951  !
952  ! -- Record the storage rate if requested
953  if (ibinun /= 0) then
954  iprint = 0
955  dinact = dzero
956  !
957  ! -- sto
958  call this%dis%record_array(this%ratesto, this%iout, iprint, -ibinun, &
959  budtxt(1), cdatafmp, nvaluesp, &
960  nwidthp, editdesc, dinact)
961  !
962  ! -- dcy
963  if (this%idcy /= decay_off) &
964  call this%dis%record_array(this%ratedcy, this%iout, iprint, -ibinun, &
965  budtxt(2), cdatafmp, nvaluesp, &
966  nwidthp, editdesc, dinact)
967  !
968  ! -- srb
969  if (this%isrb /= sorption_off) &
970  call this%dis%record_array(this%ratesrb, this%iout, iprint, -ibinun, &
971  budtxt(3), cdatafmp, nvaluesp, &
972  nwidthp, editdesc, dinact)
973  !
974  ! -- dcy srb
975  if (this%isrb /= sorption_off .and. this%idcy /= decay_off) &
976  call this%dis%record_array(this%ratedcys, this%iout, iprint, -ibinun, &
977  budtxt(4), cdatafmp, nvaluesp, &
978  nwidthp, editdesc, dinact)
979  end if

◆ mst_srb_term()

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

Subroutine to calculate sorption terms

Parameters
[in]isrbsorption flag 1, 2, 3 are linear, freundlich, and langmuir
[in]volfracmvolume fraction of mobile domain (fhat_m)
[in]rhobmbulk density of mobile domain (rhob_m)
[in]vcellvolume of cell
[in]tledone over time step length
[in]cnewconcentration at end of this time step
[in]coldconcentration at end of last time step
[in]swnewcell saturation at end of this time step
[in]swoldcell saturation at end of last time step
[in]const1distribution coefficient or freundlich or langmuir constant
[in]const2zero, freundlich exponent, or langmuir sorption sites
[out]ratecalculated sorption rate
[out]hcofvaldiagonal contribution to solution coefficient matrix
[out]rhsvalcontribution to solution right-hand-side

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

383  ! -- dummy
384  integer(I4B), intent(in) :: isrb !< sorption flag 1, 2, 3 are linear, freundlich, and langmuir
385  real(DP), intent(in) :: volfracm !< volume fraction of mobile domain (fhat_m)
386  real(DP), intent(in) :: rhobm !< bulk density of mobile domain (rhob_m)
387  real(DP), intent(in) :: vcell !< volume of cell
388  real(DP), intent(in) :: tled !< one over time step length
389  real(DP), intent(in) :: cnew !< concentration at end of this time step
390  real(DP), intent(in) :: cold !< concentration at end of last time step
391  real(DP), intent(in) :: swnew !< cell saturation at end of this time step
392  real(DP), intent(in) :: swold !< cell saturation at end of last time step
393  real(DP), intent(in) :: const1 !< distribution coefficient or freundlich or langmuir constant
394  real(DP), intent(in) :: const2 !< zero, freundlich exponent, or langmuir sorption sites
395  real(DP), intent(out), optional :: rate !< calculated sorption rate
396  real(DP), intent(out), optional :: hcofval !< diagonal contribution to solution coefficient matrix
397  real(DP), intent(out), optional :: rhsval !< contribution to solution right-hand-side
398  ! -- local
399  real(DP) :: term
400  real(DP) :: derv
401  real(DP) :: cbarnew
402  real(DP) :: cbarold
403  real(DP) :: cavg
404  real(DP) :: cbaravg
405  real(DP) :: swavg
406  !
407  ! -- Calculate based on type of sorption
408  if (isrb == sorption_linear) then
409  ! -- linear
410  term = -volfracm * rhobm * vcell * tled * const1
411  if (present(hcofval)) hcofval = term * swnew
412  if (present(rhsval)) rhsval = term * swold * cold
413  if (present(rate)) rate = term * swnew * cnew - term * swold * cold
414  else
415  !
416  ! -- calculate average aqueous concentration
417  cavg = dhalf * (cold + cnew)
418  !
419  ! -- set values based on isotherm
420  select case (isrb)
421  case (sorption_freund)
422  ! -- freundlich
423  cbarnew = get_freundlich_conc(cnew, const1, const2)
424  cbarold = get_freundlich_conc(cold, const1, const2)
425  derv = get_freundlich_derivative(cavg, const1, const2)
426  case (sorption_lang)
427  ! -- langmuir
428  cbarnew = get_langmuir_conc(cnew, const1, const2)
429  cbarold = get_langmuir_conc(cold, const1, const2)
430  derv = get_langmuir_derivative(cavg, const1, const2)
431  end select
432  !
433  ! -- calculate hcof, rhs, and rate for freundlich and langmuir
434  term = -volfracm * rhobm * vcell * tled
435  cbaravg = (cbarold + cbarnew) * dhalf
436  swavg = (swnew + swold) * dhalf
437  if (present(hcofval)) then
438  hcofval = term * derv * swavg
439  end if
440  if (present(rhsval)) then
441  rhsval = term * derv * swavg * cold - term * cbaravg * (swnew - swold)
442  end if
443  if (present(rate)) then
444  rate = term * derv * swavg * (cnew - cold) &
445  + term * cbaravg * (swnew - swold)
446  end if
447  end if
Here is the call graph for this function:
Here is the caller graph for this function:

◆ source_data()

subroutine gwtmstmodule::source_data ( class(gwtmsttype)  this)

Method to source data for the package.

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

1283  ! -- modules
1284  use memorymanagermodule, only: get_isize, mem_reallocate
1285  use memorymanagerextmodule, only: mem_set_value
1286  use gwtmstinputmodule, only: gwtmstparamfoundtype
1287  ! -- dummy
1288  class(GwtMsttype) :: this
1289  ! -- locals
1290  character(len=LINELENGTH) :: errmsg
1291  type(GwtMstParamFoundType) :: found
1292  integer(I4B) :: n, asize
1293  integer(I4B), dimension(:), pointer, contiguous :: map
1294  !
1295  ! -- set map to convert user input data into reduced data
1296  map => null()
1297  if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
1298  !
1299  ! -- reallocate
1300  if (this%isrb == sorption_off) then
1301  call get_isize('BULK_DENSITY', this%input_mempath, asize)
1302  if (asize > 0) &
1303  call mem_reallocate(this%bulk_density, this%dis%nodes, &
1304  'BULK_DENSITY', this%memoryPath)
1305  call get_isize('DISTCOEF', this%input_mempath, asize)
1306  if (asize > 0) &
1307  call mem_reallocate(this%distcoef, this%dis%nodes, 'DISTCOEF', &
1308  this%memoryPath)
1309  end if
1310  if (this%idcy == decay_off) then
1311  call get_isize('DECAY', this%input_mempath, asize)
1312  if (asize > 0) &
1313  call mem_reallocate(this%decay, this%dis%nodes, 'DECAY', this%memoryPath)
1314  end if
1315  call get_isize('DECAY_SORBED', this%input_mempath, asize)
1316  if (asize > 0) then
1317  call mem_reallocate(this%decay_sorbed, this%dis%nodes, &
1318  'DECAY_SORBED', this%memoryPath)
1319  end if
1320  if (this%isrb == sorption_off .or. this%isrb == sorption_linear) then
1321  call get_isize('SP2', this%input_mempath, asize)
1322  if (asize > 0) &
1323  call mem_reallocate(this%sp2, this%dis%nodes, 'SP2', this%memoryPath)
1324  end if
1325  !
1326  ! -- update defaults with memory sourced values
1327  call mem_set_value(this%porosity, 'POROSITY', this%input_mempath, map, &
1328  found%porosity)
1329  call mem_set_value(this%decay, 'DECAY', this%input_mempath, map, &
1330  found%decay)
1331  call mem_set_value(this%decay_sorbed, 'DECAY_SORBED', this%input_mempath, &
1332  map, found%decay_sorbed)
1333  call mem_set_value(this%bulk_density, 'BULK_DENSITY', this%input_mempath, &
1334  map, found%bulk_density)
1335  call mem_set_value(this%distcoef, 'DISTCOEF', this%input_mempath, map, &
1336  found%distcoef)
1337  call mem_set_value(this%sp2, 'SP2', this%input_mempath, map, &
1338  found%sp2)
1339 
1340  ! -- log options
1341  if (this%iout > 0) then
1342  call this%log_data(found)
1343  end if
1344 
1345  ! -- Check for required porosity
1346  if (.not. found%porosity) then
1347  write (errmsg, '(a)') 'POROSITY not specified in GRIDDATA block.'
1348  call store_error(errmsg)
1349  end if
1350 
1351  ! -- Check for required sorption variables
1352  if (this%isrb == sorption_linear .or. this%isrb == sorption_freund .or. &
1353  this%isrb == sorption_lang) then
1354  if (.not. found%bulk_density) then
1355  write (errmsg, '(a)') 'Sorption is active but BULK_DENSITY &
1356  &not specified. BULK_DENSITY must be specified in GRIDDATA block.'
1357  call store_error(errmsg)
1358  end if
1359  if (.not. found%distcoef) then
1360  write (errmsg, '(a)') 'Sorption is active but distribution &
1361  &coefficient not specified. DISTCOEF must be specified in &
1362  &GRIDDATA block.'
1363  call store_error(errmsg)
1364  end if
1365  if (this%isrb == sorption_freund .or. this%isrb == sorption_lang) then
1366  if (.not. found%sp2) then
1367  write (errmsg, '(a)') 'Freundlich or langmuir sorption is active &
1368  &but SP2 not specified. SP2 must be specified in &
1369  &GRIDDATA block.'
1370  call store_error(errmsg)
1371  end if
1372  end if
1373  else
1374  if (found%bulk_density) then
1375  write (warnmsg, '(a)') 'Sorption is not active but &
1376  &BULK_DENSITY was specified. BULK_DENSITY will have no affect on &
1377  &simulation results.'
1378  call store_warning(warnmsg)
1379  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1380  end if
1381  if (found%distcoef) then
1382  write (warnmsg, '(a)') 'Sorption is not active but &
1383  &distribution coefficient was specified. DISTCOEF will have &
1384  &no affect on simulation results.'
1385  call store_warning(warnmsg)
1386  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1387  end if
1388  if (found%sp2) then
1389  write (warnmsg, '(a)') 'Sorption is not active but &
1390  &SP2 was specified. SP2 will have &
1391  &no affect on simulation results.'
1392  call store_warning(warnmsg)
1393  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1394  end if
1395  end if
1396 
1397  ! -- Check for required decay/production rate coefficients
1398  if (this%idcy /= decay_off) then
1399  if (.not. found%decay) then
1400  write (errmsg, '(a)') 'First or zero order decay is &
1401  &active but the first rate coefficient is not specified. DECAY &
1402  &must be specified in GRIDDATA block.'
1403  call store_error(errmsg)
1404  end if
1405  if (.not. found%decay_sorbed) then
1406  !
1407  ! -- If DECAY_SORBED not specified and sorption is active, then
1408  ! terminate with an error
1409  if (this%isrb == sorption_linear .or. this%isrb == sorption_freund .or. &
1410  this%isrb == sorption_lang) then
1411  write (errmsg, '(a)') 'DECAY_SORBED not provided in GRIDDATA &
1412  &block but decay and sorption are active. Specify DECAY_SORBED &
1413  &in GRIDDATA block.'
1414  call store_error(errmsg)
1415  end if
1416  end if
1417  else
1418  if (found%decay) then
1419  write (warnmsg, '(a)') 'First- or zero-order decay &
1420  &is not active but decay was specified. DECAY will &
1421  &have no affect on simulation results.'
1422  call store_warning(warnmsg)
1423  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1424  end if
1425  if (found%decay_sorbed) then
1426  write (warnmsg, '(a)') 'First- or zero-order decay &
1427  &is not active but DECAY_SORBED was specified. &
1428  &DECAY_SORBED will have no affect on simulation results.'
1429  call store_warning(warnmsg)
1430  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1431  end if
1432  end if
1433 
1434  ! -- terminate if errors
1435  if (count_errors() > 0) then
1436  call store_error_filename(this%input_fname)
1437  end if
1438 
1439  ! -- initialize thetam from porosity
1440  do n = 1, size(this%porosity)
1441  this%thetam(n) = this%porosity(n)
1442  end do
memorymanagerextmodule::mem_set_value
Definition: MemoryManagerExt.f90:14
memorymanagermodule::mem_reallocate
Definition: MemoryManager.f90:79
memorymanagermodule::get_isize
subroutine, public get_isize(name, mem_path, isize)
@ brief Get the number of elements for this variable
Definition: MemoryManager.f90:265
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 1175 of file gwt-mst.f90.

1176  ! -- modules
1177  use constantsmodule, only: lenvarname
1178  use openspecmodule, only: access, form
1179  use inputoutputmodule, only: getunit, openfile
1180  use memorymanagerextmodule, only: mem_set_value
1181  use gwtmstinputmodule, only: gwtmstparamfoundtype
1182  ! -- dummy
1183  class(GwtMstType) :: this
1184  ! -- locals
1185  type(GwtMstParamFoundType) :: found
1186  character(len=LENVARNAME), dimension(3) :: sorption_method = &
1187  &[character(len=LENVARNAME) :: 'LINEAR', 'FREUNDLICH', 'LANGMUIR']
1188  character(len=linelength) :: fname
1189  !
1190  ! -- update defaults with memory sourced values
1191  call mem_set_value(this%ipakcb, 'SAVE_FLOWS', this%input_mempath, &
1192  found%save_flows)
1193  call mem_set_value(this%idcy, 'ORDER1_DECAY', this%input_mempath, &
1194  found%order1_decay)
1195  call mem_set_value(this%idcy, 'ORDER0_DECAY', this%input_mempath, &
1196  found%order0_decay)
1197  call mem_set_value(this%isrb, 'SORPTION', this%input_mempath, &
1198  sorption_method, found%sorption)
1199  call mem_set_value(fname, 'SORBATEFILE', this%input_mempath, &
1200  found%sorbatefile)
1201 
1202  ! -- found side effects
1203  if (found%save_flows) this%ipakcb = -1
1204  if (found%order1_decay) this%idcy = decay_first_order
1205  if (found%order0_decay) this%idcy = decay_zero_order
1206  if (found%sorption) then
1207  if (this%isrb == izero) then
1208  call store_error('Unknown sorption type was specified. &
1209  &Sorption must be specified as LINEAR, &
1210  &FREUNDLICH, or LANGMUIR.')
1211  call store_error_filename(this%input_fname)
1212  end if
1213  end if
1214  if (found%sorbatefile) then
1215  this%ioutsorbate = getunit()
1216  call openfile(this%ioutsorbate, this%iout, fname, 'DATA(BINARY)', &
1217  form, access, 'REPLACE', mode_opt=mnormal)
1218  end if
1219  !
1220  ! -- log options
1221  if (this%iout > 0) then
1222  call this%log_options(found, fname)
1223  end if
constantsmodule::lenvarname
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
inputoutputmodule::getunit
integer(i4b) function, public getunit()
Get a free unit number.
Definition: InputOutput.f90:170
inputoutputmodule::openfile
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
openspecmodule::access
character(len=20) access
Definition: OpenSpec.f90:7
openspecmodule::form
character(len=20) form
Definition: OpenSpec.f90:7
Here is the call graph for this function:

Variable Documentation

◆ budtxt

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

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

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

◆ nbditems

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

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

27  integer(I4B), parameter :: NBDITEMS = 4