27 character(len=LENBUDTXT),
dimension(NBDITEMS) ::
budtxt
28 data budtxt/
' STORAGE-AQUEOUS',
' DECAY-AQUEOUS', &
29 ' STORAGE-SORBED',
' DECAY-SORBED'/
40 real(dp),
dimension(:),
pointer,
contiguous :: porosity => null()
41 real(dp),
dimension(:),
pointer,
contiguous :: thetam => null()
42 real(dp),
dimension(:),
pointer,
contiguous :: volfracim => null()
43 real(dp),
dimension(:),
pointer,
contiguous :: ratesto => null()
46 integer(I4B),
pointer :: idcy => null()
47 real(dp),
dimension(:),
pointer,
contiguous :: decay => null()
48 real(dp),
dimension(:),
pointer,
contiguous :: decay_sorbed => null()
49 real(dp),
dimension(:),
pointer,
contiguous :: ratedcy => null()
50 real(dp),
dimension(:),
pointer,
contiguous :: decaylast => null()
51 real(dp),
dimension(:),
pointer,
contiguous :: decayslast => null()
54 integer(I4B),
pointer :: isrb => null()
55 real(dp),
dimension(:),
pointer,
contiguous :: bulk_density => null()
56 real(dp),
dimension(:),
pointer,
contiguous :: distcoef => null()
57 real(dp),
dimension(:),
pointer,
contiguous :: sp2 => null()
58 real(dp),
dimension(:),
pointer,
contiguous :: ratesrb => null()
59 real(dp),
dimension(:),
pointer,
contiguous :: ratedcys => null()
62 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
97 subroutine mst_cr(mstobj, name_model, inunit, iout, fmi)
100 character(len=*),
intent(in) :: name_model
101 integer(I4B),
intent(in) :: inunit
102 integer(I4B),
intent(in) :: iout
109 call mstobj%set_names(1, name_model,
'MST',
'MST')
112 call mstobj%allocate_scalars()
115 mstobj%inunit = inunit
120 call mstobj%parser%Initialize(mstobj%inunit, mstobj%iout)
136 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
139 character(len=*),
parameter :: fmtmst = &
140 "(1x,/1x,'MST -- MOBILE STORAGE AND TRANSFER PACKAGE, VERSION 1, &
141 &7/29/2020 INPUT READ FROM UNIT ', i0, //)"
144 write (this%iout, fmtmst) this%inunit
147 call this%read_options()
151 this%ibound => ibound
154 call this%allocate_arrays(dis%nodes)
157 call this%read_data()
168 subroutine mst_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, &
173 integer,
intent(in) :: nodes
174 real(DP),
intent(in),
dimension(nodes) :: cold
175 integer(I4B),
intent(in) :: nja
177 integer(I4B),
intent(in),
dimension(nja) :: idxglo
178 real(DP),
intent(inout),
dimension(nodes) :: rhs
179 real(DP),
intent(in),
dimension(nodes) :: cnew
180 integer(I4B),
intent(in) :: kiter
184 call this%mst_fc_sto(nodes, cold, nja, matrix_sln, idxglo, rhs)
187 if (this%idcy /= 0)
then
188 call this%mst_fc_dcy(nodes, cold, cnew, nja, matrix_sln, idxglo, &
193 if (this%isrb /= 0)
then
194 call this%mst_fc_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, cnew)
198 if (this%isrb /= 0 .and. this%idcy /= 0)
then
199 call this%mst_fc_dcy_srb(nodes, cold, nja, matrix_sln, idxglo, rhs, &
212 subroutine mst_fc_sto(this, nodes, cold, nja, matrix_sln, idxglo, rhs)
217 integer,
intent(in) :: nodes
218 real(DP),
intent(in),
dimension(nodes) :: cold
219 integer(I4B),
intent(in) :: nja
221 integer(I4B),
intent(in),
dimension(nja) :: idxglo
222 real(DP),
intent(inout),
dimension(nodes) :: rhs
224 integer(I4B) :: n, idiag
226 real(DP) :: hhcof, rrhs
227 real(DP) :: vnew, vold
233 do n = 1, this%dis%nodes
236 if (this%ibound(n) <= 0) cycle
239 vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
240 this%fmi%gwfsat(n) * this%thetam(n)
242 if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) *
delt
243 if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) *
delt
247 rrhs = -vold * tled * cold(n)
248 idiag = this%dis%con%ia(n)
249 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
250 rhs(n) = rhs(n) + rrhs
262 subroutine mst_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, &
268 integer,
intent(in) :: nodes
269 real(DP),
intent(in),
dimension(nodes) :: cold
270 real(DP),
intent(in),
dimension(nodes) :: cnew
271 integer(I4B),
intent(in) :: nja
273 integer(I4B),
intent(in),
dimension(nja) :: idxglo
274 real(DP),
intent(inout),
dimension(nodes) :: rhs
275 integer(I4B),
intent(in) :: kiter
277 integer(I4B) :: n, idiag
278 real(DP) :: hhcof, rrhs
281 real(DP) :: decay_rate
284 do n = 1, this%dis%nodes
287 if (this%ibound(n) <= 0) cycle
290 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
291 swtpdt = this%fmi%gwfsat(n)
294 idiag = this%dis%con%ia(n)
295 if (this%idcy == 1)
then
299 hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
300 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
301 elseif (this%idcy == 2)
then
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
323 subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, &
329 integer,
intent(in) :: nodes
330 real(DP),
intent(in),
dimension(nodes) :: cold
331 integer(I4B),
intent(in) :: nja
333 integer(I4B),
intent(in),
dimension(nja) :: idxglo
334 real(DP),
intent(inout),
dimension(nodes) :: rhs
335 real(DP),
intent(in),
dimension(nodes) :: cnew
337 integer(I4B) :: n, idiag
339 real(DP) :: hhcof, rrhs
340 real(DP) :: swt, swtpdt
351 do n = 1, this%dis%nodes
354 if (this%ibound(n) <= 0) cycle
357 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
358 swtpdt = this%fmi%gwfsat(n)
359 swt = this%fmi%gwfsatold(n,
delt)
360 idiag = this%dis%con%ia(n)
361 const1 = this%distcoef(n)
363 if (this%isrb > 1) const2 = this%sp2(n)
364 volfracm = this%get_volfracm(n)
365 rhobm = this%bulk_density(n)
366 call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
367 cold(n), swtpdt, swt, const1, const2, &
368 hcofval=hhcof, rhsval=rrhs)
371 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
372 rhs(n) = rhs(n) + rrhs
385 subroutine mst_srb_term(isrb, volfracm, rhobm, vcell, tled, cnew, cold, &
386 swnew, swold, const1, const2, rate, hcofval, rhsval)
389 integer(I4B),
intent(in) :: isrb
390 real(DP),
intent(in) :: volfracm
391 real(DP),
intent(in) :: rhobm
392 real(DP),
intent(in) :: vcell
393 real(DP),
intent(in) :: tled
394 real(DP),
intent(in) :: cnew
395 real(DP),
intent(in) :: cold
396 real(DP),
intent(in) :: swnew
397 real(DP),
intent(in) :: swold
398 real(DP),
intent(in) :: const1
399 real(DP),
intent(in) :: const2
400 real(DP),
intent(out),
optional :: rate
401 real(DP),
intent(out),
optional :: hcofval
402 real(DP),
intent(out),
optional :: rhsval
415 term = -volfracm * rhobm * vcell * tled * const1
416 if (
present(hcofval)) hcofval = term * swnew
417 if (
present(rhsval)) rhsval = term * swold * cold
418 if (
present(rate)) rate = term * swnew * cnew - term * swold * cold
422 cavg =
dhalf * (cold + cnew)
430 else if (isrb == 3)
then
438 term = -volfracm * rhobm * vcell * tled
439 cbaravg = (cbarold + cbarnew) *
dhalf
440 swavg = (swnew + swold) *
dhalf
441 if (
present(hcofval))
then
442 hcofval = term * derv * swavg
444 if (
present(rhsval))
then
445 rhsval = term * derv * swavg * cold - term * cbaravg * (swnew - swold)
447 if (
present(rate))
then
448 rate = term * derv * swavg * (cnew - cold) &
449 + term * cbaravg * (swnew - swold)
466 integer,
intent(in) :: nodes
467 real(DP),
intent(in),
dimension(nodes) :: cold
468 integer(I4B),
intent(in) :: nja
470 integer(I4B),
intent(in),
dimension(nja) :: idxglo
471 real(DP),
intent(inout),
dimension(nodes) :: rhs
472 real(DP),
intent(in),
dimension(nodes) :: cnew
473 integer(I4B),
intent(in) :: kiter
475 integer(I4B) :: n, idiag
476 real(DP) :: hhcof, rrhs
484 real(DP) :: decay_rate
489 do n = 1, this%dis%nodes
492 if (this%ibound(n) <= 0) cycle
497 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
498 swnew = this%fmi%gwfsat(n)
499 distcoef = this%distcoef(n)
500 idiag = this%dis%con%ia(n)
501 volfracm = this%get_volfracm(n)
502 rhobm = this%bulk_density(n)
503 term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
506 if (this%idcy == 1)
then
508 if (this%isrb == 1)
then
512 hhcof = -term * distcoef
513 else if (this%isrb == 2)
then
518 else if (this%isrb == 3)
then
524 elseif (this%idcy == 2)
then
528 if (distcoef >
dzero)
then
530 if (this%isrb == 1)
then
531 csrbold = cold(n) * distcoef
532 csrbnew = cnew(n) * distcoef
533 else if (this%isrb == 2)
then
536 else if (this%isrb == 3)
then
542 this%decayslast(n), &
543 kiter, csrbold, csrbnew,
delt)
544 this%decayslast(n) = decay_rate
545 rrhs = decay_rate * volfracm * rhobm * swnew * vcell
551 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
552 rhs(n) = rhs(n) + rrhs
565 subroutine mst_cq(this, nodes, cnew, cold, flowja)
569 integer(I4B),
intent(in) :: nodes
570 real(DP),
intent(in),
dimension(nodes) :: cnew
571 real(DP),
intent(in),
dimension(nodes) :: cold
572 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
576 call this%mst_cq_sto(nodes, cnew, cold, flowja)
579 if (this%idcy /= 0)
then
580 call this%mst_cq_dcy(nodes, cnew, cold, flowja)
584 if (this%isrb /= 0)
then
585 call this%mst_cq_srb(nodes, cnew, cold, flowja)
589 if (this%isrb /= 0 .and. this%idcy /= 0)
then
590 call this%mst_cq_dcy_srb(nodes, cnew, cold, flowja)
607 integer(I4B),
intent(in) :: nodes
608 real(DP),
intent(in),
dimension(nodes) :: cnew
609 real(DP),
intent(in),
dimension(nodes) :: cold
610 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
613 integer(I4B) :: idiag
616 real(DP) :: vnew, vold
617 real(DP) :: hhcof, rrhs
624 this%ratesto(n) =
dzero
627 if (this%ibound(n) <= 0) cycle
630 vnew = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) * &
631 this%fmi%gwfsat(n) * this%thetam(n)
633 if (this%fmi%igwfstrgss /= 0) vold = vold + this%fmi%gwfstrgss(n) *
delt
634 if (this%fmi%igwfstrgsy /= 0) vold = vold + this%fmi%gwfstrgsy(n) *
delt
638 rrhs = -vold * tled * cold(n)
639 rate = hhcof * cnew(n) - rrhs
640 this%ratesto(n) = rate
641 idiag = this%dis%con%ia(n)
642 flowja(idiag) = flowja(idiag) + rate
659 integer(I4B),
intent(in) :: nodes
660 real(DP),
intent(in),
dimension(nodes) :: cnew
661 real(DP),
intent(in),
dimension(nodes) :: cold
662 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
665 integer(I4B) :: idiag
668 real(DP) :: hhcof, rrhs
670 real(DP) :: decay_rate
678 this%ratedcy(n) =
dzero
679 if (this%ibound(n) <= 0) cycle
682 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
683 swtpdt = this%fmi%gwfsat(n)
689 if (this%idcy == 1)
then
690 hhcof = -this%decay(n) * vcell * swtpdt * this%thetam(n)
691 elseif (this%idcy == 2)
then
693 0, cold(n), cnew(n),
delt)
694 rrhs = decay_rate * vcell * swtpdt * this%thetam(n)
696 rate = hhcof * cnew(n) - rrhs
697 this%ratedcy(n) = rate
698 idiag = this%dis%con%ia(n)
699 flowja(idiag) = flowja(idiag) + rate
717 integer(I4B),
intent(in) :: nodes
718 real(DP),
intent(in),
dimension(nodes) :: cnew
719 real(DP),
intent(in),
dimension(nodes) :: cold
720 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
723 integer(I4B) :: idiag
726 real(DP) :: swt, swtpdt
740 this%ratesrb(n) =
dzero
743 if (this%ibound(n) <= 0) cycle
746 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
747 swtpdt = this%fmi%gwfsat(n)
748 swt = this%fmi%gwfsatold(n,
delt)
749 volfracm = this%get_volfracm(n)
750 rhobm = this%bulk_density(n)
751 const1 = this%distcoef(n)
753 if (this%isrb > 1) const2 = this%sp2(n)
754 call mst_srb_term(this%isrb, volfracm, rhobm, vcell, tled, cnew(n), &
755 cold(n), swtpdt, swt, const1, const2, &
757 this%ratesrb(n) = rate
758 idiag = this%dis%con%ia(n)
759 flowja(idiag) = flowja(idiag) + rate
777 integer(I4B),
intent(in) :: nodes
778 real(DP),
intent(in),
dimension(nodes) :: cnew
779 real(DP),
intent(in),
dimension(nodes) :: cold
780 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
783 integer(I4B) :: idiag
785 real(DP) :: hhcof, rrhs
795 real(DP) :: decay_rate
802 this%ratedcys(n) =
dzero
805 if (this%ibound(n) <= 0) cycle
810 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
811 swnew = this%fmi%gwfsat(n)
812 distcoef = this%distcoef(n)
813 volfracm = this%get_volfracm(n)
814 rhobm = this%bulk_density(n)
815 term = this%decay_sorbed(n) * volfracm * rhobm * swnew * vcell
818 if (this%idcy == 1)
then
820 if (this%isrb == 1)
then
824 hhcof = -term * distcoef
825 else if (this%isrb == 2)
then
830 else if (this%isrb == 3)
then
836 elseif (this%idcy == 2)
then
840 if (distcoef >
dzero)
then
841 if (this%isrb == 1)
then
842 csrbold = cold(n) * distcoef
843 csrbnew = cnew(n) * distcoef
844 else if (this%isrb == 2)
then
847 else if (this%isrb == 3)
then
852 this%decayslast(n), &
853 0, csrbold, csrbnew,
delt)
854 rrhs = decay_rate * volfracm * rhobm * swnew * vcell
859 rate = hhcof * cnew(n) - rrhs
860 this%ratedcys(n) = rate
861 idiag = this%dis%con%ia(n)
862 flowja(idiag) = flowja(idiag) + rate
875 subroutine mst_bd(this, isuppress_output, model_budget)
881 integer(I4B),
intent(in) :: isuppress_output
882 type(
budgettype),
intent(inout) :: model_budget
889 call model_budget%addentry(rin, rout,
delt,
budtxt(1), &
890 isuppress_output, rowlabel=this%packName)
893 if (this%idcy /= 0)
then
895 call model_budget%addentry(rin, rout,
delt,
budtxt(2), &
896 isuppress_output, rowlabel=this%packName)
900 if (this%isrb /= 0)
then
902 call model_budget%addentry(rin, rout,
delt,
budtxt(3), &
903 isuppress_output, rowlabel=this%packName)
907 if (this%isrb /= 0 .and. this%idcy /= 0)
then
909 call model_budget%addentry(rin, rout,
delt,
budtxt(4), &
910 isuppress_output, rowlabel=this%packName)
925 integer(I4B),
intent(in) :: icbcfl
926 integer(I4B),
intent(in) :: icbcun
928 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)
955 if (this%idcy /= 0) &
956 call this%dis%record_array(this%ratedcy, this%iout, iprint, -ibinun, &
957 budtxt(2), cdatafmp, nvaluesp, &
958 nwidthp, editdesc, dinact)
961 if (this%isrb /= 0) &
962 call this%dis%record_array(this%ratesrb, this%iout, iprint, -ibinun, &
963 budtxt(3), cdatafmp, nvaluesp, &
964 nwidthp, editdesc, dinact)
967 if (this%isrb /= 0 .and. this%idcy /= 0) &
968 call this%dis%record_array(this%ratedcys, this%iout, iprint, -ibinun, &
969 budtxt(4), cdatafmp, nvaluesp, &
970 nwidthp, editdesc, dinact)
989 if (this%inunit > 0)
then
1006 this%ibound => null()
1013 call this%NumericalPackageType%da()
1032 call this%NumericalPackageType%allocate_scalars()
1057 integer(I4B),
intent(in) :: nodes
1063 call mem_allocate(this%porosity, nodes,
'POROSITY', this%memoryPath)
1064 call mem_allocate(this%thetam, nodes,
'THETAM', this%memoryPath)
1065 call mem_allocate(this%volfracim, nodes,
'VOLFRACIM', this%memoryPath)
1066 call mem_allocate(this%ratesto, nodes,
'RATESTO', this%memoryPath)
1069 if (this%idcy == 0)
then
1070 call mem_allocate(this%ratedcy, 1,
'RATEDCY', this%memoryPath)
1071 call mem_allocate(this%decay, 1,
'DECAY', this%memoryPath)
1072 call mem_allocate(this%decaylast, 1,
'DECAYLAST', this%memoryPath)
1074 call mem_allocate(this%ratedcy, this%dis%nodes,
'RATEDCY', this%memoryPath)
1075 call mem_allocate(this%decay, nodes,
'DECAY', this%memoryPath)
1076 call mem_allocate(this%decaylast, nodes,
'DECAYLAST', this%memoryPath)
1078 if (this%idcy /= 0 .and. this%isrb /= 0)
then
1079 call mem_allocate(this%ratedcys, this%dis%nodes,
'RATEDCYS', &
1081 call mem_allocate(this%decayslast, this%dis%nodes,
'DECAYSLAST', &
1084 call mem_allocate(this%ratedcys, 1,
'RATEDCYS', this%memoryPath)
1085 call mem_allocate(this%decayslast, 1,
'DECAYSLAST', this%memoryPath)
1087 call mem_allocate(this%decay_sorbed, 1,
'DECAY_SORBED', &
1091 if (this%isrb == 0)
then
1092 call mem_allocate(this%bulk_density, 1,
'BULK_DENSITY', this%memoryPath)
1094 call mem_allocate(this%distcoef, 1,
'DISTCOEF', this%memoryPath)
1095 call mem_allocate(this%ratesrb, 1,
'RATESRB', this%memoryPath)
1097 call mem_allocate(this%bulk_density, nodes,
'BULK_DENSITY', this%memoryPath)
1098 call mem_allocate(this%distcoef, nodes,
'DISTCOEF', this%memoryPath)
1099 call mem_allocate(this%ratesrb, nodes,
'RATESRB', this%memoryPath)
1100 if (this%isrb == 1)
then
1103 call mem_allocate(this%sp2, nodes,
'SP2', this%memoryPath)
1109 this%porosity(n) =
dzero
1110 this%thetam(n) =
dzero
1111 this%volfracim(n) =
dzero
1112 this%ratesto(n) =
dzero
1114 do n = 1,
size(this%decay)
1115 this%decay(n) =
dzero
1116 this%ratedcy(n) =
dzero
1117 this%decaylast(n) =
dzero
1119 do n = 1,
size(this%bulk_density)
1120 this%bulk_density(n) =
dzero
1121 this%distcoef(n) =
dzero
1122 this%ratesrb(n) =
dzero
1124 do n = 1,
size(this%sp2)
1127 do n = 1,
size(this%ratedcys)
1128 this%ratedcys(n) =
dzero
1129 this%decayslast(n) =
dzero
1147 character(len=LINELENGTH) :: keyword, keyword2
1148 integer(I4B) :: ierr
1149 logical :: isfound, endOfBlock
1151 character(len=*),
parameter :: fmtisvflow = &
1152 "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
1153 &WHENEVER ICBCFL IS NOT ZERO.')"
1154 character(len=*),
parameter :: fmtisrb = &
1155 &
"(4x,'LINEAR SORPTION IS ACTIVE. ')"
1156 character(len=*),
parameter :: fmtfreundlich = &
1157 &
"(4x,'FREUNDLICH SORPTION IS ACTIVE. ')"
1158 character(len=*),
parameter :: fmtlangmuir = &
1159 &
"(4x,'LANGMUIR SORPTION IS ACTIVE. ')"
1160 character(len=*),
parameter :: fmtidcy1 = &
1161 &
"(4x,'FIRST-ORDER DECAY IS ACTIVE. ')"
1162 character(len=*),
parameter :: fmtidcy2 = &
1163 &
"(4x,'ZERO-ORDER DECAY IS ACTIVE. ')"
1166 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, blockrequired=.false., &
1167 supportopenclose=.true.)
1171 write (this%iout,
'(1x,a)')
'PROCESSING MOBILE STORAGE AND TRANSFER OPTIONS'
1173 call this%parser%GetNextLine(endofblock)
1174 if (endofblock)
exit
1175 call this%parser%GetStringCaps(keyword)
1176 select case (keyword)
1179 write (this%iout, fmtisvflow)
1180 case (
'SORBTION',
'SORPTION')
1182 call this%parser%GetStringCaps(keyword2)
1183 if (trim(adjustl(keyword2)) ==
'LINEAR') this%isrb = 1
1184 if (trim(adjustl(keyword2)) ==
'FREUNDLICH') this%isrb = 2
1185 if (trim(adjustl(keyword2)) ==
'LANGMUIR') this%isrb = 3
1186 select case (this%isrb)
1188 write (this%iout, fmtisrb)
1190 write (this%iout, fmtfreundlich)
1192 write (this%iout, fmtlangmuir)
1194 case (
'FIRST_ORDER_DECAY')
1196 write (this%iout, fmtidcy1)
1197 case (
'ZERO_ORDER_DECAY')
1199 write (this%iout, fmtidcy2)
1201 write (
errmsg,
'(a,a)')
'Unknown MST option: ', trim(keyword)
1203 call this%parser%StoreErrorUnit()
1206 write (this%iout,
'(1x,a)')
'END OF MOBILE STORAGE AND TRANSFER OPTIONS'
1225 character(len=LINELENGTH) :: keyword
1226 character(len=:),
allocatable :: line
1227 integer(I4B) :: istart, istop, lloc, ierr, n
1228 logical :: isfound, endOfBlock
1229 logical,
dimension(6) :: lname
1230 character(len=24),
dimension(6) :: aname
1233 data aname(1)/
' MOBILE DOMAIN POROSITY'/
1234 data aname(2)/
' BULK DENSITY'/
1235 data aname(3)/
'DISTRIBUTION COEFFICIENT'/
1236 data aname(4)/
' DECAY RATE'/
1237 data aname(5)/
' DECAY SORBED RATE'/
1238 data aname(6)/
' SECOND SORPTION PARAM'/
1245 call this%parser%GetBlock(
'GRIDDATA', isfound, ierr)
1247 write (this%iout,
'(1x,a)')
'PROCESSING GRIDDATA'
1249 call this%parser%GetNextLine(endofblock)
1250 if (endofblock)
exit
1251 call this%parser%GetStringCaps(keyword)
1252 call this%parser%GetRemainingLine(line)
1254 select case (keyword)
1256 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1257 this%parser%iuactive, this%porosity, &
1260 case (
'BULK_DENSITY')
1261 if (this%isrb == 0) &
1263 'BULK_DENSITY', trim(this%memoryPath))
1264 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1265 this%parser%iuactive, &
1266 this%bulk_density, aname(2))
1269 if (this%isrb == 0) &
1271 trim(this%memoryPath))
1272 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1273 this%parser%iuactive, this%distcoef, &
1277 if (this%idcy == 0) &
1279 trim(this%memoryPath))
1280 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1281 this%parser%iuactive, this%decay, &
1284 case (
'DECAY_SORBED')
1286 'DECAY_SORBED', trim(this%memoryPath))
1287 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1288 this%parser%iuactive, &
1289 this%decay_sorbed, aname(5))
1292 if (this%isrb < 2) &
1294 trim(this%memoryPath))
1295 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1296 this%parser%iuactive, this%sp2, &
1300 write (
errmsg,
'(a,a)')
'Unknown GRIDDATA tag: ', trim(keyword)
1302 call this%parser%StoreErrorUnit()
1305 write (this%iout,
'(1x,a)')
'END PROCESSING GRIDDATA'
1307 write (
errmsg,
'(a)')
'Required GRIDDATA block not found.'
1309 call this%parser%StoreErrorUnit()
1313 if (.not. lname(1))
then
1314 write (
errmsg,
'(a)')
'POROSITY not specified in GRIDDATA block.'
1319 if (this%isrb > 0)
then
1320 if (.not. lname(2))
then
1321 write (
errmsg,
'(a)')
'Sorption is active but BULK_DENSITY &
1322 ¬ specified. BULK_DENSITY must be specified in GRIDDATA block.'
1325 if (.not. lname(3))
then
1326 write (
errmsg,
'(a)')
'Sorption is active but distribution &
1327 &coefficient not specified. DISTCOEF must be specified in &
1331 if (this%isrb > 1)
then
1332 if (.not. lname(6))
then
1333 write (
errmsg,
'(a)')
'Freundlich or langmuir sorption is active &
1334 &but SP2 not specified. SP2 must be specified in &
1341 write (
warnmsg,
'(a)')
'Sorption is not active but &
1342 &BULK_DENSITY was specified. BULK_DENSITY will have no affect on &
1343 &simulation results.'
1345 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1348 write (
warnmsg,
'(a)')
'Sorption is not active but &
1349 &distribution coefficient was specified. DISTCOEF will have &
1350 &no affect on simulation results.'
1352 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1355 write (
warnmsg,
'(a)')
'Sorption is not active but &
1356 &SP2 was specified. SP2 will have &
1357 &no affect on simulation results.'
1359 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1364 if (this%idcy > 0)
then
1365 if (.not. lname(4))
then
1366 write (
errmsg,
'(a)')
'First or zero order decay is &
1367 &active but the first rate coefficient is not specified. DECAY &
1368 &must be specified in GRIDDATA block.'
1371 if (.not. lname(5))
then
1375 if (this%isrb > 0)
then
1376 write (
errmsg,
'(a)')
'DECAY_SORBED not provided in GRIDDATA &
1377 &block but decay and sorption are active. Specify DECAY_SORBED &
1378 &in GRIDDATA block.'
1384 write (
warnmsg,
'(a)')
'First- or zero-order decay &
1385 &is not active but decay was specified. DECAY will &
1386 &have no affect on simulation results.'
1388 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1391 write (
warnmsg,
'(a)')
'First- or zero-order decay &
1392 &is not active but DECAY_SORBED was specified. &
1393 &DECAY_SORBED will have no affect on simulation results.'
1395 write (this%iout,
'(1x,a)')
'WARNING. '//
warnmsg
1401 call this%parser%StoreErrorUnit()
1405 do n = 1,
size(this%porosity)
1406 this%thetam(n) = this%porosity(n)
1423 real(DP),
dimension(:),
intent(in) :: volfracim
1428 do n = 1, this%dis%nodes
1429 this%volfracim(n) = this%volfracim(n) + volfracim(n)
1434 do n = 1, this%dis%nodes
1435 this%thetam(n) = this%get_volfracm(n) * this%porosity(n)
1451 integer(I4B),
intent(in) :: node
1453 real(dp) :: volfracm
1455 volfracm =
done - this%volfracim(node)
1468 real(dp),
intent(in) :: conc
1469 real(dp),
intent(in) :: kf
1470 real(dp),
intent(in) :: a
1474 if (conc >
dzero)
then
1489 real(dp),
intent(in) :: conc
1490 real(dp),
intent(in) :: kl
1491 real(dp),
intent(in) :: sbar
1495 if (conc >
dzero)
then
1496 cbar = (kl * sbar * conc) / (
done + kl * conc)
1510 real(dp),
intent(in) :: conc
1511 real(dp),
intent(in) :: kf
1512 real(dp),
intent(in) :: a
1516 if (conc >
dzero)
then
1517 derv = kf * a * conc**(a -
done)
1531 real(dp),
intent(in) :: conc
1532 real(dp),
intent(in) :: kl
1533 real(dp),
intent(in) :: sbar
1537 if (conc >
dzero)
then
1538 derv = (kl * sbar) / (
done + kl * conc)**
dtwo
1554 cold, cnew, delt)
result(decay_rate)
1556 real(dp),
intent(in) :: decay_rate_usr
1557 real(dp),
intent(in) :: decay_rate_last
1558 integer(I4B),
intent(in) :: kiter
1559 real(dp),
intent(in) :: cold
1560 real(dp),
intent(in) :: cnew
1561 real(dp),
intent(in) :: delt
1563 real(dp) :: decay_rate
1566 if (decay_rate_usr <
dzero)
then
1569 decay_rate = decay_rate_usr
1575 if (kiter == 1)
then
1576 decay_rate = min(decay_rate_usr, cold / delt)
1578 decay_rate = decay_rate_last
1579 if (cnew <
dzero)
then
1580 decay_rate = decay_rate_last + cnew / delt
1581 else if (cnew > cold)
then
1582 decay_rate = decay_rate_last + cnew / delt
1584 decay_rate = min(decay_rate_usr, decay_rate)
1586 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
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dzero
real constant zero
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 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
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
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_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
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