29 character(len=LENBUDTXT),
dimension(2) ::
budtxt = & !< text labels for budget terms
30 &[
' STO-SS',
' STO-SY']
33 integer(I4B),
pointer :: istor_coef => null()
34 integer(I4B),
pointer :: iconf_ss => null()
35 integer(I4B),
pointer :: iorig_ss => null()
36 integer(I4B),
pointer :: iss => null()
37 integer(I4B),
pointer :: iusesy => null()
38 integer(I4B),
dimension(:),
pointer,
contiguous :: iconvert => null()
39 real(dp),
dimension(:),
pointer,
contiguous :: ss => null()
40 real(dp),
dimension(:),
pointer,
contiguous :: sy => null()
41 real(dp),
dimension(:),
pointer,
contiguous :: strgss => null()
42 real(dp),
dimension(:),
pointer,
contiguous :: strgsy => null()
43 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
44 real(dp),
pointer :: satomega => null()
45 integer(I4B),
pointer :: integratechanges => null()
46 integer(I4B),
pointer :: intvs => null()
48 real(dp),
dimension(:),
pointer,
contiguous,
private :: oldss => null()
49 real(dp),
dimension(:),
pointer,
contiguous,
private :: oldsy => null()
75 subroutine sto_cr(stoobj, name_model, inunit, iout)
78 character(len=*),
intent(in) :: name_model
79 integer(I4B),
intent(in) :: inunit
80 integer(I4B),
intent(in) :: iout
86 call stoobj%set_names(1, name_model,
'STO',
'STO')
89 call stoobj%allocate_scalars()
92 stoobj%inunit = inunit
96 call stoobj%parser%Initialize(stoobj%inunit, stoobj%iout)
114 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
117 character(len=*),
parameter :: fmtsto = &
118 "(1x,/1x,'STO -- STORAGE PACKAGE, VERSION 1, 5/19/2014', &
119 &' INPUT READ FROM UNIT ', i0, //)"
122 write (this%iout, fmtsto) this%inunit
126 this%ibound => ibound
132 call this%allocate_arrays(dis%nodes)
138 call this%read_options()
141 call this%read_data()
144 if (this%intvs /= 0)
then
145 call this%tvs%ar(this%dis)
165 logical :: isfound, readss, readsy, endOfBlock
166 character(len=16) :: css(0:1)
167 character(len=LINELENGTH) :: line, keyword
169 character(len=*),
parameter :: fmtlsp = &
170 &
"(1X,/1X,'REUSING ',A,' FROM LAST STRESS PERIOD')"
171 character(len=*),
parameter :: fmtblkerr = &
172 &
"('Error. Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
174 data css(0)/
' TRANSIENT'/
175 data css(1)/
' STEADY-STATE'/
179 if (this%integratechanges /= 0)
then
180 call this%save_old_ss_sy()
184 if (this%ionper <
kper)
then
187 call this%parser%GetBlock(
'PERIOD', isfound, ierr, &
188 supportopenclose=.true., &
189 blockrequired=.false.)
193 call this%read_check_ionper()
199 this%ionper =
nper + 1
202 call this%parser%GetCurrentLine(line)
203 write (
errmsg, fmtblkerr) adjustl(trim(line))
205 call this%parser%StoreErrorUnit()
216 if (this%ionper ==
kper)
then
217 write (this%iout,
'(//,1x,a)')
'PROCESSING STORAGE PERIOD DATA'
219 call this%parser%GetNextLine(endofblock)
221 call this%parser%GetStringCaps(keyword)
222 select case (keyword)
223 case (
'STEADY-STATE')
228 write (
errmsg,
'(a,a)')
'Unknown STORAGE data tag: ', &
231 call this%parser%StoreErrorUnit()
234 write (this%iout,
'(1x,a)')
'END PROCESSING STORAGE PERIOD DATA'
239 write (this%iout,
'(//1X,A,I0,A,A,/)') &
240 'STRESS PERIOD ',
kper,
' IS ', trim(adjustl(css(this%iss)))
243 if (this%intvs /= 0)
then
263 if (this%integratechanges /= 0 .and.
kstp > 1)
then
264 call this%save_old_ss_sy()
268 if (this%intvs /= 0)
then
281 subroutine sto_fc(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
286 integer(I4B),
intent(in) :: kiter
287 real(DP),
intent(in),
dimension(:) :: hold
288 real(DP),
intent(in),
dimension(:) :: hnew
290 integer(I4B),
intent(in),
dimension(:) :: idxglo
291 real(DP),
intent(inout),
dimension(:) :: rhs
294 integer(I4B) :: idiag
311 character(len=*),
parameter :: fmtsperror = &
312 &
"('DETECTED TIME STEP LENGTH OF ZERO. GWF STORAGE PACKAGE CANNOT BE ', &
313 &'USED UNLESS DELT IS NON-ZERO.')"
316 if (this%iss /= 0)
return
320 write (
errmsg, fmtsperror)
328 do n = 1, this%dis%nodes
329 idiag = this%dis%con%ia(n)
330 if (this%ibound(n) < 1) cycle
337 if (this%iconvert(n) == 0)
then
346 sc1 =
sscapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
349 if (this%integratechanges /= 0)
then
353 sc1old =
sscapacity(this%istor_coef, tp, bt, this%dis%area(n), &
355 rho1old = sc1old * tled
363 call ssterms(this%iconvert(n), this%iorig_ss, this%iconf_ss, tp, bt, &
364 rho1, rho1old, snnew, snold, hnew(n), hold(n), &
368 call matrix_sln%add_value_pos(idxglo(idiag), aterm)
369 rhs(n) = rhs(n) + rhsterm
372 if (this%iconvert(n) /= 0)
then
376 sc2 =
sycapacity(this%dis%area(n), this%sy(n))
379 if (this%integratechanges /= 0)
then
383 sc2old =
sycapacity(this%dis%area(n), this%oldsy(n))
384 rho2old = sc2old * tled
392 call syterms(tp, bt, rho2, rho2old, snnew, snold, &
396 call matrix_sln%add_value_pos(idxglo(idiag), aterm)
397 rhs(n) = rhs(n) + rhsterm
411 subroutine sto_fn(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
416 integer(I4B),
intent(in) :: kiter
417 real(DP),
intent(in),
dimension(:) :: hold
418 real(DP),
intent(in),
dimension(:) :: hnew
420 integer(I4B),
intent(in),
dimension(:) :: idxglo
421 real(DP),
intent(inout),
dimension(:) :: rhs
424 integer(I4B) :: idiag
440 if (this%iss /= 0)
return
446 do n = 1, this%dis%nodes
447 idiag = this%dis%con%ia(n)
448 if (this%ibound(n) <= 0) cycle
460 sc1 =
sscapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
461 sc2 =
sycapacity(this%dis%area(n), this%sy(n))
467 if (this%iconvert(n) /= 0)
then
473 if (this%iconf_ss == 0)
then
474 if (this%iorig_ss == 0)
then
475 drterm = -rho1 * derv * (h - bt) + rho1 * tthk * snnew * derv
477 drterm = -(rho1 * derv * h)
479 call matrix_sln%add_value_pos(idxglo(idiag), drterm)
480 rhs(n) = rhs(n) + drterm * h
486 if (snnew <
done)
then
488 if (snnew >
dzero)
then
489 rterm = -rho2 * tthk * snnew
490 drterm = -rho2 * tthk * derv
491 call matrix_sln%add_value_pos(idxglo(idiag), drterm + rho2)
492 rhs(n) = rhs(n) - rterm + drterm * h + rho2 * bt
508 subroutine sto_cq(this, flowja, hnew, hold)
513 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
514 real(DP),
dimension(:),
contiguous,
intent(in) :: hnew
515 real(DP),
dimension(:),
contiguous,
intent(in) :: hold
518 integer(I4B) :: idiag
537 do n = 1, this%dis%nodes
538 this%strgss(n) =
dzero
539 this%strgsy(n) =
dzero
543 if (this%iss == 0)
then
549 do n = 1, this%dis%nodes
550 if (this%ibound(n) <= 0) cycle
556 if (this%iconvert(n) == 0)
then
565 sc1 =
sscapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
568 if (this%integratechanges /= 0)
then
572 sc1old =
sscapacity(this%istor_coef, tp, bt, this%dis%area(n), &
574 rho1old = sc1old * tled
582 call ssterms(this%iconvert(n), this%iorig_ss, this%iconf_ss, tp, bt, &
583 rho1, rho1old, snnew, snold, hnew(n), hold(n), &
584 aterm, rhsterm, rate)
587 this%strgss(n) = rate
590 idiag = this%dis%con%ia(n)
591 flowja(idiag) = flowja(idiag) + rate
595 if (this%iconvert(n) /= 0)
then
598 sc2 =
sycapacity(this%dis%area(n), this%sy(n))
601 if (this%integratechanges /= 0)
then
605 sc2old =
sycapacity(this%dis%area(n), this%oldsy(n))
606 rho2old = sc2old * tled
614 call syterms(tp, bt, rho2, rho2old, snnew, snold, &
615 aterm, rhsterm, rate)
618 this%strgsy(n) = rate
621 idiag = this%dis%con%ia(n)
622 flowja(idiag) = flowja(idiag) + rate
636 subroutine sto_bd(this, isuppress_output, model_budget)
642 integer(I4B),
intent(in) :: isuppress_output
643 type(
budgettype),
intent(inout) :: model_budget
650 call model_budget%addentry(rin, rout,
delt,
budtxt(1), &
651 isuppress_output,
' STORAGE')
654 if (this%iusesy == 1)
then
656 call model_budget%addentry(rin, rout,
delt,
budtxt(2), &
657 isuppress_output,
' STORAGE')
672 integer(I4B),
intent(in) :: icbcfl
673 integer(I4B),
intent(in) :: icbcun
675 integer(I4B) :: ibinun
676 integer(I4B) :: iprint, nvaluesp, nwidthp
677 character(len=1) :: cdatafmp =
' ', editdesc =
' '
681 if (this%ipakcb < 0)
then
683 elseif (this%ipakcb == 0)
then
688 if (icbcfl == 0) ibinun = 0
691 if (ibinun /= 0)
then
696 call this%dis%record_array(this%strgss, this%iout, iprint, -ibinun, &
697 budtxt(1), cdatafmp, nvaluesp, &
698 nwidthp, editdesc, dinact)
701 if (this%iusesy == 1)
then
702 call this%dis%record_array(this%strgsy, this%iout, iprint, -ibinun, &
703 budtxt(2), cdatafmp, nvaluesp, &
704 nwidthp, editdesc, dinact)
724 if (this%intvs /= 0)
then
726 deallocate (this%tvs)
730 if (this%inunit > 0)
then
738 if (
associated(this%oldss))
then
741 if (
associated(this%oldsy))
then
756 call this%NumericalPackageType%da()
775 call this%NumericalPackageType%allocate_scalars()
778 call mem_allocate(this%istor_coef,
'ISTOR_COEF', this%memoryPath)
779 call mem_allocate(this%iconf_ss,
'ICONF_SS', this%memoryPath)
780 call mem_allocate(this%iorig_ss,
'IORIG_SS', this%memoryPath)
781 call mem_allocate(this%iusesy,
'IUSESY', this%memoryPath)
782 call mem_allocate(this%satomega,
'SATOMEGA', this%memoryPath)
783 call mem_allocate(this%integratechanges,
'INTEGRATECHANGES', &
792 this%satomega =
dzero
793 this%integratechanges = 0
810 integer(I4B),
intent(in) :: nodes
815 call mem_allocate(this%iconvert, nodes,
'ICONVERT', this%memoryPath)
816 call mem_allocate(this%ss, nodes,
'SS', this%memoryPath)
817 call mem_allocate(this%sy, nodes,
'SY', this%memoryPath)
818 call mem_allocate(this%strgss, nodes,
'STRGSS', this%memoryPath)
819 call mem_allocate(this%strgsy, nodes,
'STRGSY', this%memoryPath)
827 this%strgss(n) =
dzero
828 this%strgsy(n) =
dzero
829 if (this%integratechanges /= 0)
then
830 this%oldss(n) =
dzero
831 if (this%iusesy /= 0)
then
832 this%oldsy(n) =
dzero
851 character(len=LINELENGTH) :: keyword, fname
853 logical :: isfound, endOfBlock
855 character(len=*),
parameter :: fmtisvflow = &
856 "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
857 &WHENEVER ICBCFL IS NOT ZERO.')"
858 character(len=*),
parameter :: fmtflow = &
859 &
"(4x, 'FLOWS WILL BE SAVED TO FILE: ', a, /4x, 'OPENED ON UNIT: ', I7)"
860 character(len=*),
parameter :: fmtorigss = &
861 "(4X,'ORIGINAL_SPECIFIC_STORAGE OPTION:',/, &
862 &1X,'The original specific storage formulation will be used')"
863 character(len=*),
parameter :: fmtstoc = &
864 "(4X,'STORAGECOEFFICIENT OPTION:',/, &
865 &1X,'Read storage coefficient rather than specific storage')"
866 character(len=*),
parameter :: fmtconfss = &
867 "(4X,'SS_CONFINED_ONLY OPTION:',/, &
868 &1X,'Specific storage changes only occur under confined conditions')"
871 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, &
872 supportopenclose=.true., blockrequired=.false.)
876 write (this%iout,
'(1x,a)')
'PROCESSING STORAGE OPTIONS'
878 call this%parser%GetNextLine(endofblock)
880 call this%parser%GetStringCaps(keyword)
881 select case (keyword)
884 write (this%iout, fmtisvflow)
885 case (
'STORAGECOEFFICIENT')
887 write (this%iout, fmtstoc)
888 case (
'SS_CONFINED_ONLY')
891 write (this%iout, fmtconfss)
893 if (this%intvs /= 0)
then
894 errmsg =
'Multiple TVS6 keywords detected in OPTIONS block.'// &
895 ' Only one TVS6 entry allowed.'
898 call this%parser%GetStringCaps(keyword)
899 if (trim(adjustl(keyword)) /=
'FILEIN')
then
900 errmsg =
'TVS6 keyword must be followed by "FILEIN" '// &
904 call this%parser%GetString(fname)
906 call openfile(this%intvs, this%iout, fname,
'TVS')
907 call tvs_cr(this%tvs, this%name_model, this%intvs, this%iout)
913 case (
'DEV_ORIGINAL_SPECIFIC_STORAGE')
915 write (this%iout, fmtorigss)
916 case (
'DEV_OLDSTORAGEFORMULATION')
917 call this%parser%DevOpt()
920 write (this%iout, fmtconfss)
922 write (
errmsg,
'(a,a)')
'Unknown STO option: ', &
927 write (this%iout,
'(1x,a)')
'END OF STORAGE OPTIONS'
931 if (this%inewton > 0)
then
949 character(len=LINELENGTH) :: keyword
950 character(len=:),
allocatable :: line
951 character(len=LINELENGTH) :: cellstr
952 integer(I4B) :: istart, istop, lloc, ierr
953 logical :: isfound, endOfBlock
958 character(len=24),
dimension(4) :: aname
962 data aname(1)/
' ICONVERT'/
963 data aname(2)/
' SPECIFIC STORAGE'/
964 data aname(3)/
' SPECIFIC YIELD'/
965 data aname(4)/
' STORAGE COEFFICIENT'/
975 call this%parser%GetBlock(
'GRIDDATA', isfound, ierr)
977 write (this%iout,
'(1x,a)')
'PROCESSING GRIDDATA'
979 call this%parser%GetNextLine(endofblock)
981 call this%parser%GetStringCaps(keyword)
982 call this%parser%GetRemainingLine(line)
984 select case (keyword)
986 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
987 this%parser%iuactive, this%iconvert, &
991 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
992 this%parser%iuactive, this%ss, &
996 call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
997 this%parser%iuactive, this%sy, &
1001 write (
errmsg,
'(a,a)')
'Unknown GRIDDATA tag: ', &
1004 call this%parser%StoreErrorUnit()
1007 write (this%iout,
'(1x,a)')
'END PROCESSING GRIDDATA'
1009 write (
errmsg,
'(a)')
'Required GRIDDATA block not found.'
1011 call this%parser%StoreErrorUnit()
1015 if (.not. readiconv)
then
1016 write (
errmsg,
'(a, a, a)')
'Error in GRIDDATA block: ', &
1017 trim(adjustl(aname(1))),
' not found.'
1021 do n = 1, this%dis%nodes
1022 if (this%iconvert(n) /= 0)
then
1031 if (.not. readss)
then
1032 write (
errmsg,
'(a, a, a)')
'Error in GRIDDATA block: ', &
1033 trim(adjustl(aname(2))),
' not found.'
1038 if (.not. readsy .and. isconv)
then
1039 write (
errmsg,
'(a, a, a)')
'Error in GRIDDATA block: ', &
1040 trim(adjustl(aname(3))),
' not found.'
1045 call this%parser%StoreErrorUnit()
1049 do n = 1, this%dis%nodes
1050 if (this%ss(n) <
dzero)
then
1051 call this%dis%noder_to_string(n, cellstr)
1052 write (
errmsg,
'(a,2(1x,a),1x,g0,1x,a)') &
1053 'Error in SS DATA: SS value in cell', trim(adjustl(cellstr)), &
1054 'is less than zero (', this%ss(n),
').'
1058 if (this%sy(n) <
dzero)
then
1059 call this%dis%noder_to_string(n, cellstr)
1060 write (
errmsg,
'(a,2(1x,a),1x,g0,1x,a)') &
1061 'Error in SY DATA: SY value in cell', trim(adjustl(cellstr)), &
1062 'is less than zero (', this%sy(n),
').'
1087 if (.not.
associated(this%oldss))
then
1088 call mem_allocate(this%oldss, this%dis%nodes,
'OLDSS', this%memoryPath)
1090 if (this%iusesy == 1 .and. .not.
associated(this%oldsy))
then
1091 call mem_allocate(this%oldsy, this%dis%nodes,
'OLDSY', this%memoryPath)
1095 do n = 1, this%dis%nodes
1096 this%oldss(n) = this%ss(n)
1100 if (this%iusesy == 1)
then
1101 do n = 1, this%dis%nodes
1102 this%oldsy(n) = this%sy(n)
This module contains block parser methods.
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 dem4
real constant 1e-4
real(dp), parameter dem6
real constant 1e-6
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
This module contains the storage package methods.
subroutine sto_bd(this, isuppress_output, model_budget)
@ brief Model budget calculation for package
subroutine read_data(this)
@ brief Read data for package
subroutine, public sto_cr(stoobj, name_model, inunit, iout)
@ brief Create a new package object
subroutine sto_ar(this, dis, ibound)
@ brief Allocate and read method for package
subroutine sto_da(this)
@ brief Deallocate package memory
subroutine sto_fc(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
@ brief Fill A and right-hand side for the package
subroutine read_options(this)
@ brief Read options for package
subroutine sto_fn(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
@ brief Fill Newton-Raphson terms in A and right-hand side for the package
subroutine sto_rp(this)
@ brief Read and prepare method for package
subroutine sto_cq(this, flowja, hnew, hold)
@ brief Calculate flows for package
subroutine save_old_ss_sy(this)
@ brief Save old storage property values
subroutine sto_save_model_flows(this, icbcfl, icbcun)
@ brief Save model flows for package
subroutine allocate_scalars(this)
@ brief Allocate scalars
subroutine allocate_arrays(this, nodes)
@ brief Allocate package arrays
subroutine sto_ad(this)
@ brief Advance the package
character(len=lenbudtxt), dimension(2) budtxt
This module contains stateless storage subroutines and functions.
pure real(dp) function, public sycapacity(area, sy)
Calculate the specific yield capacity.
pure subroutine, public syterms(top, bot, rho2, rho2old, snnew, snold, aterm, rhsterm, rate)
Calculate the specific yield storage terms.
pure subroutine, public ssterms(iconvert, iorig_ss, iconf_ss, top, bot, rho1, rho1old, snnew, snold, hnew, hold, aterm, rhsterm, rate)
Calculate the specific storage terms.
pure real(dp) function, public sscapacity(istor_coef, top, bot, area, ss)
Calculate the specific storage capacity.
This module defines variable data types.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
This module contains the base numerical package type.
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.
character(len=maxcharlen) errmsg
error message string
real(dp) function squadraticsaturation(top, bot, x, eps)
@ brief sQuadraticSaturation
real(dp) function slinearsaturation(top, bot, x)
@ brief sLinearSaturation
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function
real(dp) function sqsaturation(top, bot, x, c1, c2)
@ brief sQSaturation
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number
real(dp), pointer, public delt
length of the current time step
integer(i4b), pointer, public nper
number of stress period
This module contains the time-varying storage package methods.
subroutine, public tvs_cr(tvs, name_model, inunit, iout)
Create a new TvsType object.
Derived type for the Budget object.