32 character(len=LENFTYPE) ::
ftype =
'IST'
33 character(len=LENPACKAGENAME) ::
text =
' IMMOBILE DOMAIN'
35 character(len=LENBUDTXT),
dimension(NBDITEMS) ::
budtxt
36 data budtxt/
' STORAGE-AQUEOUS',
' STORAGE-SORBED', &
37 ' DECAY-AQUEOUS',
' DECAY-SORBED', &
55 integer(I4B),
pointer :: icimout => null()
56 integer(I4B),
pointer :: ibudgetout => null()
57 integer(I4B),
pointer :: ibudcsv => null()
58 integer(I4B),
pointer :: idcy => null()
59 integer(I4B),
pointer :: isrb => null()
60 integer(I4B),
pointer :: kiter => null()
61 real(dp),
dimension(:),
pointer,
contiguous :: cim => null()
62 real(dp),
dimension(:),
pointer,
contiguous :: cimnew => null()
63 real(dp),
dimension(:),
pointer,
contiguous :: cimold => null()
64 real(dp),
dimension(:),
pointer,
contiguous :: zetaim => null()
65 real(dp),
dimension(:),
pointer,
contiguous :: porosity => null()
66 real(dp),
dimension(:),
pointer,
contiguous :: volfrac => null()
67 real(dp),
dimension(:),
pointer,
contiguous :: bulk_density => null()
68 real(dp),
dimension(:),
pointer,
contiguous :: distcoef => null()
69 real(dp),
dimension(:),
pointer,
contiguous :: decay => null()
70 real(dp),
dimension(:),
pointer,
contiguous :: decaylast => null()
71 real(dp),
dimension(:),
pointer,
contiguous :: decayslast => null()
72 real(dp),
dimension(:),
pointer,
contiguous :: decay_sorbed => null()
73 real(dp),
dimension(:),
pointer,
contiguous :: strg => null()
74 real(dp),
dimension(2, NBDITEMS) :: budterm
107 subroutine ist_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
110 class(
bndtype),
pointer :: packobj
111 integer(I4B),
intent(in) :: id
112 integer(I4B),
intent(in) :: ibcnum
113 integer(I4B),
intent(in) :: inunit
114 integer(I4B),
intent(in) :: iout
115 character(len=*),
intent(in) :: namemodel
116 character(len=*),
intent(in) :: pakname
127 call packobj%set_names(ibcnum, namemodel, pakname,
ftype)
131 call packobj%allocate_scalars()
134 call packobj%pack_initialize()
137 packobj%inunit = inunit
140 packobj%ibcnum = ibcnum
167 call this%ist_allocate_arrays()
171 call this%ocd%init_dbl(
'CIM', this%cimnew, this%dis,
'PRINT LAST ', &
172 'COLUMNS 10 WIDTH 11 DIGITS 4 GENERAL ', &
176 call this%read_data()
179 do n = 1, this%dis%nodes
180 this%cimnew(n) = this%cim(n)
184 call this%mst%addto_volfracim(this%volfrac)
187 call budget_cr(this%budget, this%memoryPath)
188 call this%budget%budget_df(
nbditems,
'MASS',
'M', bdzone=this%packName)
189 call this%budget%set_ibudcsv(this%ibudcsv)
193 if (this%idcy /= this%mst%idcy)
then
194 call store_error(
'DECAY must be activated consistently between the &
195 &MST and IST Packages. Activate or deactivate DECAY for both &
198 if (this%isrb /= this%mst%isrb)
then
199 call store_error(
'Sorption is active for the IST Package but it is not &
200 &compatible with the sorption option selected for the MST Package. &
201 &If sorption is active for the IST Package, then SORPTION LINEAR must &
202 &be specified in the options block of the MST Package.')
205 call this%parser%StoreErrorUnit()
242 call this%BndType%bnd_ad()
250 do n = 1, this%dis%nodes
251 this%cimold(n) = this%cimnew(n)
254 do n = 1, this%dis%nodes
255 this%cimnew(n) = this%cimold(n)
267 subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln)
272 real(DP),
dimension(:),
intent(inout) :: rhs
273 integer(I4B),
dimension(:),
intent(in) :: ia
274 integer(I4B),
dimension(:),
intent(in) :: idxglo
277 integer(I4B) :: n, idiag
279 real(DP) :: hhcof, rrhs
280 real(DP) :: swt, swtpdt
285 real(DP) :: volfracim
287 real(DP) :: lambda1im
288 real(DP) :: lambda2im
293 real(DP) :: cimsrbold
294 real(DP) :: cimsrbnew
295 real(DP),
dimension(9) :: ddterm
299 this%kiter = this%kiter + 1
302 do n = 1, this%dis%nodes
305 if (this%ibound(n) <= 0) cycle
308 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
309 swtpdt = this%fmi%gwfsat(n)
310 swt = this%fmi%gwfsatold(n,
delt)
311 thetaim = this%get_thetaim(n)
315 zetaim = this%zetaim(n)
327 if (this%idcy == 1) lambda1im = this%decay(n)
328 if (this%idcy == 2)
then
330 this%kiter, this%cimold(n), &
331 this%cimnew(n),
delt)
332 this%decaylast(n) = gamma1im
336 if (this%isrb > 0)
then
337 kd = this%distcoef(n)
338 volfracim = this%volfrac(n)
339 rhobim = this%bulk_density(n)
340 if (this%idcy == 1) lambda2im = this%decay_sorbed(n)
341 if (this%idcy == 2)
then
342 cimsrbold = this%cimold(n) * kd
343 cimsrbnew = this%cimnew(n) * kd
345 this%decayslast(n), &
346 this%kiter, cimsrbold, &
348 this%decayslast(n) = gamma2im
354 volfracim, rhobim, kd, lambda1im, lambda2im, &
355 gamma1im, gamma2im, zetaim, ddterm, f)
356 cimold = this%cimold(n)
360 call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
361 rhs(n) = rhs(n) + rrhs
380 real(DP),
dimension(:),
intent(in) :: x
381 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
382 integer(I4B),
optional,
intent(in) :: iadv
384 integer(I4B) :: idiag
387 real(DP) :: swt, swtpdt
388 real(DP) :: hhcof, rrhs
393 real(DP) :: volfracim
395 real(DP) :: lambda1im
396 real(DP) :: lambda2im
402 real(DP) :: cimsrbold
403 real(DP) :: cimsrbnew
404 real(DP),
dimension(9) :: ddterm
408 this%budterm(:, :) =
dzero
411 do n = 1, this%dis%nodes
416 if (this%ibound(n) > 0)
then
419 vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
420 swtpdt = this%fmi%gwfsat(n)
421 swt = this%fmi%gwfsatold(n,
delt)
422 thetaim = this%get_thetaim(n)
425 zetaim = this%zetaim(n)
438 if (this%idcy == 1) lambda1im = this%decay(n)
439 if (this%idcy == 2)
then
441 this%cimold(n), this%cimnew(n),
delt)
443 if (this%isrb > 0)
then
444 kd = this%distcoef(n)
445 volfracim = this%volfrac(n)
446 rhobim = this%bulk_density(n)
447 if (this%idcy == 1) lambda2im = this%decay_sorbed(n)
448 if (this%idcy == 2)
then
449 cimsrbold = this%cimold(n) * kd
450 cimsrbnew = this%cimnew(n) * kd
452 this%decayslast(n), &
460 volfracim, rhobim, kd, lambda1im, lambda2im, &
461 gamma1im, gamma2im, zetaim, ddterm, f)
462 cimold = this%cimold(n)
466 rate = hhcof * x(n) - rrhs
478 idiag = this%dis%con%ia(n)
479 flowja(idiag) = flowja(idiag) + rate
482 this%cimnew(n) = cimnew
500 type(
budgettype),
intent(inout) :: model_budget
504 integer(I4B) :: isuppress_output
507 call model_budget%addentry(ratin, ratout,
delt, this%text, &
508 isuppress_output, this%packName)
523 integer(I4B),
intent(in) :: icbcfl
524 integer(I4B),
intent(in) :: ibudfl
525 integer(I4B),
intent(in) :: icbcun
526 integer(I4B),
dimension(:),
optional,
intent(in) :: imap
529 integer(I4B) :: ibinun
530 integer(I4B) :: nbound
535 if (this%ipakcb < 0)
then
537 elseif (this%ipakcb == 0)
then
542 if (icbcfl == 0) ibinun = 0
547 if (ibinun /= 0)
then
548 nbound = this%dis%nodes
550 call this%dis%record_srcdst_list_header(this%text, this%name_model, &
551 this%name_model, this%name_model, &
552 this%packName, naux, this%auxname, &
553 ibinun, nbound, this%iout)
557 do n = 1, this%dis%nodes
561 if (this%ibound(n) > 0)
then
568 if (ibinun /= 0)
then
569 call this%dis%record_mf6_list_entry(ibinun, n, n, rate, &
570 naux, this%auxvar(:, n), &
589 integer(I4B),
intent(in) :: idvsave
590 integer(I4B),
intent(in) :: idvprint
592 integer(I4B) :: ipflg
593 integer(I4B) :: ibinun
600 if (idvsave == 0) ibinun = 0
601 if (ibinun /= 0)
then
603 iprint_opt=0, isav_opt=ibinun)
607 if (idvprint /= 0)
then
609 iprint_opt=idvprint, isav_opt=0)
625 integer(I4B),
intent(in) :: kstp
626 integer(I4B),
intent(in) :: kper
627 integer(I4B),
intent(in) :: iout
628 integer(I4B),
intent(in) :: ibudfl
630 integer(I4B) :: isuppress_output = 0
633 call this%budget%reset()
634 call this%budget%addentry(this%budterm,
delt,
budtxt, isuppress_output)
637 call this%budget%finalize_step(
delt)
638 if (ibudfl /= 0)
then
639 call this%budget%budget_ot(kstp, kper, iout)
643 call this%budget%writecsv(
totim)
659 if (this%inunit > 0)
then
686 call this%budget%budget_da()
687 deallocate (this%budget)
688 call this%ocd%ocd_da()
689 deallocate (this%ocd)
692 call this%BndType%bnd_da()
712 call this%BndType%allocate_scalars()
715 call mem_allocate(this%icimout,
'ICIMOUT', this%memoryPath)
716 call mem_allocate(this%ibudgetout,
'IBUDGETOUT', this%memoryPath)
717 call mem_allocate(this%ibudcsv,
'IBUDCSV', this%memoryPath)
753 call this%BndType%allocate_arrays()
756 call mem_allocate(this%strg, this%dis%nodes,
'STRG', this%memoryPath)
757 call mem_allocate(this%cim, this%dis%nodes,
'CIM', this%memoryPath)
758 call mem_allocate(this%cimnew, this%dis%nodes,
'CIMNEW', this%memoryPath)
759 call mem_allocate(this%cimold, this%dis%nodes,
'CIMOLD', this%memoryPath)
760 call mem_allocate(this%porosity, this%dis%nodes,
'POROSITY', this%memoryPath)
761 call mem_allocate(this%zetaim, this%dis%nodes,
'ZETAIM', this%memoryPath)
762 call mem_allocate(this%volfrac, this%dis%nodes,
'VOLFRAC', this%memoryPath)
763 if (this%isrb == 0)
then
764 call mem_allocate(this%bulk_density, 1,
'BULK_DENSITY', this%memoryPath)
765 call mem_allocate(this%distcoef, 1,
'DISTCOEF', this%memoryPath)
767 call mem_allocate(this%bulk_density, this%dis%nodes,
'BULK_DENSITY', &
769 call mem_allocate(this%distcoef, this%dis%nodes,
'DISTCOEF', &
772 if (this%idcy == 0)
then
773 call mem_allocate(this%decay, 1,
'DECAY', this%memoryPath)
774 call mem_allocate(this%decaylast, 1,
'DECAYLAST', this%memoryPath)
776 call mem_allocate(this%decay, this%dis%nodes,
'DECAY', this%memoryPath)
777 call mem_allocate(this%decaylast, this%dis%nodes,
'DECAYLAST', &
780 if (this%isrb == 0 .and. this%idcy == 0)
then
781 call mem_allocate(this%decayslast, 1,
'DECAYSLAST', this%memoryPath)
783 call mem_allocate(this%decayslast, this%dis%nodes,
'DECAYSLAST', &
786 call mem_allocate(this%decay_sorbed, 1,
'DECAY_SORBED', this%memoryPath)
789 do n = 1, this%dis%nodes
792 this%cimnew(n) =
dzero
793 this%cimold(n) =
dzero
794 this%porosity(n) =
dzero
795 this%zetaim(n) =
dzero
796 this%volfrac(n) =
dzero
798 do n = 1,
size(this%decay)
799 this%decay(n) =
dzero
800 this%decaylast(n) =
dzero
802 do n = 1,
size(this%decayslast)
803 this%decayslast(n) =
dzero
807 this%ocd%dis => this%dis
827 character(len=LINELENGTH) :: errmsg, keyword
828 character(len=LINELENGTH) :: fname
829 character(len=:),
allocatable :: keyword2
831 logical :: isfound, endOfBlock
834 character(len=*),
parameter :: fmtisvflow = &
835 "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
836 &WHENEVER ICBCFL IS NOT ZERO.')"
837 character(len=*),
parameter :: fmtisrb = &
838 &
"(4x,'LINEAR SORPTION IS SELECTED. ')"
839 character(len=*),
parameter :: fmtidcy1 = &
840 &
"(4x,'FIRST-ORDER DECAY IS ACTIVE. ')"
841 character(len=*),
parameter :: fmtidcy2 = &
842 &
"(4x,'ZERO-ORDER DECAY IS ACTIVE. ')"
843 character(len=*),
parameter :: fmtistbin = &
844 "(4x, 'IST ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, &
845 &/4x, 'OPENED ON UNIT: ', I0)"
848 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, blockrequired=.false., &
849 supportopenclose=.true.)
853 write (this%iout,
'(1x,a)')
'PROCESSING IMMOBILE STORAGE AND TRANSFER &
856 call this%parser%GetNextLine(endofblock)
858 call this%parser%GetStringCaps(keyword)
859 select case (keyword)
862 write (this%iout, fmtisvflow)
864 call this%parser%GetRemainingLine(keyword2)
865 call this%ocd%set_option(keyword2, this%inunit, this%iout)
867 call this%parser%GetStringCaps(keyword)
868 if (keyword ==
'FILEOUT')
then
869 call this%parser%GetString(fname)
871 call openfile(this%ibudgetout, this%iout, fname,
'DATA(BINARY)', &
873 write (this%iout, fmtistbin)
'BUDGET', trim(adjustl(fname)), &
878 &be followed by FILEOUT')
881 call this%parser%GetStringCaps(keyword)
882 if (keyword ==
'FILEOUT')
then
883 call this%parser%GetString(fname)
885 call openfile(this%ibudcsv, this%iout, fname,
'CSV', &
886 filstat_opt=
'REPLACE')
887 write (this%iout, fmtistbin)
'BUDGET CSV', trim(adjustl(fname)), &
890 call store_error(
'Optional BUDGETCSV keyword must be followed by &
893 case (
'SORBTION',
'SORPTION')
895 write (this%iout, fmtisrb)
896 case (
'FIRST_ORDER_DECAY')
898 write (this%iout, fmtidcy1)
899 case (
'ZERO_ORDER_DECAY')
901 write (this%iout, fmtidcy2)
903 write (errmsg,
'(a,a)')
'Unknown IST option: ', &
906 call this%parser%StoreErrorUnit()
909 write (this%iout,
'(1x,a)')
'END OF IMMOBILE STORAGE AND TRANSFER &
945 character(len=LINELENGTH) :: errmsg, keyword
946 character(len=:),
allocatable :: line
947 integer(I4B) :: istart, istop, lloc, ierr
948 logical :: isfound, endOfBlock
949 logical,
dimension(8) :: lname
950 character(len=24),
dimension(8) :: aname
953 data aname(1)/
' BULK DENSITY'/
954 data aname(2)/
'DISTRIBUTION COEFFICIENT'/
955 data aname(3)/
' DECAY RATE'/
956 data aname(4)/
' DECAY SORBED RATE'/
957 data aname(5)/
' INITIAL IMMOBILE CONC'/
958 data aname(6)/
' FIRST ORDER TRANS RATE'/
959 data aname(7)/
'IMMOBILE DOMAIN POROSITY'/
960 data aname(8)/
'IMMOBILE VOLUME FRACTION'/
967 call this%parser%GetBlock(
'GRIDDATA', isfound, ierr)
969 write (this%iout,
'(1x,a)')
'PROCESSING GRIDDATA'
971 call this%parser%GetNextLine(endofblock)
973 call this%parser%GetStringCaps(keyword)
974 call this%parser%GetRemainingLine(line)
976 select case (keyword)
977 case (
'BULK_DENSITY')
978 if (this%isrb == 0) &
980 'BULK_DENSITY', trim(this%memoryPath))
981 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
982 this%parser%iuactive, &
983 this%bulk_density, aname(1))
986 if (this%isrb == 0) &
988 trim(this%memoryPath))
989 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
990 this%parser%iuactive, this%distcoef, &
994 if (this%idcy == 0) &
996 trim(this%memoryPath))
997 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
998 this%parser%iuactive, this%decay, &
1001 case (
'DECAY_SORBED')
1003 'DECAY_SORBED', trim(this%memoryPath))
1004 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1005 this%parser%iuactive, &
1006 this%decay_sorbed, aname(4))
1009 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1010 this%parser%iuactive, this%cim, &
1014 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1015 this%parser%iuactive, this%zetaim, &
1019 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1020 this%parser%iuactive, this%porosity, &
1024 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
1025 this%parser%iuactive, this%volfrac, &
1029 write (errmsg,
'(a)') &
1030 'THETAIM is no longer supported. See Chapter 9 in &
1031 &mf6suptechinfo.pdf for revised parameterization of mobile and &
1032 &immobile domain simulations.'
1034 call this%parser%StoreErrorUnit()
1036 write (errmsg,
'(a,a)')
'Unknown GRIDDATA tag: ', trim(keyword)
1038 call this%parser%StoreErrorUnit()
1041 write (this%iout,
'(1x,a)')
'END PROCESSING GRIDDATA'
1043 write (errmsg,
'(a)')
'Required GRIDDATA block not found.'
1045 call this%parser%StoreErrorUnit()
1049 if (this%isrb > 0)
then
1050 if (.not. lname(1))
then
1051 write (errmsg,
'(a)')
'Sorption is active but BULK_DENSITY &
1052 ¬ specified. BULK_DENSITY must be specified in griddata block.'
1055 if (.not. lname(2))
then
1056 write (errmsg,
'(a)')
'Sorption is active but distribution &
1057 &coefficient not specified. DISTCOEF must be specified in &
1063 write (this%iout,
'(1x,a)')
'Warning. Sorption is not active but &
1064 &BULK_DENSITY was specified. BULK_DENSITY will have no affect on &
1065 &simulation results.'
1068 write (this%iout,
'(1x,a)')
'Warning. Sorption is not active but &
1069 &distribution coefficient was specified. DISTCOEF will have &
1070 &no affect on simulation results.'
1075 if (this%idcy > 0)
then
1076 if (.not. lname(3))
then
1077 write (errmsg,
'(a)')
'First- or zero-order decay is &
1078 &active but the first rate coefficient was not specified. &
1079 &Decay must be specified in GRIDDATA block.'
1082 if (.not. lname(4))
then
1086 if (this%isrb > 0)
then
1087 write (errmsg,
'(a)')
'DECAY_SORBED not provided in GRIDDATA &
1088 &block but decay and sorption are active. Specify DECAY_SORBED &
1089 &in GRIDDATA block.'
1095 write (this%iout,
'(1x,a)')
'Warning. First- or zero-order decay &
1096 &is not active but DECAY was specified. DECAY will &
1097 &have no affect on simulation results.'
1100 write (this%iout,
'(1x,a)')
'Warning. First- or zero-order decay &
1101 &is not active but DECAY_SORBED was specified. &
1102 &DECAY_SORBED will have no affect on simulation &
1109 if (.not. lname(5))
then
1110 write (this%iout,
'(1x,a)')
'Warning. Dual domain is active but &
1111 &initial immobile domain concentration was not specified. &
1112 &Setting CIM to zero.'
1114 if (.not. lname(6))
then
1115 write (errmsg,
'(a)')
'Dual domain is active but dual &
1116 &domain mass transfer rate (ZETAIM) was not specified. ZETAIM &
1117 &must be specified in GRIDDATA block.'
1120 if (.not. lname(7))
then
1121 write (errmsg,
'(a)')
'Dual domain is active but &
1122 &immobile domain POROSITY was not specified. POROSITY &
1123 &must be specified in GRIDDATA block.'
1126 if (.not. lname(8))
then
1127 write (errmsg,
'(a)')
'Dual domain is active but &
1128 &immobile domain VOLFRAC was not specified. VOLFRAC &
1129 &must be specified in GRIDDATA block. This is a new &
1130 &requirement for MODFLOW versions later than version &
1137 call this%parser%StoreErrorUnit()
1153 integer(I4B),
intent(in) :: node
1157 thetaim = this%volfrac(node) * this%porosity(node)
1172 volfracim, rhobim, kd, lambda1im, lambda2im, &
1173 gamma1im, gamma2im, zetaim, ddterm, f)
1175 real(DP),
intent(in) :: thetaim
1176 real(DP),
intent(in) :: vcell
1177 real(DP),
intent(in) :: delt
1178 real(DP),
intent(in) :: swtpdt
1179 real(DP),
intent(in) :: volfracim
1180 real(DP),
intent(in) :: rhobim
1181 real(DP),
intent(in) :: kd
1182 real(DP),
intent(in) :: lambda1im
1183 real(DP),
intent(in) :: lambda2im
1184 real(DP),
intent(in) :: gamma1im
1185 real(DP),
intent(in) :: gamma2im
1186 real(DP),
intent(in) :: zetaim
1187 real(DP),
dimension(:),
intent(inout) :: ddterm
1188 real(DP),
intent(inout) :: f
1199 ddterm(1) = thetaim * vcell * tled
1200 ddterm(2) = thetaim * vcell * tled
1201 ddterm(3) = volfracim * rhobim * vcell * kd * tled
1202 ddterm(4) = volfracim * rhobim * vcell * kd * tled
1203 ddterm(5) = thetaim * lambda1im * vcell
1204 ddterm(6) = lambda2im * volfracim * rhobim * kd * vcell
1205 ddterm(7) = thetaim * gamma1im * vcell
1206 ddterm(8) = gamma2im * volfracim * rhobim * vcell
1207 ddterm(9) = vcell * swtpdt * zetaim
1210 f = ddterm(1) + ddterm(3) + ddterm(5) + ddterm(6) + ddterm(9)
1224 real(DP),
dimension(:),
intent(in) :: ddterm
1225 real(DP),
intent(in) :: f
1226 real(DP),
intent(in) :: cimold
1227 real(DP),
intent(inout) :: hcof
1228 real(DP),
intent(inout) :: rhs
1231 hcof = ddterm(9)**2 / f - ddterm(9)
1235 rhs = (ddterm(2) + ddterm(4)) * cimold - ddterm(7) - ddterm(8)
1236 rhs = rhs * ddterm(9) / f
1250 real(dp),
dimension(:),
intent(in) :: ddterm
1251 real(dp),
intent(in) :: f
1252 real(dp),
intent(in) :: cimold
1253 real(dp),
intent(in) :: cnew
1259 cimnew = (ddterm(2) + ddterm(4)) * cimold + ddterm(9) * cnew - ddterm(7) &
1276 real(DP),
dimension(:, :),
intent(inout) :: budterm
1277 real(DP),
dimension(:),
intent(in) :: ddterm
1278 real(DP),
intent(in) :: cimnew
1279 real(DP),
intent(in) :: cimold
1280 real(DP),
intent(in) :: cnew
1281 integer(I4B),
intent(in) :: idcy
1288 rate = -ddterm(1) * cimnew + ddterm(2) * cimold
1289 if (rate >
dzero)
then
1290 budterm(1, i) = budterm(1, i) + rate
1292 budterm(2, i) = budterm(2, i) - rate
1297 rate = -ddterm(3) * cimnew + ddterm(4) * cimold
1298 if (rate >
dzero)
then
1299 budterm(1, i) = budterm(1, i) + rate
1301 budterm(2, i) = budterm(2, i) - rate
1308 rate = -ddterm(5) * cimnew
1309 else if (idcy == 2)
then
1314 if (rate >
dzero)
then
1315 budterm(1, i) = budterm(1, i) + rate
1317 budterm(2, i) = budterm(2, i) - rate
1323 rate = -ddterm(6) * cimnew
1324 else if (idcy == 2)
then
1329 if (rate >
dzero)
then
1330 budterm(1, i) = budterm(1, i) + rate
1332 budterm(2, i) = budterm(2, i) - rate
1337 rate = ddterm(9) * cnew - ddterm(9) * cimnew
1338 if (rate >
dzero)
then
1339 budterm(1, i) = budterm(1, i) + rate
1341 budterm(2, i) = budterm(2, i) - rate
This module contains the base boundary package.
This module contains the BudgetModule.
subroutine, public budget_cr(this, name_model)
@ brief Create a new budget object
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
integer(i4b), parameter lenpackagename
maximum length of the package name
real(dp), parameter dhnoflo
real no flow constant
integer(i4b), parameter lenftype
maximum length of a package type (DIS, WEL, OC, etc.)
real(dp), parameter dzero
real constant zero
integer(i4b), parameter lenbudtxt
maximum length of a budget component names
real(dp), parameter done
real constant 1
– @ brief Immobile Storage and Transfer (IST) Module
subroutine ist_da(this)
@ brief Deallocate package memory
subroutine ist_rp(this)
@ brief Read and prepare method for package
real(dp) function get_ddconc(ddterm, f, cimold, cnew)
@ brief Calculate the immobile domain concentration
subroutine, public ist_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, fmi, mst)
@ brief Create a new package object
real(dp) function get_thetaim(this, node)
@ brief Return thetaim
subroutine get_hcofrhs(ddterm, f, cimold, hcof, rhs)
@ brief Calculate the hcof and rhs terms for immobile domain
integer(i4b), parameter nbditems
subroutine ist_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
@ brief Output model flow terms.
subroutine ist_read_dimensions(this)
@ brief Read dimensions for package
subroutine allocate_scalars(this)
@ brief Allocate package scalars
character(len=lenftype) ftype
subroutine ist_ar(this)
@ brief Allocate and read method for package
subroutine ist_ot_bdsummary(this, kstp, kper, iout, ibudfl)
@ brief Output IST package budget summary.
subroutine ist_ot_dv(this, idvsave, idvprint)
@ brief Output immobile domain concentration.
subroutine get_ddterm(thetaim, vcell, delt, swtpdt, volfracim, rhobim, kd, lambda1im, lambda2im, gamma1im, gamma2im, zetaim, ddterm, f)
@ brief Calculate immobile domain equation terms
subroutine ist_allocate_arrays(this)
@ brief Allocate package arrays
subroutine read_data(this)
@ brief Read data for package
subroutine ist_cq(this, x, flowja, iadv)
@ brief Calculate package flows.
subroutine ist_ad(this)
@ brief Advance the ist package
subroutine accumulate_budterm(budterm, ddterm, cimnew, cimold, cnew, idcy)
@ brief Calculate the immobile domain budget terms
subroutine ist_bd(this, model_budget)
@ brief Add package flows to model budget.
subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln)
@ brief Fill coefficient method for package
character(len=lenbudtxt), dimension(nbditems) budtxt
subroutine read_options(this)
@ brief Read options for package
character(len=lenpackagename) text
– @ brief Mobile Storage and Transfer (MST) Module
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
This module defines variable data types.
This module contains the OutputControlDataModule.
subroutine, public ocd_cr(ocdobj)
@ brief Create OutputControlDataType
This module contains simulation methods.
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.
integer(i4b) ifailedstepretry
current retry for this time step
logical(lgp), pointer, public endofperiod
flag indicating end of stress period
real(dp), pointer, public totim
time relative to start of simulation
integer(i4b), pointer, public kstp
current time step number
real(dp), pointer, public delt
length of the current time step
Derived type for the Budget object.
@ brief Immobile storage and transfer
@ brief Mobile storage and transfer
@ brief OutputControlDataType