MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
gwt-mst.f90
Go to the documentation of this file.
1 !> -- @ brief Mobile Storage and Transfer (MST) Module
2 !!
3 !! The GwtMstModule contains the GwtMstType, which is the
4 !! derived type responsible for adding the effects of
5 !! 1. Changes in dissolved solute mass
6 !! 2. Decay of dissolved solute mass
7 !! 3. Sorption
8 !! 4. Decay of sorbed solute mass
9 !<
11 
12  use kindmodule, only: dp, i4b
16  use simmodule, only: store_error, count_errors, &
20  use basedismodule, only: disbasetype
21  use tspfmimodule, only: tspfmitype
22 
23  implicit none
24  public :: gwtmsttype
25  public :: mst_cr
26  !
27  integer(I4B), parameter :: nbditems = 4
28  character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt
29  data budtxt/' STORAGE-AQUEOUS', ' DECAY-AQUEOUS', &
30  ' STORAGE-SORBED', ' DECAY-SORBED'/
31 
32  !> @brief Enumerator that defines the decay options
33  !<
34  ENUM, BIND(C)
35  ENUMERATOR :: decay_off = 0 !< Decay (or production) of mass inactive (default)
36  ENUMERATOR :: decay_first_order = 1 !< First-order decay
37  ENUMERATOR :: decay_zero_order = 2 !< Zeroth-order decay
38  ENUMERATOR :: sorption_off = 0 !< Sorption is inactive (default)
39  ENUMERATOR :: sorption_linear = 1 !< Linear sorption between aqueous and solid phases
40  ENUMERATOR :: sorption_freund = 2 !< Freundlich sorption between aqueous and solid phases
41  ENUMERATOR :: sorption_lang = 3 !< Langmuir sorption between aqueous and solid phases
42  END ENUM
43 
44  !> @ brief Mobile storage and transfer
45  !!
46  !! Data and methods for handling changes in solute storage,
47  !! decay of dissolved solute mass, sorption, and decay of
48  !! sorbed mass.
49  !<
50  type, extends(numericalpackagetype) :: gwtmsttype
51  !
52  ! -- storage
53  real(dp), dimension(:), pointer, contiguous :: porosity => null() !< mobile porosity defined as volume mobile voids per volume of mobile domain
54  real(dp), dimension(:), pointer, contiguous :: thetam => null() !< mobile porosity defined as volume mobile voids per volume of aquifer
55  real(dp), dimension(:), pointer, contiguous :: volfracim => null() !< sum of all immobile domain volume fractions
56  real(dp), dimension(:), pointer, contiguous :: ratesto => null() !< rate of mobile storage
57  !
58  ! -- decay
59  integer(I4B), pointer :: idcy => null() !< order of decay rate (0:none, 1:first, 2:zero)
60  real(dp), dimension(:), pointer, contiguous :: decay => null() !< first or zero order decay rate (aqueous)
61  real(dp), dimension(:), pointer, contiguous :: decay_sorbed => null() !< first or zero order decay rate (sorbed)
62  real(dp), dimension(:), pointer, contiguous :: ratedcy => null() !< rate of decay
63  real(dp), dimension(:), pointer, contiguous :: decaylast => null() !< decay rate used for last iteration (needed for zero order decay)
64  real(dp), dimension(:), pointer, contiguous :: decayslast => null() !< sorbed decay rate used for last iteration (needed for zero order decay)
65  !
66  ! -- sorption
67  integer(I4B), pointer :: isrb => null() !< sorption active flag (0:off, 1:linear, 2:freundlich, 3:langmuir)
68  integer(I4B), pointer :: ioutsorbate => null() !< unit number for sorbate concentration output
69  real(dp), dimension(:), pointer, contiguous :: bulk_density => null() !< bulk density of mobile domain; mass of mobile domain solid per aquifer volume
70  real(dp), dimension(:), pointer, contiguous :: distcoef => null() !< kd distribution coefficient
71  real(dp), dimension(:), pointer, contiguous :: sp2 => null() !< second sorption parameter
72  real(dp), dimension(:), pointer, contiguous :: ratesrb => null() !< rate of sorption
73  real(dp), dimension(:), pointer, contiguous :: ratedcys => null() !< rate of sorbed mass decay
74  real(dp), dimension(:), pointer, contiguous :: csrb => null() !< sorbate concentration
75  !
76  ! -- misc
77  integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound
78  type(tspfmitype), pointer :: fmi => null() !< pointer to fmi object
79 
80  contains
81 
82  procedure :: mst_ar
83  procedure :: mst_fc
84  procedure :: mst_fc_sto
85  procedure :: mst_fc_dcy
86  procedure :: mst_fc_srb
87  procedure :: mst_fc_dcy_srb
88  procedure :: mst_cq
89  procedure :: mst_cq_sto
90  procedure :: mst_cq_dcy
91  procedure :: mst_cq_srb
92  procedure :: mst_cq_dcy_srb
93  procedure :: mst_calc_csrb
94  procedure :: mst_bd
95  procedure :: mst_ot_flow
96  procedure :: mst_ot_dv
97  procedure :: mst_da
98  procedure :: allocate_scalars
99  procedure :: addto_volfracim
100  procedure :: get_volfracm
101  procedure, private :: allocate_arrays
102  procedure, private :: read_options
103  procedure, private :: read_data
104 
105  end type gwtmsttype
106 
107 contains
108 
109  !> @ brief Create a new package object
110  !!
111  !! Create a new MST object
112  !<
113  subroutine mst_cr(mstobj, name_model, inunit, iout, fmi)
114  ! -- dummy
115  type(gwtmsttype), pointer :: mstobj !< unallocated new mst object to create
116  character(len=*), intent(in) :: name_model !< name of the model
117  integer(I4B), intent(in) :: inunit !< unit number of WEL package input file
118  integer(I4B), intent(in) :: iout !< unit number of model listing file
119  type(tspfmitype), intent(in), target :: fmi !< fmi package for this GWT model
120  !
121  ! -- Create the object
122  allocate (mstobj)
123  !
124  ! -- create name and memory path
125  call mstobj%set_names(1, name_model, 'MST', 'MST')
126  !
127  ! -- Allocate scalars
128  call mstobj%allocate_scalars()
129  !
130  ! -- Set variables
131  mstobj%inunit = inunit
132  mstobj%iout = iout
133  mstobj%fmi => fmi
134  !
135  ! -- Initialize block parser
136  call mstobj%parser%Initialize(mstobj%inunit, mstobj%iout)
137  end subroutine mst_cr
138 
139  !> @ brief Allocate and read method for package
140  !!
141  !! Method to allocate and read static data for the package.
142  !<
143  subroutine mst_ar(this, dis, ibound)
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 UNIT ', i0, //)"
154  !
155  ! --print a message identifying the mobile storage and transfer package.
156  write (this%iout, fmtmst) this%inunit
157  !
158  ! -- Read options
159  call this%read_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  ! -- read the data block
169  call this%read_data()
170  end subroutine mst_ar
171 
172  !> @ brief Fill coefficient method for package
173  !!
174  !! Method to calculate and fill coefficients for the package.
175  !<
176  subroutine mst_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, &
177  rhs, kiter)
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
209  end subroutine mst_fc
210 
211  !> @ brief Fill storage coefficient method for package
212  !!
213  !! Method to calculate and fill storage coefficients for the package.
214  !<
215  subroutine mst_fc_sto(this, nodes, cold, nja, matrix_sln, idxglo, rhs)
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
255  end subroutine mst_fc_sto
256 
257  !> @ brief Fill decay coefficient method for package
258  !!
259  !! Method to calculate and fill decay coefficients for the package.
260  !<
261  subroutine mst_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, &
262  idxglo, rhs, kiter)
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  hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
300  call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
301  case (decay_zero_order)
302  !
303  ! -- Call function to get zero-order decay rate, which may be changed
304  ! from the user-specified rate to prevent negative concentrations
305  decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), &
306  kiter, cold(n), cnew(n), delt)
307  this%decaylast(n) = decay_rate
308  rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
309  rhs(n) = rhs(n) + rrhs
310  end select
311  !
312  end do
313  end subroutine mst_fc_dcy
314 
315  !> @ brief Fill sorption coefficient method for package
316  !!
317  !! Method to calculate and fill sorption coefficients for the package.
318  !<
319  subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, &
320  cnew)
321  ! -- modules
322  use tdismodule, only: delt
323  ! -- dummy
324  class(gwtmsttype) :: this !< GwtMstType object
325  integer, intent(in) :: nodes !< number of nodes
326  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
327  integer(I4B), intent(in) :: nja !< number of GWT connections
328  class(matrixbasetype), pointer :: matrix_sln !< solution coefficient matrix
329  integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
330  real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
331  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
332  ! -- local
333  integer(I4B) :: n, idiag
334  real(DP) :: tled
335  real(DP) :: hhcof, rrhs
336  real(DP) :: swt, swtpdt
337  real(DP) :: vcell
338  real(DP) :: const1
339  real(DP) :: const2
340  real(DP) :: volfracm
341  real(DP) :: rhobm
342  !
343  ! -- set variables
344  tled = done / delt
345  !
346  ! -- loop through and calculate sorption contribution to hcof and rhs
347  do n = 1, this%dis%nodes
348  !
349  ! -- skip if transport inactive
350  if (this%ibound(n) <= 0) cycle
351  !
352  ! -- assign variables
353  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
354  swtpdt = this%fmi%gwfsat(n)
355  swt = this%fmi%gwfsatold(n, delt)
356  idiag = this%dis%con%ia(n)
357  const1 = this%distcoef(n)
358  const2 = 0.
359  if (this%isrb == sorption_freund .or. this%isrb == sorption_lang) then
360  const2 = this%sp2(n)
361  end if
362  volfracm = this%get_volfracm(n)
363  rhobm = this%bulk_density(n)
364  call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
365  cold(n), swtpdt, swt, const1, const2, &
366  hcofval=hhcof, rhsval=rrhs)
367  !
368  ! -- Add hhcof to diagonal and rrhs to right-hand side
369  call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
370  rhs(n) = rhs(n) + rrhs
371  !
372  end do
373  end subroutine mst_fc_srb
374 
375  !> @ brief Calculate sorption terms
376  !!
377  !! Subroutine to calculate sorption terms
378  !<
379  subroutine mst_srb_term(isrb, volfracm, rhobm, vcell, tled, cnew, cold, &
380  swnew, swold, const1, const2, rate, hcofval, rhsval)
381  ! -- dummy
382  integer(I4B), intent(in) :: isrb !< sorption flag 1, 2, 3 are linear, freundlich, and langmuir
383  real(DP), intent(in) :: volfracm !< volume fraction of mobile domain (fhat_m)
384  real(DP), intent(in) :: rhobm !< bulk density of mobile domain (rhob_m)
385  real(DP), intent(in) :: vcell !< volume of cell
386  real(DP), intent(in) :: tled !< one over time step length
387  real(DP), intent(in) :: cnew !< concentration at end of this time step
388  real(DP), intent(in) :: cold !< concentration at end of last time step
389  real(DP), intent(in) :: swnew !< cell saturation at end of this time step
390  real(DP), intent(in) :: swold !< cell saturation at end of last time step
391  real(DP), intent(in) :: const1 !< distribution coefficient or freundlich or langmuir constant
392  real(DP), intent(in) :: const2 !< zero, freundlich exponent, or langmuir sorption sites
393  real(DP), intent(out), optional :: rate !< calculated sorption rate
394  real(DP), intent(out), optional :: hcofval !< diagonal contribution to solution coefficient matrix
395  real(DP), intent(out), optional :: rhsval !< contribution to solution right-hand-side
396  ! -- local
397  real(DP) :: term
398  real(DP) :: derv
399  real(DP) :: cbarnew
400  real(DP) :: cbarold
401  real(DP) :: cavg
402  real(DP) :: cbaravg
403  real(DP) :: swavg
404  !
405  ! -- Calculate based on type of sorption
406  if (isrb == sorption_linear) then
407  ! -- linear
408  term = -volfracm * rhobm * vcell * tled * const1
409  if (present(hcofval)) hcofval = term * swnew
410  if (present(rhsval)) rhsval = term * swold * cold
411  if (present(rate)) rate = term * swnew * cnew - term * swold * cold
412  else
413  !
414  ! -- calculate average aqueous concentration
415  cavg = dhalf * (cold + cnew)
416  !
417  ! -- set values based on isotherm
418  select case (isrb)
419  case (sorption_freund)
420  ! -- freundlich
421  cbarnew = get_freundlich_conc(cnew, const1, const2)
422  cbarold = get_freundlich_conc(cold, const1, const2)
423  derv = get_freundlich_derivative(cavg, const1, const2)
424  case (sorption_lang)
425  ! -- langmuir
426  cbarnew = get_langmuir_conc(cnew, const1, const2)
427  cbarold = get_langmuir_conc(cold, const1, const2)
428  derv = get_langmuir_derivative(cavg, const1, const2)
429  end select
430  !
431  ! -- calculate hcof, rhs, and rate for freundlich and langmuir
432  term = -volfracm * rhobm * vcell * tled
433  cbaravg = (cbarold + cbarnew) * dhalf
434  swavg = (swnew + swold) * dhalf
435  if (present(hcofval)) then
436  hcofval = term * derv * swavg
437  end if
438  if (present(rhsval)) then
439  rhsval = term * derv * swavg * cold - term * cbaravg * (swnew - swold)
440  end if
441  if (present(rate)) then
442  rate = term * derv * swavg * (cnew - cold) &
443  + term * cbaravg * (swnew - swold)
444  end if
445  end if
446  end subroutine mst_srb_term
447 
448  !> @ brief Fill sorption-decay coefficient method for package
449  !!
450  !! Method to calculate and fill sorption-decay coefficients for the package.
451  !<
452  subroutine mst_fc_dcy_srb(this, nodes, cold, nja, matrix_sln, idxglo, &
453  rhs, cnew, kiter)
454  ! -- modules
455  use tdismodule, only: delt
456  ! -- dummy
457  class(gwtmsttype) :: this !< GwtMstType object
458  integer, intent(in) :: nodes !< number of nodes
459  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
460  integer(I4B), intent(in) :: nja !< number of GWT connections
461  class(matrixbasetype), pointer :: matrix_sln !< solution coefficient matrix
462  integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global)
463  real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model
464  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
465  integer(I4B), intent(in) :: kiter !< solution outer iteration number
466  ! -- local
467  integer(I4B) :: n, idiag
468  real(DP) :: hhcof, rrhs
469  real(DP) :: vcell
470  real(DP) :: swnew
471  real(DP) :: distcoef
472  real(DP) :: volfracm
473  real(DP) :: rhobm
474  real(DP) :: term
475  real(DP) :: csrb
476  real(DP) :: decay_rate
477  real(DP) :: csrbold
478  real(DP) :: csrbnew
479  !
480  ! -- loop through and calculate sorption contribution to hcof and rhs
481  do n = 1, this%dis%nodes
482  !
483  ! -- skip if transport inactive
484  if (this%ibound(n) <= 0) cycle
485  !
486  ! -- set variables
487  hhcof = dzero
488  rrhs = dzero
489  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
490  swnew = this%fmi%gwfsat(n)
491  distcoef = this%distcoef(n)
492  idiag = this%dis%con%ia(n)
493  volfracm = this%get_volfracm(n)
494  rhobm = this%bulk_density(n)
495  term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
496  !
497  ! -- add sorbed mass decay rate terms to accumulators
498  select case (this%idcy)
499  case (decay_first_order)
500  !
501  select case (this%isrb)
502  case (sorption_linear)
503  !
504  ! -- first order decay rate is a function of concentration, so add
505  ! to left hand side
506  hhcof = -term * distcoef
507  case (sorption_freund)
508  !
509  ! -- nonlinear Freundlich sorption, so add to RHS
510  csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
511  rrhs = term * csrb
512  case (sorption_lang)
513  !
514  ! -- nonlinear Lanmuir sorption, so add to RHS
515  csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
516  rrhs = term * csrb
517  end select
518  case (decay_zero_order)
519  !
520  ! -- call function to get zero-order decay rate, which may be changed
521  ! from the user-specified rate to prevent negative concentrations
522  if (distcoef > dzero) then
523  select case (this%isrb)
524  case (sorption_linear)
525  csrbold = cold(n) * distcoef
526  csrbnew = cnew(n) * distcoef
527  case (sorption_freund)
528  csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n))
529  csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
530  case (sorption_lang)
531  csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n))
532  csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
533  end select
534  !
535  decay_rate = get_zero_order_decay(this%decay_sorbed(n), &
536  this%decayslast(n), &
537  kiter, csrbold, csrbnew, delt)
538  this%decayslast(n) = decay_rate
539  rrhs = decay_rate * volfracm * rhobm * swnew * vcell
540  end if
541  end select
542  !
543  ! -- Add hhcof to diagonal and rrhs to right-hand side
544  call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
545  rhs(n) = rhs(n) + rrhs
546  !
547  end do
548  end subroutine mst_fc_dcy_srb
549 
550  !> @ brief Calculate flows for package
551  !!
552  !! Method to calculate flows for the package.
553  !<
554  subroutine mst_cq(this, nodes, cnew, cold, flowja)
555  ! -- dummy
556  class(gwtmsttype) :: this !< GwtMstType object
557  integer(I4B), intent(in) :: nodes !< number of nodes
558  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
559  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
560  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
561  !
562  ! - storage
563  call this%mst_cq_sto(nodes, cnew, cold, flowja)
564  !
565  ! -- decay
566  if (this%idcy /= decay_off) then
567  call this%mst_cq_dcy(nodes, cnew, cold, flowja)
568  end if
569  !
570  ! -- sorption
571  if (this%isrb /= sorption_off) then
572  call this%mst_cq_srb(nodes, cnew, cold, flowja)
573  end if
574  !
575  ! -- decay sorbed
576  if (this%isrb /= sorption_off .and. this%idcy /= decay_off) then
577  call this%mst_cq_dcy_srb(nodes, cnew, cold, flowja)
578  end if
579  !
580  ! -- calculate csrb
581  if (this%isrb /= sorption_off) then
582  call this%mst_calc_csrb(cnew)
583  end if
584  end subroutine mst_cq
585 
586  !> @ brief Calculate storage terms for package
587  !!
588  !! Method to calculate storage terms for the package.
589  !<
590  subroutine mst_cq_sto(this, nodes, cnew, cold, flowja)
591  ! -- modules
592  use tdismodule, only: delt
593  ! -- dummy
594  class(gwtmsttype) :: this !< GwtMstType object
595  integer(I4B), intent(in) :: nodes !< number of nodes
596  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
597  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
598  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
599  ! -- local
600  integer(I4B) :: n
601  integer(I4B) :: idiag
602  real(DP) :: rate
603  real(DP) :: tled
604  real(DP) :: vnew, vold
605  real(DP) :: hhcof, rrhs
606  !
607  ! -- initialize
608  tled = done / delt
609  !
610  ! -- Calculate storage change
611  do n = 1, nodes
612  this%ratesto(n) = dzero
613  !
614  ! -- skip if transport inactive
615  if (this%ibound(n) <= 0) cycle
616  !
617  ! -- calculate new and old water volumes
618  vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
619  this%fmi%gwfsat(n) * this%thetam(n)
620  vold = vnew
621  if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) * delt
622  if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) * delt
623  !
624  ! -- calculate rate
625  hhcof = -vnew * tled
626  rrhs = -vold * tled * cold(n)
627  rate = hhcof * cnew(n) - rrhs
628  this%ratesto(n) = rate
629  idiag = this%dis%con%ia(n)
630  flowja(idiag) = flowja(idiag) + rate
631  end do
632  end subroutine mst_cq_sto
633 
634  !> @ brief Calculate decay terms for package
635  !!
636  !! Method to calculate decay terms for the package.
637  !<
638  subroutine mst_cq_dcy(this, nodes, cnew, cold, flowja)
639  ! -- modules
640  use tdismodule, only: delt
641  ! -- dummy
642  class(gwtmsttype) :: this !< GwtMstType object
643  integer(I4B), intent(in) :: nodes !< number of nodes
644  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
645  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
646  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
647  ! -- local
648  integer(I4B) :: n
649  integer(I4B) :: idiag
650  real(DP) :: rate
651  real(DP) :: swtpdt
652  real(DP) :: hhcof, rrhs
653  real(DP) :: vcell
654  real(DP) :: decay_rate
655  !
656  ! -- initialize
657  !
658  ! -- Calculate decay change
659  do n = 1, nodes
660  !
661  ! -- skip if transport inactive
662  this%ratedcy(n) = dzero
663  if (this%ibound(n) <= 0) cycle
664  !
665  ! -- calculate new and old water volumes
666  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
667  swtpdt = this%fmi%gwfsat(n)
668  !
669  ! -- calculate decay gains and losses
670  rate = dzero
671  hhcof = dzero
672  rrhs = dzero
673  if (this%idcy == decay_first_order) then
674  hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
675  elseif (this%idcy == decay_zero_order) then
676  decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), &
677  0, cold(n), cnew(n), delt)
678  rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
679  end if
680  rate = hhcof * cnew(n) - rrhs
681  this%ratedcy(n) = rate
682  idiag = this%dis%con%ia(n)
683  flowja(idiag) = flowja(idiag) + rate
684  !
685  end do
686  end subroutine mst_cq_dcy
687 
688  !> @ brief Calculate sorption terms for package
689  !!
690  !! Method to calculate sorption terms for the package.
691  !<
692  subroutine mst_cq_srb(this, nodes, cnew, cold, flowja)
693  ! -- modules
694  use tdismodule, only: delt
695  ! -- dummy
696  class(gwtmsttype) :: this !< GwtMstType object
697  integer(I4B), intent(in) :: nodes !< number of nodes
698  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
699  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
700  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
701  ! -- local
702  integer(I4B) :: n
703  integer(I4B) :: idiag
704  real(DP) :: rate
705  real(DP) :: tled
706  real(DP) :: swt, swtpdt
707  real(DP) :: vcell
708  real(DP) :: volfracm
709  real(DP) :: rhobm
710  real(DP) :: const1
711  real(DP) :: const2
712  !
713  ! -- initialize
714  tled = done / delt
715  !
716  ! -- Calculate sorption change
717  do n = 1, nodes
718  !
719  ! -- initialize rates
720  this%ratesrb(n) = dzero
721  !
722  ! -- skip if transport inactive
723  if (this%ibound(n) <= 0) cycle
724  !
725  ! -- assign variables
726  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
727  swtpdt = this%fmi%gwfsat(n)
728  swt = this%fmi%gwfsatold(n, delt)
729  volfracm = this%get_volfracm(n)
730  rhobm = this%bulk_density(n)
731  const1 = this%distcoef(n)
732  const2 = 0.
733  if (this%isrb == sorption_freund .or. this%isrb == sorption_lang) then
734  const2 = this%sp2(n)
735  end if
736  call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
737  cold(n), swtpdt, swt, const1, const2, &
738  rate=rate)
739  this%ratesrb(n) = rate
740  idiag = this%dis%con%ia(n)
741  flowja(idiag) = flowja(idiag) + rate
742  !
743  end do
744  end subroutine mst_cq_srb
745 
746  !> @ brief Calculate decay-sorption terms for package
747  !!
748  !! Method to calculate decay-sorption terms for the package.
749  !<
750  subroutine mst_cq_dcy_srb(this, nodes, cnew, cold, flowja)
751  ! -- modules
752  use tdismodule, only: delt
753  ! -- dummy
754  class(gwtmsttype) :: this !< GwtMstType object
755  integer(I4B), intent(in) :: nodes !< number of nodes
756  real(DP), intent(in), dimension(nodes) :: cnew !< concentration at end of this time step
757  real(DP), intent(in), dimension(nodes) :: cold !< concentration at end of last time step
758  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes
759  ! -- local
760  integer(I4B) :: n
761  integer(I4B) :: idiag
762  real(DP) :: rate
763  real(DP) :: hhcof, rrhs
764  real(DP) :: vcell
765  real(DP) :: swnew
766  real(DP) :: distcoef
767  real(DP) :: volfracm
768  real(DP) :: rhobm
769  real(DP) :: term
770  real(DP) :: csrb
771  real(DP) :: csrbnew
772  real(DP) :: csrbold
773  real(DP) :: decay_rate
774  !
775  ! -- Calculate sorbed decay change
776  ! This routine will only be called if sorption and decay are active
777  do n = 1, nodes
778  !
779  ! -- initialize rates
780  this%ratedcys(n) = dzero
781  !
782  ! -- skip if transport inactive
783  if (this%ibound(n) <= 0) cycle
784  !
785  ! -- set variables
786  hhcof = dzero
787  rrhs = dzero
788  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
789  swnew = this%fmi%gwfsat(n)
790  distcoef = this%distcoef(n)
791  volfracm = this%get_volfracm(n)
792  rhobm = this%bulk_density(n)
793  term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
794  !
795  ! -- add sorbed mass decay rate terms to accumulators
796  select case (this%idcy)
797  case (decay_first_order)
798  !
799  select case (this%isrb)
800  case (sorption_linear)
801  !
802  ! -- first order decay rate is a function of concentration, so add
803  ! to left hand side
804  hhcof = -term * distcoef
805  case (sorption_freund)
806  !
807  ! -- nonlinear Freundlich sorption, so add to RHS
808  csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
809  rrhs = term * csrb
810  case (sorption_lang)
811  !
812  ! -- nonlinear Lanmuir sorption, so add to RHS
813  csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
814  rrhs = term * csrb
815  end select
816  case (decay_zero_order)
817  !
818  ! -- Call function to get zero-order decay rate, which may be changed
819  ! from the user-specified rate to prevent negative concentrations
820  if (distcoef > dzero) then
821  select case (this%isrb)
822  case (sorption_linear)
823  csrbold = cold(n) * distcoef
824  csrbnew = cnew(n) * distcoef
825  case (sorption_freund)
826  csrbold = get_freundlich_conc(cold(n), distcoef, this%sp2(n))
827  csrbnew = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
828  case (sorption_lang)
829  csrbold = get_langmuir_conc(cold(n), distcoef, this%sp2(n))
830  csrbnew = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
831  end select
832  decay_rate = get_zero_order_decay(this%decay_sorbed(n), &
833  this%decayslast(n), &
834  0, csrbold, csrbnew, delt)
835  rrhs = decay_rate * volfracm * rhobm * swnew * vcell
836  end if
837  end select
838  !
839  ! -- calculate rate
840  rate = hhcof * cnew(n) - rrhs
841  this%ratedcys(n) = rate
842  idiag = this%dis%con%ia(n)
843  flowja(idiag) = flowja(idiag) + rate
844  !
845  end do
846  end subroutine mst_cq_dcy_srb
847 
848  !> @ brief Calculate sorbed concentration
849  !<
850  subroutine mst_calc_csrb(this, cnew)
851  ! -- dummy
852  class(gwtmsttype) :: this !< GwtMstType object
853  real(DP), intent(in), dimension(:) :: cnew !< concentration at end of this time step
854  ! -- local
855  integer(I4B) :: n
856  real(DP) :: distcoef
857  real(DP) :: csrb
858 
859  ! Calculate sorbed concentration
860  do n = 1, size(cnew)
861  csrb = dzero
862  if (this%ibound(n) > 0) then
863  distcoef = this%distcoef(n)
864  select case (this%isrb)
865  case (sorption_linear)
866  csrb = cnew(n) * distcoef
867  case (sorption_freund)
868  csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
869  case (sorption_lang)
870  csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
871  end select
872  end if
873  this%csrb(n) = csrb
874  end do
875 
876  end subroutine mst_calc_csrb
877 
878  !> @ brief Calculate budget terms for package
879  !<
880  subroutine mst_bd(this, isuppress_output, model_budget)
881  ! -- modules
882  use tdismodule, only: delt
884  ! -- dummy
885  class(gwtmsttype) :: this !< GwtMstType object
886  integer(I4B), intent(in) :: isuppress_output !< flag to suppress output
887  type(budgettype), intent(inout) :: model_budget !< model budget object
888  ! -- local
889  real(DP) :: rin
890  real(DP) :: rout
891  !
892  ! -- sto
893  call rate_accumulator(this%ratesto, rin, rout)
894  call model_budget%addentry(rin, rout, delt, budtxt(1), &
895  isuppress_output, rowlabel=this%packName)
896  !
897  ! -- dcy
898  if (this%idcy /= decay_off) then
899  call rate_accumulator(this%ratedcy, rin, rout)
900  call model_budget%addentry(rin, rout, delt, budtxt(2), &
901  isuppress_output, rowlabel=this%packName)
902  end if
903  !
904  ! -- srb
905  if (this%isrb /= sorption_off) then
906  call rate_accumulator(this%ratesrb, rin, rout)
907  call model_budget%addentry(rin, rout, delt, budtxt(3), &
908  isuppress_output, rowlabel=this%packName)
909  end if
910  !
911  ! -- srb dcy
912  if (this%isrb /= sorption_off .and. this%idcy /= decay_off) then
913  call rate_accumulator(this%ratedcys, rin, rout)
914  call model_budget%addentry(rin, rout, delt, budtxt(4), &
915  isuppress_output, rowlabel=this%packName)
916  end if
917  end subroutine mst_bd
918 
919  !> @ brief Output flow terms for package
920  !!
921  !! Method to output terms for the package.
922  !<
923  subroutine mst_ot_flow(this, icbcfl, icbcun)
924  ! -- dummy
925  class(gwtmsttype) :: this !< GwtMstType object
926  integer(I4B), intent(in) :: icbcfl !< flag and unit number for cell-by-cell output
927  integer(I4B), intent(in) :: icbcun !< flag indication if cell-by-cell data should be saved
928  ! -- local
929  integer(I4B) :: ibinun
930  integer(I4B) :: iprint, nvaluesp, nwidthp
931  character(len=1) :: cdatafmp = ' ', editdesc = ' '
932  real(DP) :: dinact
933  !
934  ! -- Set unit number for binary output
935  if (this%ipakcb < 0) then
936  ibinun = icbcun
937  elseif (this%ipakcb == 0) then
938  ibinun = 0
939  else
940  ibinun = this%ipakcb
941  end if
942  if (icbcfl == 0) ibinun = 0
943  !
944  ! -- Record the storage rate if requested
945  if (ibinun /= 0) then
946  iprint = 0
947  dinact = dzero
948  !
949  ! -- sto
950  call this%dis%record_array(this%ratesto, this%iout, iprint, -ibinun, &
951  budtxt(1), cdatafmp, nvaluesp, &
952  nwidthp, editdesc, dinact)
953  !
954  ! -- dcy
955  if (this%idcy /= decay_off) &
956  call this%dis%record_array(this%ratedcy, this%iout, iprint, -ibinun, &
957  budtxt(2), cdatafmp, nvaluesp, &
958  nwidthp, editdesc, dinact)
959  !
960  ! -- srb
961  if (this%isrb /= sorption_off) &
962  call this%dis%record_array(this%ratesrb, this%iout, iprint, -ibinun, &
963  budtxt(3), cdatafmp, nvaluesp, &
964  nwidthp, editdesc, dinact)
965  !
966  ! -- dcy srb
967  if (this%isrb /= sorption_off .and. this%idcy /= decay_off) &
968  call this%dis%record_array(this%ratedcys, this%iout, iprint, -ibinun, &
969  budtxt(4), cdatafmp, nvaluesp, &
970  nwidthp, editdesc, dinact)
971  end if
972  end subroutine mst_ot_flow
973 
974  !> @brief Save sorbate concentration array to binary file
975  !<
976  subroutine mst_ot_dv(this, idvsave)
977  ! -- dummy
978  class(gwtmsttype) :: this
979  integer(I4B), intent(in) :: idvsave
980  ! -- local
981  character(len=1) :: cdatafmp = ' ', editdesc = ' '
982  integer(I4B) :: ibinun
983  integer(I4B) :: iprint
984  integer(I4B) :: nvaluesp
985  integer(I4B) :: nwidthp
986  real(DP) :: dinact
987 
988  ! Set unit number for sorbate output
989  if (this%ioutsorbate /= 0) then
990  ibinun = 1
991  else
992  ibinun = 0
993  end if
994  if (idvsave == 0) ibinun = 0
995 
996  ! save sorbate concentration array
997  if (ibinun /= 0) then
998  iprint = 0
999  dinact = dhnoflo
1000  if (this%ioutsorbate /= 0) then
1001  ibinun = this%ioutsorbate
1002  call this%dis%record_array(this%csrb, this%iout, iprint, ibinun, &
1003  ' SORBATE', cdatafmp, nvaluesp, &
1004  nwidthp, editdesc, dinact)
1005  end if
1006  end if
1007 
1008  end subroutine mst_ot_dv
1009 
1010  !> @ brief Deallocate
1011  !!
1012  !! Method to deallocate memory for the package.
1013  !<
1014  subroutine mst_da(this)
1015  ! -- modules
1017  ! -- dummy
1018  class(gwtmsttype) :: this !< GwtMstType object
1019  !
1020  ! -- Deallocate arrays if package was active
1021  if (this%inunit > 0) then
1022  call mem_deallocate(this%porosity)
1023  call mem_deallocate(this%thetam)
1024  call mem_deallocate(this%volfracim)
1025  call mem_deallocate(this%ratesto)
1026  call mem_deallocate(this%idcy)
1027  call mem_deallocate(this%decay)
1028  call mem_deallocate(this%decay_sorbed)
1029  call mem_deallocate(this%ratedcy)
1030  call mem_deallocate(this%decaylast)
1031  call mem_deallocate(this%decayslast)
1032  call mem_deallocate(this%isrb)
1033  call mem_deallocate(this%ioutsorbate)
1034  call mem_deallocate(this%bulk_density)
1035  call mem_deallocate(this%distcoef)
1036  call mem_deallocate(this%sp2)
1037  call mem_deallocate(this%ratesrb)
1038  call mem_deallocate(this%csrb)
1039  call mem_deallocate(this%ratedcys)
1040  this%ibound => null()
1041  this%fmi => null()
1042  end if
1043  !
1044  ! -- Scalars
1045  !
1046  ! -- deallocate parent
1047  call this%NumericalPackageType%da()
1048  end subroutine mst_da
1049 
1050  !> @ brief Allocate scalar variables for package
1051  !!
1052  !! Method to allocate scalar variables for the package.
1053  !<
1054  subroutine allocate_scalars(this)
1055  ! -- modules
1057  ! -- dummy
1058  class(gwtmsttype) :: this !< GwtMstType object
1059  !
1060  ! -- Allocate scalars in NumericalPackageType
1061  call this%NumericalPackageType%allocate_scalars()
1062  !
1063  ! -- Allocate
1064  call mem_allocate(this%isrb, 'ISRB', this%memoryPath)
1065  call mem_allocate(this%ioutsorbate, 'IOUTSORBATE', this%memoryPath)
1066  call mem_allocate(this%idcy, 'IDCY', this%memoryPath)
1067  !
1068  ! -- Initialize
1069  this%isrb = izero
1070  this%ioutsorbate = 0
1071  this%idcy = izero
1072  end subroutine allocate_scalars
1073 
1074  !> @ brief Allocate arrays for package
1075  !!
1076  !! Method to allocate arrays for the package.
1077  !<
1078  subroutine allocate_arrays(this, nodes)
1079  ! -- modules
1081  use constantsmodule, only: dzero
1082  ! -- dummy
1083  class(gwtmsttype) :: this !< GwtMstType object
1084  integer(I4B), intent(in) :: nodes !< number of nodes
1085  ! -- local
1086  integer(I4B) :: n
1087  !
1088  ! -- Allocate
1089  ! -- sto
1090  call mem_allocate(this%porosity, nodes, 'POROSITY', this%memoryPath)
1091  call mem_allocate(this%thetam, nodes, 'THETAM', this%memoryPath)
1092  call mem_allocate(this%volfracim, nodes, 'VOLFRACIM', this%memoryPath)
1093  call mem_allocate(this%ratesto, nodes, 'RATESTO', this%memoryPath)
1094  !
1095  ! -- dcy
1096  if (this%idcy == decay_off) then
1097  call mem_allocate(this%ratedcy, 1, 'RATEDCY', this%memoryPath)
1098  call mem_allocate(this%decay, 1, 'DECAY', this%memoryPath)
1099  call mem_allocate(this%decaylast, 1, 'DECAYLAST', this%memoryPath)
1100  else
1101  call mem_allocate(this%ratedcy, this%dis%nodes, 'RATEDCY', this%memoryPath)
1102  call mem_allocate(this%decay, nodes, 'DECAY', this%memoryPath)
1103  call mem_allocate(this%decaylast, nodes, 'DECAYLAST', this%memoryPath)
1104  end if
1105  if (this%idcy /= decay_off .and. this%isrb /= sorption_off) then
1106  call mem_allocate(this%ratedcys, this%dis%nodes, 'RATEDCYS', &
1107  this%memoryPath)
1108  call mem_allocate(this%decayslast, this%dis%nodes, 'DECAYSLAST', &
1109  this%memoryPath)
1110  else
1111  call mem_allocate(this%ratedcys, 1, 'RATEDCYS', this%memoryPath)
1112  call mem_allocate(this%decayslast, 1, 'DECAYSLAST', this%memoryPath)
1113  end if
1114  call mem_allocate(this%decay_sorbed, 1, 'DECAY_SORBED', &
1115  this%memoryPath)
1116  !
1117  ! -- srb
1118  if (this%isrb == sorption_off) then
1119  call mem_allocate(this%bulk_density, 1, 'BULK_DENSITY', this%memoryPath)
1120  call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
1121  call mem_allocate(this%distcoef, 1, 'DISTCOEF', this%memoryPath)
1122  call mem_allocate(this%ratesrb, 1, 'RATESRB', this%memoryPath)
1123  call mem_allocate(this%csrb, 1, 'CSRB', this%memoryPath)
1124  else
1125  call mem_allocate(this%bulk_density, nodes, 'BULK_DENSITY', this%memoryPath)
1126  call mem_allocate(this%distcoef, nodes, 'DISTCOEF', this%memoryPath)
1127  call mem_allocate(this%ratesrb, nodes, 'RATESRB', this%memoryPath)
1128  call mem_allocate(this%csrb, nodes, 'CSRB', this%memoryPath)
1129  if (this%isrb == sorption_linear) then
1130  call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
1131  else
1132  call mem_allocate(this%sp2, nodes, 'SP2', this%memoryPath)
1133  end if
1134  end if
1135  !
1136  ! -- Initialize
1137  do n = 1, nodes
1138  this%porosity(n) = dzero
1139  this%thetam(n) = dzero
1140  this%volfracim(n) = dzero
1141  this%ratesto(n) = dzero
1142  end do
1143  do n = 1, size(this%decay)
1144  this%decay(n) = dzero
1145  this%ratedcy(n) = dzero
1146  this%decaylast(n) = dzero
1147  end do
1148  do n = 1, size(this%bulk_density)
1149  this%bulk_density(n) = dzero
1150  this%distcoef(n) = dzero
1151  this%ratesrb(n) = dzero
1152  this%csrb(n) = dzero
1153  end do
1154  do n = 1, size(this%sp2)
1155  this%sp2(n) = dzero
1156  end do
1157  do n = 1, size(this%ratedcys)
1158  this%ratedcys(n) = dzero
1159  this%decayslast(n) = dzero
1160  end do
1161  end subroutine allocate_arrays
1162 
1163  !> @ brief Read options for package
1164  !!
1165  !! Method to read options for the package.
1166  !<
1167  subroutine read_options(this)
1168  ! -- modules
1169  use openspecmodule, only: access, form
1170  use inputoutputmodule, only: getunit, openfile
1171  ! -- dummy
1172  class(gwtmsttype) :: this !< GwtMstType object
1173  ! -- local
1174  character(len=LINELENGTH) :: keyword
1175  character(len=LINELENGTH) :: sorption
1176  character(len=MAXCHARLEN) :: fname
1177  integer(I4B) :: ierr
1178  logical :: isfound, endOfBlock
1179  ! -- formats
1180  character(len=*), parameter :: fmtisvflow = &
1181  "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
1182  &WHENEVER ICBCFL IS NOT ZERO.')"
1183  character(len=*), parameter :: fmtlinear = &
1184  &"(4x,'LINEAR SORPTION IS ACTIVE. ')"
1185  character(len=*), parameter :: fmtfreundlich = &
1186  &"(4x,'FREUNDLICH SORPTION IS ACTIVE. ')"
1187  character(len=*), parameter :: fmtlangmuir = &
1188  &"(4x,'LANGMUIR SORPTION IS ACTIVE. ')"
1189  character(len=*), parameter :: fmtidcy1 = &
1190  &"(4x,'FIRST-ORDER DECAY IS ACTIVE. ')"
1191  character(len=*), parameter :: fmtidcy2 = &
1192  &"(4x,'ZERO-ORDER DECAY IS ACTIVE. ')"
1193  character(len=*), parameter :: fmtfileout = &
1194  "(4x,'MST ',1x,a,1x,' WILL BE SAVED TO FILE: ',a,/4x,&
1195  &'OPENED ON UNIT: ',I7)"
1196  !
1197  ! -- get options block
1198  call this%parser%GetBlock('OPTIONS', isfound, ierr, blockrequired=.false., &
1199  supportopenclose=.true.)
1200  !
1201  ! -- parse options block if detected
1202  if (isfound) then
1203  write (this%iout, '(1x,a)') 'PROCESSING MOBILE STORAGE AND TRANSFER OPTIONS'
1204  do
1205  call this%parser%GetNextLine(endofblock)
1206  if (endofblock) exit
1207  call this%parser%GetStringCaps(keyword)
1208  select case (keyword)
1209  case ('SAVE_FLOWS')
1210  this%ipakcb = -1
1211  write (this%iout, fmtisvflow)
1212  case ('SORBTION')
1213  call store_error('SORBTION is not a valid option. Use &
1214  &SORPTION instead.')
1215  call this%parser%StoreErrorUnit()
1216  case ('SORPTION')
1217  call this%parser%GetStringCaps(sorption)
1218  select case (sorption)
1219  case ('LINEAR', '')
1220  this%isrb = sorption_linear
1221  write (this%iout, fmtlinear)
1222  case ('FREUNDLICH')
1223  this%isrb = sorption_freund
1224  write (this%iout, fmtfreundlich)
1225  case ('LANGMUIR')
1226  this%isrb = sorption_lang
1227  write (this%iout, fmtlangmuir)
1228  case default
1229  call store_error('Unknown sorption type was specified ' &
1230  & //'('//trim(adjustl(sorption))//').'// &
1231  &' Sorption must be specified as LINEAR, &
1232  &FREUNDLICH, or LANGMUIR.')
1233  call this%parser%StoreErrorUnit()
1234  end select
1235  case ('FIRST_ORDER_DECAY')
1236  this%idcy = decay_first_order
1237  write (this%iout, fmtidcy1)
1238  case ('ZERO_ORDER_DECAY')
1239  this%idcy = decay_zero_order
1240  write (this%iout, fmtidcy2)
1241  case ('SORBATE')
1242  call this%parser%GetStringCaps(keyword)
1243  if (keyword == 'FILEOUT') then
1244  call this%parser%GetString(fname)
1245  this%ioutsorbate = getunit()
1246  call openfile(this%ioutsorbate, this%iout, fname, 'DATA(BINARY)', &
1247  form, access, 'REPLACE', mode_opt=mnormal)
1248  write (this%iout, fmtfileout) &
1249  'SORBATE', fname, this%ioutsorbate
1250  else
1251  errmsg = 'Optional SORBATE keyword must be '// &
1252  'followed by FILEOUT.'
1253  call store_error(errmsg)
1254  end if
1255  case default
1256  write (errmsg, '(a,a)') 'Unknown MST option: ', trim(keyword)
1257  call store_error(errmsg)
1258  call this%parser%StoreErrorUnit()
1259  end select
1260  end do
1261  write (this%iout, '(1x,a)') 'END OF MOBILE STORAGE AND TRANSFER OPTIONS'
1262  end if
1263  end subroutine read_options
1264 
1265  !> @ brief Read data for package
1266  !!
1267  !! Method to read data for the package.
1268  !<
1269  subroutine read_data(this)
1270  ! -- modules
1271  use constantsmodule, only: linelength
1273  ! -- dummy
1274  class(gwtmsttype) :: this !< GwtMstType object
1275  ! -- local
1276  character(len=LINELENGTH) :: keyword
1277  character(len=:), allocatable :: line
1278  integer(I4B) :: istart, istop, lloc, ierr, n
1279  logical :: isfound, endOfBlock
1280  logical, dimension(6) :: lname
1281  character(len=24), dimension(6) :: aname
1282  ! -- data
1283  data aname(1)/' MOBILE DOMAIN POROSITY'/
1284  data aname(2)/' BULK DENSITY'/
1285  data aname(3)/'DISTRIBUTION COEFFICIENT'/
1286  data aname(4)/' DECAY RATE'/
1287  data aname(5)/' DECAY SORBED RATE'/
1288  data aname(6)/' SECOND SORPTION PARAM'/
1289  !
1290  ! -- initialize
1291  isfound = .false.
1292  lname(:) = .false.
1293  !
1294  ! -- get griddata block
1295  call this%parser%GetBlock('GRIDDATA', isfound, ierr)
1296  if (isfound) then
1297  write (this%iout, '(1x,a)') 'PROCESSING GRIDDATA'
1298  do
1299  call this%parser%GetNextLine(endofblock)
1300  if (endofblock) exit
1301  call this%parser%GetStringCaps(keyword)
1302  call this%parser%GetRemainingLine(line)
1303  lloc = 1
1304  select case (keyword)
1305  case ('POROSITY')
1306  call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1307  this%parser%iuactive, this%porosity, &
1308  aname(1))
1309  lname(1) = .true.
1310  case ('BULK_DENSITY')
1311  if (this%isrb == sorption_off) &
1312  call mem_reallocate(this%bulk_density, this%dis%nodes, &
1313  'BULK_DENSITY', trim(this%memoryPath))
1314  call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1315  this%parser%iuactive, &
1316  this%bulk_density, aname(2))
1317  lname(2) = .true.
1318  case ('DISTCOEF')
1319  if (this%isrb == sorption_off) &
1320  call mem_reallocate(this%distcoef, this%dis%nodes, 'DISTCOEF', &
1321  trim(this%memoryPath))
1322  call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1323  this%parser%iuactive, this%distcoef, &
1324  aname(3))
1325  lname(3) = .true.
1326  case ('DECAY')
1327  if (this%idcy == decay_off) &
1328  call mem_reallocate(this%decay, this%dis%nodes, 'DECAY', &
1329  trim(this%memoryPath))
1330  call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1331  this%parser%iuactive, this%decay, &
1332  aname(4))
1333  lname(4) = .true.
1334  case ('DECAY_SORBED')
1335  call mem_reallocate(this%decay_sorbed, this%dis%nodes, &
1336  'DECAY_SORBED', trim(this%memoryPath))
1337  call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1338  this%parser%iuactive, &
1339  this%decay_sorbed, aname(5))
1340  lname(5) = .true.
1341  case ('SP2')
1342  if (this%isrb == sorption_off .or. this%isrb == sorption_linear) &
1343  call mem_reallocate(this%sp2, this%dis%nodes, 'SP2', &
1344  trim(this%memoryPath))
1345  call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1346  this%parser%iuactive, this%sp2, &
1347  aname(6))
1348  lname(6) = .true.
1349  case default
1350  write (errmsg, '(a,a)') 'Unknown GRIDDATA tag: ', trim(keyword)
1351  call store_error(errmsg)
1352  call this%parser%StoreErrorUnit()
1353  end select
1354  end do
1355  write (this%iout, '(1x,a)') 'END PROCESSING GRIDDATA'
1356  else
1357  write (errmsg, '(a)') 'Required GRIDDATA block not found.'
1358  call store_error(errmsg)
1359  call this%parser%StoreErrorUnit()
1360  end if
1361  !
1362  ! -- Check for required porosity
1363  if (.not. lname(1)) then
1364  write (errmsg, '(a)') 'POROSITY not specified in GRIDDATA block.'
1365  call store_error(errmsg)
1366  end if
1367  !
1368  ! -- Check for required sorption variables
1369  if (this%isrb == sorption_linear .or. this%isrb == sorption_freund .or. &
1370  this%isrb == sorption_lang) then
1371  if (.not. lname(2)) then
1372  write (errmsg, '(a)') 'Sorption is active but BULK_DENSITY &
1373  &not specified. BULK_DENSITY must be specified in GRIDDATA block.'
1374  call store_error(errmsg)
1375  end if
1376  if (.not. lname(3)) then
1377  write (errmsg, '(a)') 'Sorption is active but distribution &
1378  &coefficient not specified. DISTCOEF must be specified in &
1379  &GRIDDATA block.'
1380  call store_error(errmsg)
1381  end if
1382  if (this%isrb == sorption_freund .or. this%isrb == sorption_lang) then
1383  if (.not. lname(6)) then
1384  write (errmsg, '(a)') 'Freundlich or langmuir sorption is active &
1385  &but SP2 not specified. SP2 must be specified in &
1386  &GRIDDATA block.'
1387  call store_error(errmsg)
1388  end if
1389  end if
1390  else
1391  if (lname(2)) then
1392  write (warnmsg, '(a)') 'Sorption is not active but &
1393  &BULK_DENSITY was specified. BULK_DENSITY will have no affect on &
1394  &simulation results.'
1395  call store_warning(warnmsg)
1396  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1397  end if
1398  if (lname(3)) then
1399  write (warnmsg, '(a)') 'Sorption is not active but &
1400  &distribution coefficient was specified. DISTCOEF will have &
1401  &no affect on simulation results.'
1402  call store_warning(warnmsg)
1403  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1404  end if
1405  if (lname(6)) then
1406  write (warnmsg, '(a)') 'Sorption is not active but &
1407  &SP2 was specified. SP2 will have &
1408  &no affect on simulation results.'
1409  call store_warning(warnmsg)
1410  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1411  end if
1412  end if
1413  !
1414  ! -- Check for required decay/production rate coefficients
1415  if (this%idcy /= decay_off) then
1416  if (.not. lname(4)) then
1417  write (errmsg, '(a)') 'First or zero order decay is &
1418  &active but the first rate coefficient is not specified. DECAY &
1419  &must be specified in GRIDDATA block.'
1420  call store_error(errmsg)
1421  end if
1422  if (.not. lname(5)) then
1423  !
1424  ! -- If DECAY_SORBED not specified and sorption is active, then
1425  ! terminate with an error
1426  if (this%isrb == sorption_linear .or. this%isrb == sorption_freund .or. &
1427  this%isrb == sorption_lang) then
1428  write (errmsg, '(a)') 'DECAY_SORBED not provided in GRIDDATA &
1429  &block but decay and sorption are active. Specify DECAY_SORBED &
1430  &in GRIDDATA block.'
1431  call store_error(errmsg)
1432  end if
1433  end if
1434  else
1435  if (lname(4)) then
1436  write (warnmsg, '(a)') 'First- or zero-order decay &
1437  &is not active but decay was specified. DECAY will &
1438  &have no affect on simulation results.'
1439  call store_warning(warnmsg)
1440  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1441  end if
1442  if (lname(5)) then
1443  write (warnmsg, '(a)') 'First- or zero-order decay &
1444  &is not active but DECAY_SORBED was specified. &
1445  &DECAY_SORBED will have no affect on simulation results.'
1446  call store_warning(warnmsg)
1447  write (this%iout, '(1x,a)') 'WARNING. '//warnmsg
1448  end if
1449  end if
1450  !
1451  ! -- terminate if errors
1452  if (count_errors() > 0) then
1453  call this%parser%StoreErrorUnit()
1454  end if
1455  !
1456  ! -- initialize thetam from porosity
1457  do n = 1, size(this%porosity)
1458  this%thetam(n) = this%porosity(n)
1459  end do
1460  end subroutine read_data
1461 
1462  !> @ brief Add volfrac values to volfracim
1463  !!
1464  !! Method to add immobile domain volume fracions, which are stored as a
1465  !! cumulative value in volfracim.
1466  !<
1467  subroutine addto_volfracim(this, volfracim)
1468  ! -- dummy
1469  class(gwtmsttype) :: this !< GwtMstType object
1470  real(DP), dimension(:), intent(in) :: volfracim !< immobile domain volume fraction that contributes to total immobile volume fraction
1471  ! -- local
1472  integer(I4B) :: n
1473  !
1474  ! -- Add to volfracim
1475  do n = 1, this%dis%nodes
1476  this%volfracim(n) = this%volfracim(n) + volfracim(n)
1477  end do
1478  !
1479  ! -- An immobile domain is adding a volume fraction, so update thetam
1480  ! accordingly.
1481  do n = 1, this%dis%nodes
1482  this%thetam(n) = this%get_volfracm(n) * this%porosity(n)
1483  end do
1484  end subroutine addto_volfracim
1485 
1486  !> @ brief Return mobile domain volume fraction
1487  !!
1488  !! Calculate and return the volume fraction of the aquifer that is mobile
1489  !<
1490  function get_volfracm(this, node) result(volfracm)
1491  ! -- dummy
1492  class(gwtmsttype) :: this !< GwtMstType object
1493  integer(I4B), intent(in) :: node !< node number
1494  ! -- return
1495  real(dp) :: volfracm
1496  !
1497  volfracm = done - this%volfracim(node)
1498  end function get_volfracm
1499 
1500  !> @ brief Calculate sorption concentration using Freundlich
1501  !!
1502  !! Function to calculate sorption concentration using Freundlich
1503  !<
1504  function get_freundlich_conc(conc, kf, a) result(cbar)
1505  ! -- dummy
1506  real(dp), intent(in) :: conc !< solute concentration
1507  real(dp), intent(in) :: kf !< freundlich constant
1508  real(dp), intent(in) :: a !< freundlich exponent
1509  ! -- return
1510  real(dp) :: cbar
1511  !
1512  if (conc > dzero) then
1513  cbar = kf * conc**a
1514  else
1515  cbar = dzero
1516  end if
1517  end function
1518 
1519  !> @ brief Calculate sorption concentration using Langmuir
1520  !!
1521  !! Function to calculate sorption concentration using Langmuir
1522  !<
1523  function get_langmuir_conc(conc, kl, sbar) result(cbar)
1524  ! -- dummy
1525  real(dp), intent(in) :: conc !< solute concentration
1526  real(dp), intent(in) :: kl !< langmuir constant
1527  real(dp), intent(in) :: sbar !< langmuir sorption sites
1528  ! -- return
1529  real(dp) :: cbar
1530  !
1531  if (conc > dzero) then
1532  cbar = (kl * sbar * conc) / (done + kl * conc)
1533  else
1534  cbar = dzero
1535  end if
1536  end function
1537 
1538  !> @ brief Calculate sorption derivative using Freundlich
1539  !!
1540  !! Function to calculate sorption derivative using Freundlich
1541  !<
1542  function get_freundlich_derivative(conc, kf, a) result(derv)
1543  ! -- dummy
1544  real(dp), intent(in) :: conc !< solute concentration
1545  real(dp), intent(in) :: kf !< freundlich constant
1546  real(dp), intent(in) :: a !< freundlich exponent
1547  ! -- return
1548  real(dp) :: derv
1549  !
1550  if (conc > dzero) then
1551  derv = kf * a * conc**(a - done)
1552  else
1553  derv = dzero
1554  end if
1555  end function
1556 
1557  !> @ brief Calculate sorption derivative using Langmuir
1558  !!
1559  !! Function to calculate sorption derivative using Langmuir
1560  !<
1561  function get_langmuir_derivative(conc, kl, sbar) result(derv)
1562  ! -- dummy
1563  real(dp), intent(in) :: conc !< solute concentration
1564  real(dp), intent(in) :: kl !< langmuir constant
1565  real(dp), intent(in) :: sbar !< langmuir sorption sites
1566  ! -- return
1567  real(dp) :: derv
1568  !
1569  if (conc > dzero) then
1570  derv = (kl * sbar) / (done + kl * conc)**dtwo
1571  else
1572  derv = dzero
1573  end if
1574  end function
1575 
1576  !> @ brief Get effective Freundlich distribution coefficient
1577  !<
1578  function get_freundlich_kd(conc, kf, a) result(kd)
1579  ! -- dummy
1580  real(dp), intent(in) :: conc !< solute concentration
1581  real(dp), intent(in) :: kf !< freundlich constant
1582  real(dp), intent(in) :: a !< freundlich exponent
1583  ! -- return
1584  real(dp) :: kd !< effective distribution coefficient
1585  !
1586  if (conc > dzero) then
1587  kd = kf * conc**(a - done)
1588  else
1589  kd = dzero
1590  end if
1591  end function get_freundlich_kd
1592 
1593  !> @ brief Get effective Langmuir distribution coefficient
1594  !<
1595  function get_langmuir_kd(conc, kl, sbar) result(kd)
1596  ! -- dummy
1597  real(dp), intent(in) :: conc !< solute concentration
1598  real(dp), intent(in) :: kl !< langmuir constant
1599  real(dp), intent(in) :: sbar !< langmuir sorption sites
1600  ! -- return
1601  real(dp) :: kd !< effective distribution coefficient
1602  !
1603  if (conc > dzero) then
1604  kd = (kl * sbar) / (done + kl * conc)
1605  else
1606  kd = dzero
1607  end if
1608  end function get_langmuir_kd
1609 
1610  !> @ brief Calculate zero-order decay rate and constrain if necessary
1611  !!
1612  !! Function to calculate the zero-order decay rate from the user specified
1613  !! decay rate. If the decay rate is positive, then the decay rate must
1614  !! be constrained so that more mass is not removed than is available.
1615  !! Without this constraint, negative concentrations could result from
1616  !! zero-order decay.
1617  !<
1618  function get_zero_order_decay(decay_rate_usr, decay_rate_last, kiter, &
1619  cold, cnew, delt) result(decay_rate)
1620  ! -- dummy
1621  real(dp), intent(in) :: decay_rate_usr !< user-entered decay rate
1622  real(dp), intent(in) :: decay_rate_last !< decay rate used for last iteration
1623  integer(I4B), intent(in) :: kiter !< Picard iteration counter
1624  real(dp), intent(in) :: cold !< concentration at end of last time step
1625  real(dp), intent(in) :: cnew !< concentration at end of this time step
1626  real(dp), intent(in) :: delt !< length of time step
1627  ! -- return
1628  real(dp) :: decay_rate !< returned value for decay rate
1629  !
1630  ! -- Return user rate if production, otherwise constrain, if necessary
1631  if (decay_rate_usr < dzero) then
1632  !
1633  ! -- Production, no need to limit rate
1634  decay_rate = decay_rate_usr
1635  else
1636  !
1637  ! -- Need to ensure decay does not result in negative
1638  ! concentration, so reduce the rate if it would result in
1639  ! removing more mass than is in the cell.
1640  if (kiter == 1) then
1641  decay_rate = min(decay_rate_usr, cold / delt)
1642  else
1643  decay_rate = decay_rate_last
1644  if (cnew < dzero) then
1645  decay_rate = decay_rate_last + cnew / delt
1646  else if (cnew > cold) then
1647  decay_rate = decay_rate_last + cnew / delt
1648  end if
1649  decay_rate = min(decay_rate_usr, decay_rate)
1650  end if
1651  decay_rate = max(decay_rate, dzero)
1652  end if
1653  end function get_zero_order_decay
1654 
1655 end module gwtmstmodule
This module contains the BudgetModule.
Definition: Budget.f90:20
subroutine, public rate_accumulator(flow, rin, rout)
@ brief Rate accumulator subroutine
Definition: Budget.f90:632
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
@ mnormal
normal output mode
Definition: Constants.f90:206
real(dp), parameter dhnoflo
real no flow constant
Definition: Constants.f90:93
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
integer(i4b), parameter izero
integer constant zero
Definition: Constants.f90:51
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter maxcharlen
maximum length of char string
Definition: Constants.f90:47
real(dp), parameter dtwo
real constant 2
Definition: Constants.f90:79
integer(i4b), parameter lenbudtxt
maximum length of a budget component names
Definition: Constants.f90:37
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
– @ brief Mobile Storage and Transfer (MST) Module
Definition: gwt-mst.f90:10
subroutine mst_srb_term(isrb, volfracm, rhobm, vcell, tled, cnew, cold, swnew, swold, const1, const2, rate, hcofval, rhsval)
@ brief Calculate sorption terms
Definition: gwt-mst.f90:381
subroutine mst_cq_srb(this, nodes, cnew, cold, flowja)
@ brief Calculate sorption terms for package
Definition: gwt-mst.f90:693
subroutine mst_calc_csrb(this, cnew)
@ brief Calculate sorbed concentration
Definition: gwt-mst.f90:851
subroutine addto_volfracim(this, volfracim)
@ brief Add volfrac values to volfracim
Definition: gwt-mst.f90:1468
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
Definition: gwt-mst.f90:1620
character(len=lenbudtxt), dimension(nbditems) budtxt
Definition: gwt-mst.f90:28
subroutine mst_cq(this, nodes, cnew, cold, flowja)
@ brief Calculate flows for package
Definition: gwt-mst.f90:555
@ decay_zero_order
Zeroth-order decay.
Definition: gwt-mst.f90:37
@ decay_first_order
First-order decay.
Definition: gwt-mst.f90:36
@ sorption_freund
Freundlich sorption between aqueous and solid phases.
Definition: gwt-mst.f90:40
@ sorption_lang
Langmuir sorption between aqueous and solid phases.
Definition: gwt-mst.f90:41
@ sorption_linear
Linear sorption between aqueous and solid phases.
Definition: gwt-mst.f90:39
@ decay_off
Decay (or production) of mass inactive (default)
Definition: gwt-mst.f90:35
@ sorption_off
Sorption is inactive (default)
Definition: gwt-mst.f90:38
subroutine mst_bd(this, isuppress_output, model_budget)
@ brief Calculate budget terms for package
Definition: gwt-mst.f90:881
subroutine allocate_arrays(this, nodes)
@ brief Allocate arrays for package
Definition: gwt-mst.f90:1079
real(dp) function get_freundlich_derivative(conc, kf, a)
@ brief Calculate sorption derivative using Freundlich
Definition: gwt-mst.f90:1543
real(dp) function get_langmuir_kd(conc, kl, sbar)
@ brief Get effective Langmuir distribution coefficient
Definition: gwt-mst.f90:1596
subroutine mst_cq_dcy_srb(this, nodes, cnew, cold, flowja)
@ brief Calculate decay-sorption terms for package
Definition: gwt-mst.f90:751
real(dp) function get_freundlich_conc(conc, kf, a)
@ brief Calculate sorption concentration using Freundlich
Definition: gwt-mst.f90:1505
real(dp) function get_volfracm(this, node)
@ brief Return mobile domain volume fraction
Definition: gwt-mst.f90:1491
subroutine mst_ot_dv(this, idvsave)
Save sorbate concentration array to binary file.
Definition: gwt-mst.f90:977
subroutine mst_ot_flow(this, icbcfl, icbcun)
@ brief Output flow terms for package
Definition: gwt-mst.f90:924
subroutine read_options(this)
@ brief Read options for package
Definition: gwt-mst.f90:1168
subroutine mst_ar(this, dis, ibound)
@ brief Allocate and read method for package
Definition: gwt-mst.f90:144
subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
@ brief Fill sorption coefficient method for package
Definition: gwt-mst.f90:321
subroutine mst_cq_dcy(this, nodes, cnew, cold, flowja)
@ brief Calculate decay terms for package
Definition: gwt-mst.f90:639
subroutine read_data(this)
@ brief Read data for package
Definition: gwt-mst.f90:1270
subroutine allocate_scalars(this)
@ brief Allocate scalar variables for package
Definition: gwt-mst.f90:1055
subroutine mst_fc_dcy_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, cnew, kiter)
@ brief Fill sorption-decay coefficient method for package
Definition: gwt-mst.f90:454
subroutine mst_da(this)
@ brief Deallocate
Definition: gwt-mst.f90:1015
subroutine mst_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, idxglo, rhs, kiter)
@ brief Fill decay coefficient method for package
Definition: gwt-mst.f90:263
integer(i4b), parameter nbditems
Definition: gwt-mst.f90:27
subroutine mst_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, rhs, kiter)
@ brief Fill coefficient method for package
Definition: gwt-mst.f90:178
real(dp) function get_freundlich_kd(conc, kf, a)
@ brief Get effective Freundlich distribution coefficient
Definition: gwt-mst.f90:1579
subroutine mst_fc_sto(this, nodes, cold, nja, matrix_sln, idxglo, rhs)
@ brief Fill storage coefficient method for package
Definition: gwt-mst.f90:216
real(dp) function get_langmuir_derivative(conc, kl, sbar)
@ brief Calculate sorption derivative using Langmuir
Definition: gwt-mst.f90:1562
subroutine, public mst_cr(mstobj, name_model, inunit, iout, fmi)
@ brief Create a new package object
Definition: gwt-mst.f90:114
subroutine mst_cq_sto(this, nodes, cnew, cold, flowja)
@ brief Calculate storage terms for package
Definition: gwt-mst.f90:591
real(dp) function get_langmuir_conc(conc, kl, sbar)
@ brief Calculate sorption concentration using Langmuir
Definition: gwt-mst.f90:1524
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
This module defines variable data types.
Definition: kind.f90:8
This module contains the base numerical package type.
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_warning(msg, substring)
Store warning message.
Definition: Sim.f90:236
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
character(len=maxcharlen) warnmsg
warning message string
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
Derived type for the Budget object.
Definition: Budget.f90:39
@ brief Mobile storage and transfer
Definition: gwt-mst.f90:50