26 character(len=LENBUDTXT),
dimension(1) ::
budtxt = & !< text labels for budget terms
30 integer(I4B),
pointer :: iss => null()
31 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
32 real(dp),
dimension(:),
pointer,
contiguous :: qsto => null()
35 integer(I4B),
dimension(:),
pointer,
contiguous :: idcxs => null()
69 subroutine sto_cr(stoobj, name_model, inunit, iout, cxs)
72 character(len=*),
intent(in) :: name_model
73 integer(I4B),
intent(in) :: inunit
74 integer(I4B),
intent(in) :: iout
81 call stoobj%set_names(1, name_model,
'STO',
'STO')
84 call stoobj%allocate_scalars()
87 stoobj%inunit = inunit
94 call stoobj%parser%Initialize(stoobj%inunit, stoobj%iout)
110 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
113 character(len=*),
parameter :: fmtsto = &
114 "(1x,/1x,'STO -- STORAGE PACKAGE, VERSION 1, 10/27/2023', &
115 &' INPUT READ FROM UNIT ', i0, //)"
118 write (this%iout, fmtsto) this%inunit
121 call this%set_dfw_pointers()
126 this%ibound => ibound
132 call this%allocate_arrays(dis%nodes)
135 call this%read_options()
158 logical :: isfound, endOfBlock
159 character(len=16) :: css(0:1)
160 character(len=LINELENGTH) :: line, keyword
162 character(len=*),
parameter :: fmtlsp = &
163 &
"(1X,/1X,'REUSING ',A,' FROM LAST STRESS PERIOD')"
164 character(len=*),
parameter :: fmtblkerr = &
165 &
"('Error. Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
167 data css(0)/
' TRANSIENT'/
168 data css(1)/
' STEADY-STATE'/
171 if (this%ionper <
kper)
then
174 call this%parser%GetBlock(
'PERIOD', isfound, ierr, &
175 supportopenclose=.true., &
176 blockrequired=.false.)
180 call this%read_check_ionper()
186 this%ionper =
nper + 1
189 call this%parser%GetCurrentLine(line)
190 write (
errmsg, fmtblkerr) adjustl(trim(line))
192 call this%parser%StoreErrorUnit()
197 if (this%ionper ==
kper)
then
198 write (this%iout,
'(//,1x,a)')
'PROCESSING STORAGE PERIOD DATA'
200 call this%parser%GetNextLine(endofblock)
202 call this%parser%GetStringCaps(keyword)
203 select case (keyword)
204 case (
'STEADY-STATE')
209 write (
errmsg,
'(a,a)')
'Unknown STORAGE data tag: ', &
212 call this%parser%StoreErrorUnit()
215 write (this%iout,
'(1x,a)')
'END PROCESSING STORAGE PERIOD DATA'
218 write (this%iout,
'(//1X,A,I0,A,A,/)') &
219 'STRESS PERIOD ',
kper,
' IS ', trim(adjustl(css(this%iss)))
244 subroutine sto_fc(this, kiter, stage_old, stage_new, matrix_sln, idxglo, rhs)
249 integer(I4B) :: kiter
250 real(DP),
intent(inout),
dimension(:) :: stage_old
251 real(DP),
intent(inout),
dimension(:) :: stage_new
253 integer(I4B),
intent(in),
dimension(:) :: idxglo
254 real(DP),
intent(inout),
dimension(:) :: rhs
256 character(len=LINELENGTH) :: distype =
''
258 character(len=*),
parameter :: fmtsperror = &
259 &
"('Detected time step length of zero. SWF Storage Package cannot be ', &
260 &'used unless delt is non-zero.')"
263 if (this%iss /= 0)
return
267 write (
errmsg, fmtsperror)
271 call this%dis%get_dis_type(distype)
272 if (distype ==
'DISV1D')
then
273 call this%sto_fc_dis1d(kiter, stage_old, stage_new, matrix_sln, idxglo, rhs)
275 call this%sto_fc_dis2d(kiter, stage_old, stage_new, matrix_sln, idxglo, rhs)
285 subroutine sto_fc_dis1d(this, kiter, stage_old, stage_new, matrix_sln, &
290 integer(I4B) :: kiter
291 real(DP),
intent(inout),
dimension(:) :: stage_old
292 real(DP),
intent(inout),
dimension(:) :: stage_new
294 integer(I4B),
intent(in),
dimension(:) :: idxglo
295 real(DP),
intent(inout),
dimension(:) :: rhs
297 integer(I4B) :: n, idiag
300 real(DP),
dimension(:),
pointer :: reach_length
303 reach_length => this%reach_length_pointer()
306 do n = 1, this%dis%nodes
309 if (this%ibound(n) < 0) cycle
311 call this%calc_storage_dis1d(n, stage_new(n), stage_old(n), &
312 reach_length(n), qsto, derv)
315 idiag = this%dis%con%ia(n)
316 call matrix_sln%add_value_pos(idxglo(idiag), -derv)
317 rhs(n) = rhs(n) + qsto - derv * stage_new(n)
330 subroutine sto_fc_dis2d(this, kiter, stage_old, stage_new, matrix_sln, &
335 integer(I4B) :: kiter
336 real(DP),
intent(inout),
dimension(:) :: stage_old
337 real(DP),
intent(inout),
dimension(:) :: stage_new
339 integer(I4B),
intent(in),
dimension(:) :: idxglo
340 real(DP),
intent(inout),
dimension(:) :: rhs
342 integer(I4B) :: n, idiag
347 do n = 1, this%dis%nodes
350 if (this%ibound(n) < 0) cycle
353 call this%calc_storage_dis2d(n, stage_new(n), stage_old(n), &
357 idiag = this%dis%con%ia(n)
358 call matrix_sln%add_value_pos(idxglo(idiag), -derv)
359 rhs(n) = rhs(n) + qsto - derv * stage_new(n)
367 subroutine sto_cq(this, flowja, stage_new, stage_old)
370 real(DP),
intent(inout),
dimension(:) :: flowja
371 real(DP),
intent(inout),
dimension(:) :: stage_new
372 real(DP),
intent(inout),
dimension(:) :: stage_old
374 real(DP),
dimension(:),
pointer :: reach_length
376 integer(I4B) :: idiag
381 if (this%iss /= 0)
return
384 reach_length => this%reach_length_pointer()
387 do n = 1, this%dis%nodes
390 if (this%ibound(n) < 0) cycle
394 if (
associated(reach_length))
then
396 call this%calc_storage_dis1d(n, stage_new(n), stage_old(n), dx, q)
398 call this%calc_storage_dis2d(n, stage_new(n), stage_old(n), q)
401 idiag = this%dis%con%ia(n)
402 flowja(idiag) = flowja(idiag) + this%qsto(n)
416 integer(I4B),
intent(in) :: n
417 real(DP),
intent(in) :: stage_new
418 real(DP),
intent(in) :: stage_old
419 real(DP),
intent(in) :: dx
420 real(DP),
intent(inout) :: qsto
421 real(DP),
intent(inout),
optional :: derv
423 real(DP) :: depth_new
424 real(DP) :: depth_old
427 real(DP) :: cxs_area_new
428 real(DP) :: cxs_area_old
429 real(DP) :: cxs_area_eps
432 call this%dis%get_flow_width(n, n, 0, width_n, width_m)
433 depth_new = stage_new - this%dis%bot(n)
434 depth_old = stage_old - this%dis%bot(n)
435 cxs_area_new = this%cxs%get_area(this%idcxs(n), width_n, depth_new)
436 cxs_area_old = this%cxs%get_area(this%idcxs(n), width_n, depth_old)
437 qsto = (cxs_area_new - cxs_area_old) * dx /
delt
438 if (
present(derv))
then
440 cxs_area_eps = this%cxs%get_area(this%idcxs(n), width_n, depth_new + eps)
441 derv = (cxs_area_eps - cxs_area_new) * dx /
delt / eps
452 integer(I4B),
intent(in) :: n
453 real(DP),
intent(in) :: stage_new
454 real(DP),
intent(in) :: stage_old
455 real(DP),
intent(inout) :: qsto
456 real(DP),
intent(inout),
optional :: derv
459 real(DP) :: depth_new
460 real(DP) :: depth_old
461 real(DP) :: depth_eps
462 real(DP) :: volume_new
463 real(DP) :: volume_old
466 area = this%dis%get_area(n)
467 depth_new = stage_new - this%dis%bot(n)
468 depth_old = stage_old - this%dis%bot(n)
469 volume_new = area * depth_new
470 volume_old = area * depth_old
471 qsto = (volume_new - volume_old) /
delt
473 if (
present(derv))
then
475 depth_eps = depth_new + eps
476 derv = (depth_eps - depth_new) * area /
delt / eps
487 subroutine sto_bd(this, isuppress_output, model_budget)
493 integer(I4B),
intent(in) :: isuppress_output
494 type(
budgettype),
intent(inout) :: model_budget
501 call model_budget%addentry(rin, rout,
delt,
' STO', &
502 isuppress_output,
' STORAGE')
516 integer(I4B),
intent(in) :: icbcfl
517 integer(I4B),
intent(in) :: icbcun
519 integer(I4B) :: ibinun
520 integer(I4B) :: iprint, nvaluesp, nwidthp
521 character(len=1) :: cdatafmp =
' ', editdesc =
' '
525 if (this%ipakcb < 0)
then
527 elseif (this%ipakcb == 0)
then
532 if (icbcfl == 0) ibinun = 0
535 if (ibinun /= 0)
then
540 call this%dis%record_array(this%qsto, this%iout, iprint, -ibinun, &
541 budtxt(1), cdatafmp, nvaluesp, &
542 nwidthp, editdesc, dinact)
561 if (this%inunit > 0)
then
569 call this%NumericalPackageType%da()
588 call this%NumericalPackageType%allocate_scalars()
610 integer(I4B),
intent(in) :: nodes
615 call mem_allocate(this%qsto, nodes,
'STRGSS', this%memoryPath)
637 character(len=LINELENGTH) :: keyword
639 logical :: isfound, endOfBlock
641 character(len=*),
parameter :: fmtisvflow = &
642 "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
643 &WHENEVER ICBCFL IS NOT ZERO.')"
644 character(len=*),
parameter :: fmtflow = &
645 &
"(4x, 'FLOWS WILL BE SAVED TO FILE: ', a, /4x, 'OPENED ON UNIT: ', I7)"
648 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, &
649 supportopenclose=.true., blockrequired=.false.)
653 write (this%iout,
'(1x,a)')
'PROCESSING STORAGE OPTIONS'
655 call this%parser%GetNextLine(endofblock)
657 call this%parser%GetStringCaps(keyword)
658 select case (keyword)
661 write (this%iout, fmtisvflow)
663 write (
errmsg,
'(a,a)')
'Unknown STO option: ', &
668 write (this%iout,
'(1x,a)')
'END OF STORAGE OPTIONS'
685 character(len=LINELENGTH) :: keyword
686 character(len=:),
allocatable :: line
687 integer(I4B) :: lloc, ierr
688 logical :: isfound, endOfBlock
689 character(len=24),
dimension(1) :: aname
692 data aname(1)/
' XXX'/
698 call this%parser%GetBlock(
'GRIDDATA', isfound, ierr)
700 write (this%iout,
'(1x,a)')
'PROCESSING GRIDDATA'
702 call this%parser%GetNextLine(endofblock)
704 call this%parser%GetStringCaps(keyword)
705 call this%parser%GetRemainingLine(line)
707 select case (keyword)
709 write (
errmsg,
'(a,a)')
'Unknown GRIDDATA tag: ', &
712 call this%parser%StoreErrorUnit()
715 write (this%iout,
'(1x,a)')
'END PROCESSING GRIDDATA'
717 write (
errmsg,
'(a)')
'Required GRIDDATA block not found.'
719 call this%parser%StoreErrorUnit()
723 call this%parser%StoreErrorUnit()
738 character(len=LENMEMPATH) :: dfw_mem_path
741 call mem_setptr(this%idcxs,
'IDCXS', dfw_mem_path)
749 real(dp),
dimension(:),
pointer :: ptr
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
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter done
real constant 1
This module defines variable data types.
real(dp) function, public get_perturbation(x)
Calculate a numerical perturbation given the value of x.
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
This module contains the storage package methods.
subroutine sto_fc_dis2d(this, kiter, stage_old, stage_new, matrix_sln, idxglo, rhs)
@ brief Fill A and right-hand side for the package
subroutine calc_storage_dis2d(this, n, stage_new, stage_old, qsto, derv)
subroutine sto_save_model_flows(this, icbcfl, icbcun)
@ brief Save model flows for package
subroutine read_options(this)
@ brief Read options for package
subroutine sto_fc(this, kiter, stage_old, stage_new, matrix_sln, idxglo, rhs)
@ brief Fill A and right-hand side for the package
subroutine sto_bd(this, isuppress_output, model_budget)
@ brief Model budget calculation for package
subroutine sto_cq(this, flowja, stage_new, stage_old)
@ brief Calculate flows for package
subroutine sto_rp(this)
@ brief Read and prepare method for package
subroutine sto_ad(this)
@ brief Advance the package
subroutine sto_ar(this, dis, ibound)
@ brief Allocate and read method for package
subroutine read_data(this)
@ brief Read data for package
subroutine sto_da(this)
@ brief Deallocate package memory
subroutine sto_fc_dis1d(this, kiter, stage_old, stage_new, matrix_sln, idxglo, rhs)
@ brief Fill A and right-hand side for the package
real(dp) function, dimension(:), pointer reach_length_pointer(this)
subroutine allocate_arrays(this, nodes)
@ brief Allocate package arrays
subroutine calc_storage_dis1d(this, n, stage_new, stage_old, dx, qsto, derv)
subroutine, public sto_cr(stoobj, name_model, inunit, iout, cxs)
@ brief Create a new package object
subroutine set_dfw_pointers(this)
Set pointers to channel properties in DFW Package.
subroutine allocate_scalars(this)
@ brief Allocate scalars
character(len=lenbudtxt), dimension(1) budtxt
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
Derived type for the Budget object.