28 character(len=LENBUDTXT),
dimension(NBDITEMS) ::
budtxt
29 data budtxt/
' STORAGE-AQUEOUS',
' DECAY-AQUEOUS', &
30 ' STORAGE-SORBED',
' DECAY-SORBED'/
53 real(dp),
dimension(:),
pointer,
contiguous :: porosity => null()
54 real(dp),
dimension(:),
pointer,
contiguous :: thetam => null()
55 real(dp),
dimension(:),
pointer,
contiguous :: volfracim => null()
56 real(dp),
dimension(:),
pointer,
contiguous :: ratesto => null()
59 integer(I4B),
pointer :: idcy => null()
60 real(dp),
dimension(:),
pointer,
contiguous :: decay => null()
61 real(dp),
dimension(:),
pointer,
contiguous :: decay_sorbed => null()
62 real(dp),
dimension(:),
pointer,
contiguous :: ratedcy => null()
63 real(dp),
dimension(:),
pointer,
contiguous :: decaylast => null()
64 real(dp),
dimension(:),
pointer,
contiguous :: decayslast => null()
67 integer(I4B),
pointer :: isrb => null()
68 integer(I4B),
pointer :: ioutsorbate => null()
69 real(dp),
dimension(:),
pointer,
contiguous :: bulk_density => null()
70 real(dp),
dimension(:),
pointer,
contiguous :: distcoef => null()
71 real(dp),
dimension(:),
pointer,
contiguous :: sp2 => null()
72 real(dp),
dimension(:),
pointer,
contiguous :: ratesrb => null()
73 real(dp),
dimension(:),
pointer,
contiguous :: ratedcys => null()
74 real(dp),
dimension(:),
pointer,
contiguous :: csrb => null()
77 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
113 subroutine mst_cr(mstobj, name_model, inunit, iout, fmi)
116 character(len=*),
intent(in) :: name_model
117 integer(I4B),
intent(in) :: inunit
118 integer(I4B),
intent(in) :: iout
125 call mstobj%set_names(1, name_model,
'MST',
'MST')
128 call mstobj%allocate_scalars()
131 mstobj%inunit = inunit
136 call mstobj%parser%Initialize(mstobj%inunit, mstobj%iout)
148 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
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, //)"
156 write (this%iout, fmtmst) this%inunit
159 call this%read_options()
163 this%ibound => ibound
166 call this%allocate_arrays(dis%nodes)
169 call this%read_data()
176 subroutine mst_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, &
181 integer,
intent(in) :: nodes
182 real(DP),
intent(in),
dimension(nodes) :: cold
183 integer(I4B),
intent(in) :: nja
185 integer(I4B),
intent(in),
dimension(nja) :: idxglo
186 real(DP),
intent(inout),
dimension(nodes) :: rhs
187 real(DP),
intent(in),
dimension(nodes) :: cnew
188 integer(I4B),
intent(in) :: kiter
191 call this%mst_fc_sto(nodes, cold, nja, matrix_sln, idxglo, rhs)
195 call this%mst_fc_dcy(nodes, cold, cnew, nja, matrix_sln, idxglo, &
201 call this%mst_fc_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
206 call this%mst_fc_dcy_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, &
215 subroutine mst_fc_sto(this, nodes, cold, nja, matrix_sln, idxglo, rhs)
220 integer,
intent(in) :: nodes
221 real(DP),
intent(in),
dimension(nodes) :: cold
222 integer(I4B),
intent(in) :: nja
224 integer(I4B),
intent(in),
dimension(nja) :: idxglo
225 real(DP),
intent(inout),
dimension(nodes) :: rhs
227 integer(I4B) :: n, idiag
229 real(DP) :: hhcof, rrhs
230 real(DP) :: vnew, vold
236 do n = 1, this%dis%nodes
239 if (this%ibound(n) <= 0) cycle
242 vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
243 this%fmi%gwfsat(n) * this%thetam(n)
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
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
261 subroutine mst_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, &
267 integer,
intent(in) :: nodes
268 real(DP),
intent(in),
dimension(nodes) :: cold
269 real(DP),
intent(in),
dimension(nodes) :: cnew
270 integer(I4B),
intent(in) :: nja
272 integer(I4B),
intent(in),
dimension(nja) :: idxglo
273 real(DP),
intent(inout),
dimension(nodes) :: rhs
274 integer(I4B),
intent(in) :: kiter
276 integer(I4B) :: n, idiag
277 real(DP) :: hhcof, rrhs
280 real(DP) :: decay_rate
283 do n = 1, this%dis%nodes
286 if (this%ibound(n) <= 0) cycle
289 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
290 swtpdt = this%fmi%gwfsat(n)
293 idiag = this%dis%con%ia(n)
294 select case (this%idcy)
299 hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
300 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
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
319 subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, &
325 integer,
intent(in) :: nodes
326 real(DP),
intent(in),
dimension(nodes) :: cold
327 integer(I4B),
intent(in) :: nja
329 integer(I4B),
intent(in),
dimension(nja) :: idxglo
330 real(DP),
intent(inout),
dimension(nodes) :: rhs
331 real(DP),
intent(in),
dimension(nodes) :: cnew
333 integer(I4B) :: n, idiag
335 real(DP) :: hhcof, rrhs
336 real(DP) :: swt, swtpdt
347 do n = 1, this%dis%nodes
350 if (this%ibound(n) <= 0) cycle
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)
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)
369 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
370 rhs(n) = rhs(n) + rrhs
379 subroutine mst_srb_term(isrb, volfracm, rhobm, vcell, tled, cnew, cold, &
380 swnew, swold, const1, const2, rate, hcofval, rhsval)
382 integer(I4B),
intent(in) :: isrb
383 real(DP),
intent(in) :: volfracm
384 real(DP),
intent(in) :: rhobm
385 real(DP),
intent(in) :: vcell
386 real(DP),
intent(in) :: tled
387 real(DP),
intent(in) :: cnew
388 real(DP),
intent(in) :: cold
389 real(DP),
intent(in) :: swnew
390 real(DP),
intent(in) :: swold
391 real(DP),
intent(in) :: const1
392 real(DP),
intent(in) :: const2
393 real(DP),
intent(out),
optional :: rate
394 real(DP),
intent(out),
optional :: hcofval
395 real(DP),
intent(out),
optional :: rhsval
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
415 cavg =
dhalf * (cold + cnew)
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
438 if (
present(rhsval))
then
439 rhsval = term * derv * swavg * cold - term * cbaravg * (swnew - swold)
441 if (
present(rate))
then
442 rate = term * derv * swavg * (cnew - cold) &
443 + term * cbaravg * (swnew - swold)
458 integer,
intent(in) :: nodes
459 real(DP),
intent(in),
dimension(nodes) :: cold
460 integer(I4B),
intent(in) :: nja
462 integer(I4B),
intent(in),
dimension(nja) :: idxglo
463 real(DP),
intent(inout),
dimension(nodes) :: rhs
464 real(DP),
intent(in),
dimension(nodes) :: cnew
465 integer(I4B),
intent(in) :: kiter
467 integer(I4B) :: n, idiag
468 real(DP) :: hhcof, rrhs
476 real(DP) :: decay_rate
481 do n = 1, this%dis%nodes
484 if (this%ibound(n) <= 0) cycle
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
498 select case (this%idcy)
501 select case (this%isrb)
506 hhcof = -term * distcoef
522 if (distcoef >
dzero)
then
523 select case (this%isrb)
525 csrbold = cold(n) * distcoef
526 csrbnew = cnew(n) * distcoef
536 this%decayslast(n), &
537 kiter, csrbold, csrbnew,
delt)
538 this%decayslast(n) = decay_rate
539 rrhs = decay_rate * volfracm * rhobm * swnew * vcell
544 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
545 rhs(n) = rhs(n) + rrhs
554 subroutine mst_cq(this, nodes, cnew, cold, flowja)
557 integer(I4B),
intent(in) :: nodes
558 real(DP),
intent(in),
dimension(nodes) :: cnew
559 real(DP),
intent(in),
dimension(nodes) :: cold
560 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
563 call this%mst_cq_sto(nodes, cnew, cold, flowja)
567 call this%mst_cq_dcy(nodes, cnew, cold, flowja)
572 call this%mst_cq_srb(nodes, cnew, cold, flowja)
577 call this%mst_cq_dcy_srb(nodes, cnew, cold, flowja)
582 call this%mst_calc_csrb(cnew)
595 integer(I4B),
intent(in) :: nodes
596 real(DP),
intent(in),
dimension(nodes) :: cnew
597 real(DP),
intent(in),
dimension(nodes) :: cold
598 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
601 integer(I4B) :: idiag
604 real(DP) :: vnew, vold
605 real(DP) :: hhcof, rrhs
612 this%ratesto(n) =
dzero
615 if (this%ibound(n) <= 0) cycle
618 vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
619 this%fmi%gwfsat(n) * this%thetam(n)
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
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
643 integer(I4B),
intent(in) :: nodes
644 real(DP),
intent(in),
dimension(nodes) :: cnew
645 real(DP),
intent(in),
dimension(nodes) :: cold
646 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
649 integer(I4B) :: idiag
652 real(DP) :: hhcof, rrhs
654 real(DP) :: decay_rate
662 this%ratedcy(n) =
dzero
663 if (this%ibound(n) <= 0) cycle
666 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
667 swtpdt = this%fmi%gwfsat(n)
674 hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
677 0, cold(n), cnew(n),
delt)
678 rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
680 rate = hhcof * cnew(n) - rrhs
681 this%ratedcy(n) = rate
682 idiag = this%dis%con%ia(n)
683 flowja(idiag) = flowja(idiag) + rate
697 integer(I4B),
intent(in) :: nodes
698 real(DP),
intent(in),
dimension(nodes) :: cnew
699 real(DP),
intent(in),
dimension(nodes) :: cold
700 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
703 integer(I4B) :: idiag
706 real(DP) :: swt, swtpdt
720 this%ratesrb(n) =
dzero
723 if (this%ibound(n) <= 0) cycle
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)
736 call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
737 cold(n), swtpdt, swt, const1, const2, &
739 this%ratesrb(n) = rate
740 idiag = this%dis%con%ia(n)
741 flowja(idiag) = flowja(idiag) + rate
755 integer(I4B),
intent(in) :: nodes
756 real(DP),
intent(in),
dimension(nodes) :: cnew
757 real(DP),
intent(in),
dimension(nodes) :: cold
758 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
761 integer(I4B) :: idiag
763 real(DP) :: hhcof, rrhs
773 real(DP) :: decay_rate
780 this%ratedcys(n) =
dzero
783 if (this%ibound(n) <= 0) cycle
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
796 select case (this%idcy)
799 select case (this%isrb)
804 hhcof = -term * distcoef
820 if (distcoef >
dzero)
then
821 select case (this%isrb)
823 csrbold = cold(n) * distcoef
824 csrbnew = cnew(n) * distcoef
833 this%decayslast(n), &
834 0, csrbold, csrbnew,
delt)
835 rrhs = decay_rate * volfracm * rhobm * swnew * vcell
840 rate = hhcof * cnew(n) - rrhs
841 this%ratedcys(n) = rate
842 idiag = this%dis%con%ia(n)
843 flowja(idiag) = flowja(idiag) + rate
853 real(DP),
intent(in),
dimension(:) :: cnew
862 if (this%ibound(n) > 0)
then
863 distcoef = this%distcoef(n)
864 select case (this%isrb)
866 csrb = cnew(n) * distcoef
880 subroutine mst_bd(this, isuppress_output, model_budget)
886 integer(I4B),
intent(in) :: isuppress_output
887 type(
budgettype),
intent(inout) :: model_budget
894 call model_budget%addentry(rin, rout,
delt,
budtxt(1), &
895 isuppress_output, rowlabel=this%packName)
900 call model_budget%addentry(rin, rout,
delt,
budtxt(2), &
901 isuppress_output, rowlabel=this%packName)
907 call model_budget%addentry(rin, rout,
delt,
budtxt(3), &
908 isuppress_output, rowlabel=this%packName)
914 call model_budget%addentry(rin, rout,
delt,
budtxt(4), &
915 isuppress_output, rowlabel=this%packName)
926 integer(I4B),
intent(in) :: icbcfl
927 integer(I4B),
intent(in) :: icbcun
929 integer(I4B) :: ibinun
930 integer(I4B) :: iprint, nvaluesp, nwidthp
931 character(len=1) :: cdatafmp =
' ', editdesc =
' '
935 if (this%ipakcb < 0)
then
937 elseif (this%ipakcb == 0)
then
942 if (icbcfl == 0) ibinun = 0
945 if (ibinun /= 0)
then
950 call this%dis%record_array(this%ratesto, this%iout, iprint, -ibinun, &
951 budtxt(1), cdatafmp, nvaluesp, &
952 nwidthp, editdesc, dinact)
956 call this%dis%record_array(this%ratedcy, this%iout, iprint, -ibinun, &
957 budtxt(2), cdatafmp, nvaluesp, &
958 nwidthp, editdesc, dinact)
962 call this%dis%record_array(this%ratesrb, this%iout, iprint, -ibinun, &
963 budtxt(3), cdatafmp, nvaluesp, &
964 nwidthp, editdesc, dinact)
968 call this%dis%record_array(this%ratedcys, this%iout, iprint, -ibinun, &
969 budtxt(4), cdatafmp, nvaluesp, &
970 nwidthp, editdesc, dinact)
979 integer(I4B),
intent(in) :: idvsave
981 character(len=1) :: cdatafmp =
' ', editdesc =
' '
982 integer(I4B) :: ibinun
983 integer(I4B) :: iprint
984 integer(I4B) :: nvaluesp
985 integer(I4B) :: nwidthp
989 if (this%ioutsorbate /= 0)
then
994 if (idvsave == 0) ibinun = 0
997 if (ibinun /= 0)
then
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)
1021 if (this%inunit > 0)
then
1040 this%ibound => null()
1047 call this%NumericalPackageType%da()
1061 call this%NumericalPackageType%allocate_scalars()
1065 call mem_allocate(this%ioutsorbate,
'IOUTSORBATE', this%memoryPath)
1070 this%ioutsorbate = 0
1084 integer(I4B),
intent(in) :: nodes
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)
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)
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)
1106 call mem_allocate(this%ratedcys, this%dis%nodes,
'RATEDCYS', &
1108 call mem_allocate(this%decayslast, this%dis%nodes,
'DECAYSLAST', &
1111 call mem_allocate(this%ratedcys, 1,
'RATEDCYS', this%memoryPath)
1112 call mem_allocate(this%decayslast, 1,
'DECAYSLAST', this%memoryPath)
1114 call mem_allocate(this%decay_sorbed, 1,
'DECAY_SORBED', &
1119 call mem_allocate(this%bulk_density, 1,
'BULK_DENSITY', 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)
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)
1132 call mem_allocate(this%sp2, nodes,
'SP2', this%memoryPath)
1138 this%porosity(n) =
dzero
1139 this%thetam(n) =
dzero
1140 this%volfracim(n) =
dzero
1141 this%ratesto(n) =
dzero
1143 do n = 1,
size(this%decay)
1144 this%decay(n) =
dzero
1145 this%ratedcy(n) =
dzero
1146 this%decaylast(n) =
dzero
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
1154 do n = 1,
size(this%sp2)
1157 do n = 1,
size(this%ratedcys)
1158 this%ratedcys(n) =
dzero
1159 this%decayslast(n) =
dzero
1174 character(len=LINELENGTH) :: keyword
1175 character(len=LINELENGTH) :: sorption
1176 character(len=MAXCHARLEN) :: fname
1177 integer(I4B) :: ierr
1178 logical :: isfound, endOfBlock
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)"
1198 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, blockrequired=.false., &
1199 supportopenclose=.true.)
1203 write (this%iout,
'(1x,a)')
'PROCESSING MOBILE STORAGE AND TRANSFER OPTIONS'
1205 call this%parser%GetNextLine(endofblock)
1206 if (endofblock)
exit
1207 call this%parser%GetStringCaps(keyword)
1208 select case (keyword)
1211 write (this%iout, fmtisvflow)
1213 call store_error(
'SORBTION is not a valid option. Use &
1214 &SORPTION instead.')
1215 call this%parser%StoreErrorUnit()
1217 call this%parser%GetStringCaps(sorption)
1218 select case (sorption)
1221 write (this%iout, fmtlinear)
1224 write (this%iout, fmtfreundlich)
1227 write (this%iout, fmtlangmuir)
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()
1235 case (
'FIRST_ORDER_DECAY')
1237 write (this%iout, fmtidcy1)
1238 case (
'ZERO_ORDER_DECAY')
1240 write (this%iout, fmtidcy2)
1242 call this%parser%GetStringCaps(keyword)
1243 if (keyword ==
'FILEOUT')
then
1244 call this%parser%GetString(fname)
1246 call openfile(this%ioutsorbate, this%iout, fname,
'DATA(BINARY)', &
1248 write (this%iout, fmtfileout) &
1249 'SORBATE', fname, this%ioutsorbate
1251 errmsg =
'Optional SORBATE keyword must be '// &
1252 'followed by FILEOUT.'
1256 write (
errmsg,
'(a,a)')
'Unknown MST option: ', trim(keyword)
1258 call this%parser%StoreErrorUnit()
1261 write (this%iout,
'(1x,a)')
'END OF MOBILE STORAGE AND TRANSFER OPTIONS'
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
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'/
1295 call this%parser%GetBlock(
'GRIDDATA', isfound, ierr)
1297 write (this%iout,
'(1x,a)')
'PROCESSING GRIDDATA'
1299 call this%parser%GetNextLine(endofblock)
1300 if (endofblock)
exit
1301 call this%parser%GetStringCaps(keyword)
1302 call this%parser%GetRemainingLine(line)
1304 select case (keyword)
1306 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1307 this%parser%iuactive, this%porosity, &
1310 case (
'BULK_DENSITY')
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))
1321 trim(this%memoryPath))
1322 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1323 this%parser%iuactive, this%distcoef, &
1329 trim(this%memoryPath))
1330 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1331 this%parser%iuactive, this%decay, &
1334 case (
'DECAY_SORBED')
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))
1344 trim(this%memoryPath))
1345 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1346 this%parser%iuactive, this%sp2, &
1350 write (
errmsg,
'(a,a)')
'Unknown GRIDDATA tag: ', trim(keyword)
1352 call this%parser%StoreErrorUnit()
1355 write (this%iout,
'(1x,a)')
'END PROCESSING GRIDDATA'
1357 write (
errmsg,
'(a)')
'Required GRIDDATA block not found.'
1359 call this%parser%StoreErrorUnit()
1363 if (.not. lname(1))
then
1364 write (
errmsg,
'(a)')
'POROSITY not specified in GRIDDATA block.'
1371 if (.not. lname(2))
then
1372 write (
errmsg,
'(a)')
'Sorption is active but BULK_DENSITY &
1373 ¬ specified. BULK_DENSITY must be specified in GRIDDATA block.'
1376 if (.not. lname(3))
then
1377 write (
errmsg,
'(a)')
'Sorption is active but distribution &
1378 &coefficient not specified. DISTCOEF must be specified in &
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 &
1392 write (
warnmsg,
'(a)')
'Sorption is not active but &
1393 &BULK_DENSITY was specified. BULK_DENSITY will have no affect on &
1394 &simulation results.'
1396 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1399 write (
warnmsg,
'(a)')
'Sorption is not active but &
1400 &distribution coefficient was specified. DISTCOEF will have &
1401 &no affect on simulation results.'
1403 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1406 write (
warnmsg,
'(a)')
'Sorption is not active but &
1407 &SP2 was specified. SP2 will have &
1408 &no affect on simulation results.'
1410 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
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.'
1422 if (.not. lname(5))
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.'
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.'
1440 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
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.'
1447 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1453 call this%parser%StoreErrorUnit()
1457 do n = 1,
size(this%porosity)
1458 this%thetam(n) = this%porosity(n)
1470 real(DP),
dimension(:),
intent(in) :: volfracim
1475 do n = 1, this%dis%nodes
1476 this%volfracim(n) = this%volfracim(n) + volfracim(n)
1481 do n = 1, this%dis%nodes
1482 this%thetam(n) = this%get_volfracm(n) * this%porosity(n)
1493 integer(I4B),
intent(in) :: node
1495 real(dp) :: volfracm
1497 volfracm =
done - this%volfracim(node)
1506 real(dp),
intent(in) :: conc
1507 real(dp),
intent(in) :: kf
1508 real(dp),
intent(in) :: a
1512 if (conc >
dzero)
then
1525 real(dp),
intent(in) :: conc
1526 real(dp),
intent(in) :: kl
1527 real(dp),
intent(in) :: sbar
1531 if (conc >
dzero)
then
1532 cbar = (kl * sbar * conc) / (
done + kl * conc)
1544 real(dp),
intent(in) :: conc
1545 real(dp),
intent(in) :: kf
1546 real(dp),
intent(in) :: a
1550 if (conc >
dzero)
then
1551 derv = kf * a * conc**(a -
done)
1563 real(dp),
intent(in) :: conc
1564 real(dp),
intent(in) :: kl
1565 real(dp),
intent(in) :: sbar
1569 if (conc >
dzero)
then
1570 derv = (kl * sbar) / (
done + kl * conc)**
dtwo
1580 real(dp),
intent(in) :: conc
1581 real(dp),
intent(in) :: kf
1582 real(dp),
intent(in) :: a
1586 if (conc >
dzero)
then
1587 kd = kf * conc**(a -
done)
1597 real(dp),
intent(in) :: conc
1598 real(dp),
intent(in) :: kl
1599 real(dp),
intent(in) :: sbar
1603 if (conc >
dzero)
then
1604 kd = (kl * sbar) / (
done + kl * conc)
1619 cold, cnew, delt)
result(decay_rate)
1621 real(dp),
intent(in) :: decay_rate_usr
1622 real(dp),
intent(in) :: decay_rate_last
1623 integer(I4B),
intent(in) :: kiter
1624 real(dp),
intent(in) :: cold
1625 real(dp),
intent(in) :: cnew
1626 real(dp),
intent(in) :: delt
1628 real(dp) :: decay_rate
1631 if (decay_rate_usr <
dzero)
then
1634 decay_rate = decay_rate_usr
1640 if (kiter == 1)
then
1641 decay_rate = min(decay_rate_usr, cold / delt)
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
1649 decay_rate = min(decay_rate_usr, decay_rate)
1651 decay_rate = max(decay_rate,
dzero)
This module contains the BudgetModule.
subroutine, public rate_accumulator(flow, rin, rout)
@ brief Rate accumulator subroutine
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
@ mnormal
normal output mode
real(dp), parameter dhnoflo
real no flow constant
real(dp), parameter dhalf
real constant 1/2
integer(i4b), parameter izero
integer constant zero
real(dp), parameter dzero
real constant zero
integer(i4b), parameter maxcharlen
maximum length of char string
real(dp), parameter dtwo
real constant 2
integer(i4b), parameter lenbudtxt
maximum length of a budget component names
real(dp), parameter done
real constant 1
– @ brief Mobile Storage and Transfer (MST) Module
subroutine mst_srb_term(isrb, volfracm, rhobm, vcell, tled, cnew, cold, swnew, swold, const1, const2, rate, hcofval, rhsval)
@ brief Calculate sorption terms
subroutine mst_cq_srb(this, nodes, cnew, cold, flowja)
@ brief Calculate sorption terms for package
subroutine mst_calc_csrb(this, cnew)
@ brief Calculate sorbed concentration
subroutine addto_volfracim(this, volfracim)
@ brief Add volfrac values to volfracim
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
character(len=lenbudtxt), dimension(nbditems) budtxt
subroutine mst_cq(this, nodes, cnew, cold, flowja)
@ brief Calculate flows for package
@ decay_zero_order
Zeroth-order decay.
@ decay_first_order
First-order decay.
@ sorption_freund
Freundlich sorption between aqueous and solid phases.
@ sorption_lang
Langmuir sorption between aqueous and solid phases.
@ sorption_linear
Linear sorption between aqueous and solid phases.
@ decay_off
Decay (or production) of mass inactive (default)
@ sorption_off
Sorption is inactive (default)
subroutine mst_bd(this, isuppress_output, model_budget)
@ brief Calculate budget terms for package
subroutine allocate_arrays(this, nodes)
@ brief Allocate arrays for package
real(dp) function get_freundlich_derivative(conc, kf, a)
@ brief Calculate sorption derivative using Freundlich
real(dp) function get_langmuir_kd(conc, kl, sbar)
@ brief Get effective Langmuir distribution coefficient
subroutine mst_cq_dcy_srb(this, nodes, cnew, cold, flowja)
@ brief Calculate decay-sorption terms for package
real(dp) function get_freundlich_conc(conc, kf, a)
@ brief Calculate sorption concentration using Freundlich
real(dp) function get_volfracm(this, node)
@ brief Return mobile domain volume fraction
subroutine mst_ot_dv(this, idvsave)
Save sorbate concentration array to binary file.
subroutine mst_ot_flow(this, icbcfl, icbcun)
@ brief Output flow terms for package
subroutine read_options(this)
@ brief Read options for package
subroutine mst_ar(this, dis, ibound)
@ brief Allocate and read method for package
subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
@ brief Fill sorption coefficient method for package
subroutine mst_cq_dcy(this, nodes, cnew, cold, flowja)
@ brief Calculate decay terms for package
subroutine read_data(this)
@ brief Read data for package
subroutine allocate_scalars(this)
@ brief Allocate scalar variables for package
subroutine mst_fc_dcy_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, cnew, kiter)
@ brief Fill sorption-decay coefficient method for package
subroutine mst_da(this)
@ brief Deallocate
subroutine mst_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, idxglo, rhs, kiter)
@ brief Fill decay coefficient method for package
integer(i4b), parameter nbditems
subroutine mst_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, rhs, kiter)
@ brief Fill coefficient method for package
real(dp) function get_freundlich_kd(conc, kf, a)
@ brief Get effective Freundlich distribution coefficient
subroutine mst_fc_sto(this, nodes, cold, nja, matrix_sln, idxglo, rhs)
@ brief Fill storage coefficient method for package
real(dp) function get_langmuir_derivative(conc, kl, sbar)
@ brief Calculate sorption derivative using Langmuir
subroutine, public mst_cr(mstobj, name_model, inunit, iout, fmi)
@ brief Create a new package object
subroutine mst_cq_sto(this, nodes, cnew, cold, flowja)
@ brief Calculate storage terms for package
real(dp) function get_langmuir_conc(conc, kl, sbar)
@ brief Calculate sorption concentration using Langmuir
This module defines variable data types.
This module contains the base numerical package type.
This module contains simulation methods.
subroutine, public store_warning(msg, substring)
Store warning message.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
This module contains simulation variables.
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
Derived type for the Budget object.
@ brief Mobile storage and transfer