46 character(len=LENFTYPE) ::
ftype =
'SFR'
47 character(len=LENPACKAGENAME) ::
text =
' SFR'
57 character(len=16),
dimension(:),
pointer,
contiguous :: csfrbudget => null()
58 character(len=16),
dimension(:),
pointer,
contiguous :: cauxcbc => null()
59 character(len=LENBOUNDNAME),
dimension(:),
pointer, &
60 contiguous :: sfrname => null()
62 integer(I4B),
pointer :: iprhed => null()
63 integer(I4B),
pointer :: istageout => null()
64 integer(I4B),
pointer :: ibudgetout => null()
65 integer(I4B),
pointer :: ibudcsv => null()
66 integer(I4B),
pointer :: ipakcsv => null()
67 integer(I4B),
pointer :: idiversions => null()
68 integer(I4B),
pointer :: nconn => null()
69 integer(I4B),
pointer :: maxsfrpicard => null()
70 integer(I4B),
pointer :: maxsfrit => null()
71 integer(I4B),
pointer :: bditems => null()
72 integer(I4B),
pointer :: cbcauxitems => null()
73 integer(I4B),
pointer :: icheck => null()
74 integer(I4B),
pointer :: iconvchk => null()
75 integer(I4B),
pointer :: gwfiss => null()
76 integer(I4B),
pointer :: ianynone => null()
78 real(dp),
pointer :: unitconv => null()
79 real(dp),
pointer :: lengthconv => null()
80 real(dp),
pointer :: timeconv => null()
81 real(dp),
pointer :: dmaxchg => null()
82 real(dp),
pointer :: deps => null()
84 integer(I4B),
dimension(:),
pointer,
contiguous :: isfrorder => null()
85 integer(I4B),
dimension(:),
pointer,
contiguous :: ia => null()
86 integer(I4B),
dimension(:),
pointer,
contiguous :: ja => null()
88 real(dp),
dimension(:),
pointer,
contiguous :: qoutflow => null()
89 real(dp),
dimension(:),
pointer,
contiguous :: qextoutflow => null()
90 real(dp),
dimension(:),
pointer,
contiguous :: qauxcbc => null()
91 real(dp),
dimension(:),
pointer,
contiguous :: dbuff => null()
101 integer(I4B),
dimension(:),
pointer,
contiguous :: iboundpak => null()
102 integer(I4B),
dimension(:),
pointer,
contiguous :: igwfnode => null()
103 integer(I4B),
dimension(:),
pointer,
contiguous :: igwftopnode => null()
104 real(dp),
dimension(:),
pointer,
contiguous :: length => null()
105 real(dp),
dimension(:),
pointer,
contiguous :: width => null()
106 real(dp),
dimension(:),
pointer,
contiguous :: strtop => null()
107 real(dp),
dimension(:),
pointer,
contiguous :: bthick => null()
108 real(dp),
dimension(:),
pointer,
contiguous :: hk => null()
109 real(dp),
dimension(:),
pointer,
contiguous :: slope => null()
110 integer(I4B),
dimension(:),
pointer,
contiguous :: nconnreach => null()
111 real(dp),
dimension(:),
pointer,
contiguous :: ustrf => null()
112 real(dp),
dimension(:),
pointer,
contiguous :: ftotnd => null()
113 integer(I4B),
dimension(:),
pointer,
contiguous :: ndiv => null()
114 real(dp),
dimension(:),
pointer,
contiguous :: usflow => null()
115 real(dp),
dimension(:),
pointer,
contiguous :: dsflow => null()
116 real(dp),
dimension(:),
pointer,
contiguous :: depth => null()
117 real(dp),
dimension(:),
pointer,
contiguous :: stage => null()
118 real(dp),
dimension(:),
pointer,
contiguous :: gwflow => null()
119 real(dp),
dimension(:),
pointer,
contiguous :: simevap => null()
120 real(dp),
dimension(:),
pointer,
contiguous :: simrunoff => null()
121 real(dp),
dimension(:),
pointer,
contiguous :: stage0 => null()
122 real(dp),
dimension(:),
pointer,
contiguous :: usflow0 => null()
124 integer(I4B),
pointer :: ncrossptstot => null()
125 integer(I4B),
dimension(:),
pointer,
contiguous :: ncrosspts => null()
126 integer(I4B),
dimension(:),
pointer,
contiguous :: iacross => null()
127 real(dp),
dimension(:),
pointer,
contiguous :: station => null()
128 real(dp),
dimension(:),
pointer,
contiguous :: xsheight => null()
129 real(dp),
dimension(:),
pointer,
contiguous :: xsrough => null()
131 integer(I4B),
dimension(:),
pointer,
contiguous :: idir => null()
132 integer(I4B),
dimension(:),
pointer,
contiguous :: idiv => null()
133 real(dp),
dimension(:),
pointer,
contiguous :: qconn => null()
135 real(dp),
dimension(:),
pointer,
contiguous :: rough => null()
136 real(dp),
dimension(:),
pointer,
contiguous :: rain => null()
137 real(dp),
dimension(:),
pointer,
contiguous :: evap => null()
138 real(dp),
dimension(:),
pointer,
contiguous :: inflow => null()
139 real(dp),
dimension(:),
pointer,
contiguous :: runoff => null()
140 real(dp),
dimension(:),
pointer,
contiguous :: sstage => null()
142 real(dp),
dimension(:, :),
pointer,
contiguous :: rauxvar => null()
144 integer(I4B),
dimension(:),
pointer,
contiguous :: iadiv => null()
145 integer(I4B),
dimension(:),
pointer,
contiguous :: divreach => null()
146 character(len=10),
dimension(:),
pointer,
contiguous :: divcprior => null()
147 real(dp),
dimension(:),
pointer,
contiguous :: divflow => null()
148 real(dp),
dimension(:),
pointer,
contiguous :: divq => null()
151 integer(I4B),
pointer :: idense
152 real(dp),
dimension(:, :),
pointer,
contiguous :: denseterms => null()
155 real(dp),
dimension(:, :),
pointer,
contiguous :: viscratios => null()
235 subroutine sfr_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
239 class(
bndtype),
pointer :: packobj
240 integer(I4B),
intent(in) :: id
241 integer(I4B),
intent(in) :: ibcnum
242 integer(I4B),
intent(in) :: inunit
243 integer(I4B),
intent(in) :: iout
244 character(len=*),
intent(in) :: namemodel
245 character(len=*),
intent(in) :: pakname
247 type(
sfrtype),
pointer :: sfrobj
254 call packobj%set_names(ibcnum, namemodel, pakname,
ftype)
258 call sfrobj%sfr_allocate_scalars()
261 call packobj%pack_initialize()
263 packobj%inunit = inunit
266 packobj%ibcnum = ibcnum
287 class(
sfrtype),
intent(inout) :: this
290 call this%BndType%allocate_scalars()
293 call mem_allocate(this%iprhed,
'IPRHED', this%memoryPath)
294 call mem_allocate(this%istageout,
'ISTAGEOUT', this%memoryPath)
295 call mem_allocate(this%ibudgetout,
'IBUDGETOUT', this%memoryPath)
296 call mem_allocate(this%ibudcsv,
'IBUDCSV', this%memoryPath)
297 call mem_allocate(this%ipakcsv,
'IPAKCSV', this%memoryPath)
298 call mem_allocate(this%idiversions,
'IDIVERSIONS', this%memoryPath)
299 call mem_allocate(this%maxsfrpicard,
'MAXSFRPICARD', this%memoryPath)
300 call mem_allocate(this%maxsfrit,
'MAXSFRIT', this%memoryPath)
301 call mem_allocate(this%bditems,
'BDITEMS', this%memoryPath)
302 call mem_allocate(this%cbcauxitems,
'CBCAUXITEMS', this%memoryPath)
303 call mem_allocate(this%unitconv,
'UNITCONV', this%memoryPath)
304 call mem_allocate(this%lengthconv,
'LENGTHCONV', this%memoryPath)
305 call mem_allocate(this%timeconv,
'TIMECONV', this%memoryPath)
306 call mem_allocate(this%dmaxchg,
'DMAXCHG', this%memoryPath)
309 call mem_allocate(this%icheck,
'ICHECK', this%memoryPath)
310 call mem_allocate(this%iconvchk,
'ICONVCHK', this%memoryPath)
311 call mem_allocate(this%idense,
'IDENSE', this%memoryPath)
312 call mem_allocate(this%ianynone,
'IANYNONE', this%memoryPath)
313 call mem_allocate(this%ncrossptstot,
'NCROSSPTSTOT', this%memoryPath)
325 this%maxsfrpicard = 100
333 this%deps =
dp999 * this%dmaxchg
340 this%ncrossptstot = 0
355 class(
sfrtype),
intent(inout) :: this
361 allocate (this%csfrbudget(this%bditems))
363 'SFRNAME', this%memoryPath)
366 call mem_allocate(this%iboundpak, this%maxbound,
'IBOUNDPAK', &
368 call mem_allocate(this%igwfnode, this%maxbound,
'IGWFNODE', this%memoryPath)
369 call mem_allocate(this%igwftopnode, this%maxbound,
'IGWFTOPNODE', &
371 call mem_allocate(this%length, this%maxbound,
'LENGTH', this%memoryPath)
372 call mem_allocate(this%width, this%maxbound,
'WIDTH', this%memoryPath)
373 call mem_allocate(this%strtop, this%maxbound,
'STRTOP', this%memoryPath)
374 call mem_allocate(this%bthick, this%maxbound,
'BTHICK', this%memoryPath)
375 call mem_allocate(this%hk, this%maxbound,
'HK', this%memoryPath)
376 call mem_allocate(this%slope, this%maxbound,
'SLOPE', this%memoryPath)
377 call mem_allocate(this%nconnreach, this%maxbound,
'NCONNREACH', &
379 call mem_allocate(this%ustrf, this%maxbound,
'USTRF', this%memoryPath)
380 call mem_allocate(this%ftotnd, this%maxbound,
'FTOTND', this%memoryPath)
381 call mem_allocate(this%ndiv, this%maxbound,
'NDIV', this%memoryPath)
382 call mem_allocate(this%usflow, this%maxbound,
'USFLOW', this%memoryPath)
383 call mem_allocate(this%dsflow, this%maxbound,
'DSFLOW', this%memoryPath)
384 call mem_allocate(this%depth, this%maxbound,
'DEPTH', this%memoryPath)
385 call mem_allocate(this%stage, this%maxbound,
'STAGE', this%memoryPath)
386 call mem_allocate(this%gwflow, this%maxbound,
'GWFLOW', this%memoryPath)
387 call mem_allocate(this%simevap, this%maxbound,
'SIMEVAP', this%memoryPath)
388 call mem_allocate(this%simrunoff, this%maxbound,
'SIMRUNOFF', &
390 call mem_allocate(this%stage0, this%maxbound,
'STAGE0', this%memoryPath)
391 call mem_allocate(this%usflow0, this%maxbound,
'USFLOW0', this%memoryPath)
394 call mem_allocate(this%isfrorder, this%maxbound,
'ISFRORDER', &
396 call mem_allocate(this%ia, this%maxbound + 1,
'IA', this%memoryPath)
398 call mem_allocate(this%idir, 0,
'IDIR', this%memoryPath)
399 call mem_allocate(this%idiv, 0,
'IDIV', this%memoryPath)
400 call mem_allocate(this%qconn, 0,
'QCONN', this%memoryPath)
403 call mem_allocate(this%rough, this%maxbound,
'ROUGH', this%memoryPath)
404 call mem_allocate(this%rain, this%maxbound,
'RAIN', this%memoryPath)
405 call mem_allocate(this%evap, this%maxbound,
'EVAP', this%memoryPath)
406 call mem_allocate(this%inflow, this%maxbound,
'INFLOW', this%memoryPath)
407 call mem_allocate(this%runoff, this%maxbound,
'RUNOFF', this%memoryPath)
408 call mem_allocate(this%sstage, this%maxbound,
'SSTAGE', this%memoryPath)
411 call mem_allocate(this%rauxvar, this%naux, this%maxbound, &
412 'RAUXVAR', this%memoryPath)
415 call mem_allocate(this%iadiv, this%maxbound + 1,
'IADIV', this%memoryPath)
416 call mem_allocate(this%divreach, 0,
'DIVREACH', this%memoryPath)
417 call mem_allocate(this%divflow, 0,
'DIVFLOW', this%memoryPath)
418 call mem_allocate(this%divq, 0,
'DIVQ', this%memoryPath)
421 call mem_allocate(this%ncrosspts, this%maxbound,
'NCROSSPTS', &
423 call mem_allocate(this%iacross, this%maxbound + 1,
'IACROSS', &
425 call mem_allocate(this%station, this%ncrossptstot,
'STATION', &
427 call mem_allocate(this%xsheight, this%ncrossptstot,
'XSHEIGHT', &
429 call mem_allocate(this%xsrough, this%ncrossptstot,
'XSROUGH', &
434 do i = 1, this%maxbound
435 this%iboundpak(i) = 1
437 this%igwftopnode(i) = 0
438 this%length(i) =
dzero
439 this%width(i) =
dzero
440 this%strtop(i) =
dzero
441 this%bthick(i) =
dzero
443 this%slope(i) =
dzero
444 this%nconnreach(i) = 0
445 this%ustrf(i) =
dzero
446 this%ftotnd(i) =
dzero
448 this%usflow(i) =
dzero
449 this%dsflow(i) =
dzero
450 this%depth(i) =
dzero
451 this%stage(i) =
dzero
452 this%gwflow(i) =
dzero
453 this%simevap(i) =
dzero
454 this%simrunoff(i) =
dzero
455 this%stage0(i) =
dzero
456 this%usflow0(i) =
dzero
459 this%rough(i) =
dzero
462 this%inflow(i) =
dzero
463 this%runoff(i) =
dzero
464 this%sstage(i) =
dzero
468 this%rauxvar(j, i) =
dzero
472 this%ncrosspts(i) = 0
473 this%iacross(i + 1) = 0
477 do i = 1, this%ncrossptstot
478 this%station(i) =
dzero
479 this%xsheight(i) =
dzero
480 this%xsrough(i) =
dzero
484 this%csfrbudget(1) =
' RAINFALL'
485 this%csfrbudget(2) =
' EVAPORATION'
486 this%csfrbudget(3) =
' RUNOFF'
487 this%csfrbudget(4) =
' EXT-INFLOW'
488 this%csfrbudget(5) =
' GWF'
489 this%csfrbudget(6) =
' EXT-OUTFLOW'
490 this%csfrbudget(7) =
' FROM-MVR'
491 this%csfrbudget(8) =
' TO-MVR'
494 call mem_allocate(this%qoutflow, this%maxbound,
'QOUTFLOW', this%memoryPath)
495 call mem_allocate(this%qextoutflow, this%maxbound,
'QEXTOUTFLOW', &
497 do i = 1, this%maxbound
498 this%qoutflow(i) =
dzero
499 this%qextoutflow(i) =
dzero
503 if (this%istageout > 0)
then
504 call mem_allocate(this%dbuff, this%maxbound,
'DBUFF', this%memoryPath)
505 do i = 1, this%maxbound
506 this%dbuff(i) =
dzero
509 call mem_allocate(this%dbuff, 0,
'DBUFF', this%memoryPath)
513 allocate (this%cauxcbc(this%cbcauxitems))
516 call mem_allocate(this%qauxcbc, this%cbcauxitems,
'QAUXCBC', &
518 do i = 1, this%cbcauxitems
519 this%qauxcbc(i) =
dzero
523 this%cauxcbc(1) =
'FLOW-AREA '
526 call mem_allocate(this%denseterms, 3, 0,
'DENSETERMS', this%memoryPath)
529 call mem_allocate(this%viscratios, 2, 0,
'VISCRATIOS', this%memoryPath)
542 class(
sfrtype),
intent(inout) :: this
544 character(len=LINELENGTH) :: keyword
546 logical(LGP) :: isfound
547 logical(LGP) :: endOfBlock
553 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
554 supportopenclose=.true.)
558 write (this%iout,
'(/1x,a)') &
559 'PROCESSING '//trim(adjustl(this%text))//
' DIMENSIONS'
561 call this%parser%GetNextLine(endofblock)
563 call this%parser%GetStringCaps(keyword)
564 select case (keyword)
566 this%maxbound = this%parser%GetInteger()
567 write (this%iout,
'(4x,a,i0)')
'NREACHES = ', this%maxbound
570 'Unknown '//trim(this%text)//
' dimension: ', trim(keyword)
574 write (this%iout,
'(1x,a)') &
575 'END OF '//trim(adjustl(this%text))//
' DIMENSIONS'
577 call store_error(
'Required dimensions block not found.')
581 if (this%maxbound < 1)
then
583 'NREACHES was not specified or was specified incorrectly.'
589 call this%parser%StoreErrorUnit()
594 call this%define_listlabel()
597 this%ncrossptstot = this%maxbound
600 call this%sfr_allocate_arrays()
603 call this%sfr_read_packagedata()
606 call this%sfr_read_crossection()
609 call this%sfr_read_connectiondata()
612 call this%sfr_read_diversions()
615 call this%sfr_setup_budobj()
618 call this%sfr_setup_tableobj()
634 class(
sfrtype),
intent(inout) :: this
635 character(len=*),
intent(inout) :: option
636 logical(LGP),
intent(inout) :: found
639 character(len=MAXCHARLEN) :: fname
640 character(len=MAXCHARLEN) :: keyword
642 character(len=*),
parameter :: fmttimeconv = &
643 &
"(4x, 'TIME CONVERSION VALUE (',g0,') SPECIFIED.')"
644 character(len=*),
parameter :: fmtlengthconv = &
645 &
"(4x, 'LENGTH CONVERSION VALUE (',g0,') SPECIFIED.')"
646 character(len=*),
parameter :: fmtpicard = &
647 &
"(4x, 'MAXIMUM SFR PICARD ITERATION VALUE (',i0,') SPECIFIED.')"
648 character(len=*),
parameter :: fmtiter = &
649 &
"(4x, 'MAXIMUM SFR ITERATION VALUE (',i0,') SPECIFIED.')"
650 character(len=*),
parameter :: fmtdmaxchg = &
651 &
"(4x, 'MAXIMUM DEPTH CHANGE VALUE (',g0,') SPECIFIED.')"
652 character(len=*),
parameter :: fmtsfrbin = &
653 "(4x, 'SFR ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, /4x, &
654 &'OPENED ON UNIT: ', I0)"
661 write (this%iout,
'(4x,a)') trim(adjustl(this%text))// &
662 ' STAGES WILL BE PRINTED TO LISTING FILE.'
664 call this%parser%GetStringCaps(keyword)
665 if (keyword ==
'FILEOUT')
then
666 call this%parser%GetString(fname)
668 call openfile(this%istageout, this%iout, fname,
'DATA(BINARY)', &
670 write (this%iout, fmtsfrbin) &
671 'STAGE', trim(adjustl(fname)), this%istageout
674 &be followed by fileout.')
677 call this%parser%GetStringCaps(keyword)
678 if (keyword ==
'FILEOUT')
then
679 call this%parser%GetString(fname)
681 call openfile(this%ibudgetout, this%iout, fname,
'DATA(BINARY)', &
683 write (this%iout, fmtsfrbin) &
684 'BUDGET', trim(adjustl(fname)), this%ibudgetout
686 call store_error(
'Optional budget keyword must be '// &
687 'followed by fileout.')
690 call this%parser%GetStringCaps(keyword)
691 if (keyword ==
'FILEOUT')
then
692 call this%parser%GetString(fname)
694 call openfile(this%ibudcsv, this%iout, fname,
'CSV', &
695 filstat_opt=
'REPLACE')
696 write (this%iout, fmtsfrbin) &
697 'BUDGET CSV', trim(adjustl(fname)), this%ibudcsv
699 call store_error(
'OPTIONAL BUDGETCSV KEYWORD MUST BE FOLLOWED BY &
702 case (
'PACKAGE_CONVERGENCE')
703 call this%parser%GetStringCaps(keyword)
704 if (keyword ==
'FILEOUT')
then
705 call this%parser%GetString(fname)
707 call openfile(this%ipakcsv, this%iout, fname,
'CSV', &
708 filstat_opt=
'REPLACE', mode_opt=
mnormal)
709 write (this%iout, fmtsfrbin) &
710 'PACKAGE_CONVERGENCE', trim(adjustl(fname)), this%ipakcsv
712 call store_error(
'Optional package_convergence keyword must be '// &
713 'followed by fileout.')
715 case (
'UNIT_CONVERSION')
716 this%unitconv = this%parser%GetDouble()
720 'SETTING UNIT_CONVERSION DIRECTLY'
724 warnmsg, this%parser%GetUnit())
725 case (
'LENGTH_CONVERSION')
726 this%lengthconv = this%parser%GetDouble()
727 write (this%iout, fmtlengthconv) this%lengthconv
728 case (
'TIME_CONVERSION')
729 this%timeconv = this%parser%GetDouble()
730 write (this%iout, fmttimeconv) this%timeconv
731 case (
'MAXIMUM_PICARD_ITERATIONS')
732 this%maxsfrpicard = this%parser%GetInteger()
733 write (this%iout, fmtpicard) this%maxsfrpicard
734 case (
'MAXIMUM_ITERATIONS')
735 this%maxsfrit = this%parser%GetInteger()
736 write (this%iout, fmtiter) this%maxsfrit
737 case (
'MAXIMUM_DEPTH_CHANGE')
738 r = this%parser%GetDouble()
740 this%deps =
dp999 * r
741 write (this%iout, fmtdmaxchg) this%dmaxchg
744 write (this%iout,
'(4x,A)')
'MOVER OPTION ENABLED'
750 case (
'DEV_NO_CHECK')
751 call this%parser%DevOpt()
753 write (this%iout,
'(4x,A)')
'SFR CHECKS OF REACH GEOMETRY '// &
754 'RELATIVE TO MODEL GRID AND '// &
755 'REASONABLE PARAMETERS WILL NOT '// &
757 case (
'DEV_NO_FINAL_CHECK')
758 call this%parser%DevOpt()
760 write (this%iout,
'(4x,a)') &
761 'A FINAL CONVERGENCE CHECK OF THE CHANGE IN STREAM FLOW ROUTING &
762 &STAGES AND FLOWS WILL NOT BE MADE'
782 class(
sfrtype),
intent(inout) :: this
788 call this%obs%obs_ar()
791 call this%BndType%allocate_arrays()
794 if (this%inamedbound /= 0)
then
795 do n = 1, this%maxbound
796 this%boundname(n) = this%sfrname(n)
801 call this%copy_boundname()
804 do n = 1, this%maxbound
805 this%nodelist(n) = this%igwfnode(n)
809 call this%sfr_check_conversion()
812 call this%sfr_check_reaches()
815 call this%sfr_check_connections()
818 if (this%idiversions /= 0)
then
819 call this%sfr_check_diversions()
825 call this%parser%StoreErrorUnit()
829 if (this%imover /= 0)
then
830 allocate (this%pakmvrobj)
831 call this%pakmvrobj%ar(this%maxbound, this%maxbound, this%memoryPath)
847 class(
sfrtype),
intent(inout) :: this
849 character(len=LINELENGTH) :: text
850 character(len=LINELENGTH) :: cellid
851 character(len=10) :: cnum
852 character(len=LENBOUNDNAME) :: bndName
853 character(len=LENBOUNDNAME) :: bndNameTemp
854 character(len=LENBOUNDNAME) :: hkname
855 character(len=LENBOUNDNAME) :: manningname
856 character(len=LENBOUNDNAME) :: ustrfname
857 character(len=50),
dimension(:),
allocatable :: caux
858 integer(I4B) :: n, ierr, ival
859 logical(LGP) :: isfound
860 logical(LGP) :: endOfBlock
865 integer(I4B) :: nconzero
867 integer,
allocatable,
dimension(:) :: nboundchk
868 real(DP),
pointer :: bndElem => null()
871 allocate (nboundchk(this%maxbound))
872 do i = 1, this%maxbound
878 if (this%naux > 0)
then
879 allocate (caux(this%naux))
883 call this%parser%GetBlock(
'PACKAGEDATA', isfound, ierr, &
884 supportopenclose=.true.)
888 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
891 call this%parser%GetNextLine(endofblock)
894 n = this%parser%GetInteger()
896 if (n < 1 .or. n > this%maxbound)
then
897 write (
errmsg,
'(a,1x,a,1x,i0)') &
898 'Reach number (rno) must be greater than 0 and less', &
899 'than or equal to', this%maxbound
905 nboundchk(n) = nboundchk(n) + 1
908 call this%parser%GetCellid(this%dis%ndim, cellid, flag_string=.true.)
909 this%igwfnode(n) = this%dis%noder_from_cellid(cellid, this%inunit, &
911 flag_string=.true., &
913 this%igwftopnode(n) = this%igwfnode(n)
916 if (this%igwfnode(n) < 1)
then
917 this%ianynone = this%ianynone + 1
919 if (cellid ==
'NONE')
then
920 call this%parser%GetStringCaps(cellid)
923 write (cnum,
'(i0)') n
924 warnmsg =
'CELLID for unconnected reach '//trim(cnum)// &
925 ' specified to be NONE. Unconnected reaches '// &
926 'should be specified with a zero for each grid '// &
927 'dimension. For example, for a DIS grid a CELLID '// &
928 'of 0 0 0 should be specified for unconnected reaches'
932 warnmsg, this%parser%GetUnit())
938 this%length(n) = this%parser%GetDouble()
940 this%width(n) = this%parser%GetDouble()
942 this%slope(n) = this%parser%GetDouble()
944 this%strtop(n) = this%parser%GetDouble()
946 this%bthick(n) = this%parser%GetDouble()
948 call this%parser%GetStringCaps(hkname)
950 call this%parser%GetStringCaps(manningname)
952 ival = this%parser%GetInteger()
953 this%nconnreach(n) = ival
954 this%nconn = this%nconn + ival
956 write (
errmsg,
'(a,1x,i0,1x,a,i0,a)') &
957 'NCON for reach', n, &
958 'must be greater than or equal to 0 (', ival,
').'
960 else if (ival == 0)
then
961 nconzero = nconzero + 1
964 call this%parser%GetString(ustrfname)
966 ival = this%parser%GetInteger()
970 else if (ival < 0)
then
975 do iaux = 1, this%naux
976 call this%parser%GetString(caux(iaux))
980 write (cnum,
'(i10.10)') n
981 bndname =
'Reach'//cnum
984 if (this%inamedbound /= 0)
then
985 call this%parser%GetStringCaps(bndnametemp)
986 if (bndnametemp /=
'')
then
987 bndname = bndnametemp
991 this%sfrname(n) = bndname
996 bndelem => this%hk(n)
998 this%packName,
'BND', &
999 this%tsManager, this%iprpak, &
1005 bndelem => this%rough(n)
1007 this%packName,
'BND', &
1008 this%tsManager, this%iprpak, &
1014 bndelem => this%ustrf(n)
1016 this%packName,
'BND', &
1017 this%tsManager, this%iprpak,
'USTRF')
1020 do jj = 1, this%naux
1023 bndelem => this%rauxvar(jj, ii)
1025 this%packName,
'AUX', &
1026 this%tsManager, this%iprpak, &
1034 this%sstage(n) = this%strtop(n)
1037 write (this%iout,
'(1x,a)') &
1038 'END OF '//trim(adjustl(this%text))//
' PACKAGEDATA'
1040 call store_error(
'REQUIRED PACKAGEDATA BLOCK NOT FOUND.')
1045 do i = 1, this%maxbound
1046 if (nboundchk(i) == 0)
then
1047 write (
errmsg,
'(a,i0,1x,a)') &
1048 'Information for reach ', i,
'not specified in packagedata block.'
1050 else if (nboundchk(i) > 1)
then
1051 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1052 'Reach information specified', nboundchk(i),
'times for reach', i
1056 deallocate (nboundchk)
1059 if (nconzero > 0)
then
1060 write (
warnmsg,
'(a,1x,a,1x,a,1x,i0,1x, a)') &
1061 'SFR Package', trim(this%packName), &
1062 'has', nconzero,
'reach(es) with zero connections.'
1068 call this%parser%StoreErrorUnit()
1073 this%iacross(1) = ipos
1074 do i = 1, this%maxbound
1075 this%ncrosspts(i) = 1
1076 this%station(ipos) = this%width(i)
1077 this%xsheight(ipos) =
dzero
1078 this%xsrough(ipos) =
done
1080 this%iacross(i + 1) = ipos
1084 if (this%naux > 0)
then
1102 class(
sfrtype),
intent(inout) :: this
1104 character(len=LINELENGTH) :: keyword
1105 character(len=LINELENGTH) :: line
1106 logical(LGP) :: isfound
1107 logical(LGP) :: endOfBlock
1109 integer(I4B) :: ierr
1110 integer(I4B) :: ncrossptstot
1111 integer,
allocatable,
dimension(:) :: nboundchk
1115 call this%parser%GetBlock(
'CROSSSECTIONS', isfound, ierr, &
1116 supportopenclose=.true., &
1117 blockrequired=.false.)
1121 write (this%iout,
'(/1x,a)') &
1122 'PROCESSING '//trim(adjustl(this%text))//
' CROSSSECTIONS'
1125 allocate (nboundchk(this%maxbound))
1126 do n = 1, this%maxbound
1132 call cross_data%initialize(this%ncrossptstot, this%ncrosspts, &
1134 this%station, this%xsheight, &
1139 call this%parser%GetNextLine(endofblock)
1140 if (endofblock)
exit
1143 n = this%parser%GetInteger()
1146 if (n < 1 .or. n > this%maxbound)
then
1147 write (
errmsg,
'(a,1x,a,1x,i0)') &
1148 'SFR reach in crosssections block is less than one or greater', &
1155 nboundchk(n) = nboundchk(n) + 1
1158 call this%parser%GetStringCaps(keyword)
1159 select case (keyword)
1161 call this%parser%GetStringCaps(keyword)
1162 if (trim(adjustl(keyword)) /=
'FILEIN')
then
1163 errmsg =
'TAB6 keyword must be followed by "FILEIN" '// &
1168 call this%parser%GetString(line)
1169 call cross_data%read_table(n, this%width(n), &
1170 trim(adjustl(line)))
1172 write (
errmsg,
'(a,1x,i4,1x,a)') &
1173 'CROSS-SECTION TABLE ENTRY for REACH ', n, &
1174 'MUST INCLUDE TAB6 KEYWORD'
1180 write (this%iout,
'(1x,a)') &
1181 'END OF '//trim(adjustl(this%text))//
' CROSSSECTIONS'
1185 do n = 1, this%maxbound
1186 if (nboundchk(n) > 1)
then
1187 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1188 'Cross-section data for reach', n, &
1189 'specified', nboundchk(n),
'times.'
1196 call this%parser%StoreErrorUnit()
1200 ncrossptstot = cross_data%get_ncrossptstot()
1203 if (ncrossptstot /= this%ncrossptstot)
then
1204 this%ncrossptstot = ncrossptstot
1205 call mem_reallocate(this%station, this%ncrossptstot,
'STATION', &
1207 call mem_reallocate(this%xsheight, this%ncrossptstot,
'XSHEIGHT', &
1209 call mem_reallocate(this%xsrough, this%ncrossptstot,
'XSROUGH', &
1214 call cross_data%output(this%width, this%rough)
1217 call cross_data%pack(this%ncrossptstot, this%ncrosspts, &
1224 deallocate (nboundchk)
1225 call cross_data%destroy()
1226 deallocate (cross_data)
1227 nullify (cross_data)
1244 class(
sfrtype),
intent(inout) :: this
1246 character(len=LINELENGTH) :: line
1247 logical(LGP) :: isfound
1248 logical(LGP) :: endOfBlock
1253 integer(I4B) :: jcol
1254 integer(I4B) :: jcol2
1256 integer(I4B) :: ival
1257 integer(I4B) :: idir
1258 integer(I4B) :: ierr
1259 integer(I4B) :: nconnmax
1261 integer(I4B) :: ipos
1262 integer(I4B) :: istat
1263 integer(I4B),
dimension(:),
pointer,
contiguous :: rowmaxnnz => null()
1264 integer,
allocatable,
dimension(:) :: nboundchk
1265 integer,
allocatable,
dimension(:, :) :: iconndata
1267 integer(I4B),
dimension(:),
allocatable :: iup
1268 integer(I4B),
dimension(:),
allocatable :: order
1269 type(
dag) :: sfr_dag
1272 allocate (nboundchk(this%maxbound))
1273 do n = 1, this%maxbound
1280 allocate (rowmaxnnz(this%maxbound))
1281 do n = 1, this%maxbound
1282 ival = this%nconnreach(n)
1283 if (ival < 0) ival = 0
1284 rowmaxnnz(n) = ival + 1
1285 nja = nja + ival + 1
1286 if (ival > nconnmax)
then
1301 this%qconn(n) =
dzero
1305 allocate (iconndata(nconnmax, this%maxbound))
1308 do n = 1, this%maxbound
1318 call sparse%init(this%maxbound, this%maxbound, rowmaxnnz)
1321 call this%parser%GetBlock(
'CONNECTIONDATA', isfound, ierr, &
1322 supportopenclose=.true.)
1326 write (this%iout,
'(/1x,a)') &
1327 'PROCESSING '//trim(adjustl(this%text))//
' CONNECTIONDATA'
1329 call this%parser%GetNextLine(endofblock)
1330 if (endofblock)
exit
1333 n = this%parser%GetInteger()
1336 if (n < 1 .or. n > this%maxbound)
then
1337 write (
errmsg,
'(a,1x,a,1x,i0)') &
1338 'SFR reach in connectiondata block is less than one or greater', &
1345 nboundchk(n) = nboundchk(n) + 1
1348 call sparse%addconnection(n, n, 1)
1351 do i = 1, this%nconnreach(n)
1354 ival = this%parser%GetInteger()
1357 iconndata(i, n) = ival
1363 elseif (ival == 0)
then
1364 call store_error(
'Missing or zero connection reach in line:')
1369 if (ival > this%maxbound)
then
1370 call store_error(
'Reach number exceeds NREACHES in line:')
1375 call sparse%addconnection(n, ival, 1)
1379 write (this%iout,
'(1x,a)') &
1380 'END OF '//trim(adjustl(this%text))//
' CONNECTIONDATA'
1382 do n = 1, this%maxbound
1385 if (nboundchk(n) == 0)
then
1386 write (
errmsg,
'(a,1x,i0)') &
1387 'No connection data specified for reach', n
1389 else if (nboundchk(n) > 1)
then
1390 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1391 'Connection data for reach', n, &
1392 'specified', nboundchk(n),
'times.'
1397 call store_error(
'Required connectiondata block not found.')
1402 call this%parser%StoreErrorUnit()
1406 call sparse%filliaja(this%ia, this%ja, ierr, sort=.true.)
1410 write (
errmsg,
'(a,3(1x,a))') &
1411 'Could not fill', trim(this%packName), &
1412 'package IA and JA connection data.', &
1413 'Check connectivity data in connectiondata block.'
1418 do n = 1, this%maxbound
1419 do j = this%ia(n) + 1, this%ia(n + 1) - 1
1421 do jj = 1, this%nconnreach(n)
1422 jcol2 = iconndata(jj, n)
1423 if (abs(jcol2) == jcol)
then
1436 deallocate (rowmaxnnz)
1437 deallocate (nboundchk)
1438 deallocate (iconndata)
1441 call sparse%destroy()
1447 call sfr_dag%set_vertices(this%maxbound)
1450 fill_dag:
do n = 1, this%maxbound
1454 do j = this%ia(n) + 1, this%ia(n + 1) - 1
1455 if (this%idir(j) > 0)
then
1461 if (nup == 0) cycle fill_dag
1468 do j = this%ia(n) + 1, this%ia(n + 1) - 1
1469 if (this%idir(j) > 0)
then
1470 iup(ipos) = this%ja(j)
1476 call sfr_dag%set_edges(n, iup)
1483 call sfr_dag%toposort(order, istat)
1486 if (istat == -1)
then
1488 trim(adjustl(this%text))//
' PACKAGE ('// &
1489 trim(adjustl(this%packName))//
') cannot calculate a '// &
1490 'Directed Asyclic Graph for reach connectivity because '// &
1491 'of circular dependency. Using the reach number for '// &
1492 'solution ordering.'
1497 do n = 1, this%maxbound
1498 if (istat == 0)
then
1499 this%isfrorder(n) = order(n)
1501 this%isfrorder(n) = n
1506 call sfr_dag%destroy()
1507 if (istat == 0)
then
1524 class(
sfrtype),
intent(inout) :: this
1526 character(len=10) :: cnum
1527 character(len=10) :: cval
1530 integer(I4B) :: ierr
1531 integer(I4B) :: ival
1533 integer(I4B) :: ipos
1534 integer(I4B) :: jpos
1535 integer(I4B) :: ndiv
1536 integer(I4B) :: ndiversions
1537 integer(I4B) :: idivreach
1538 logical(LGP) :: isfound
1539 logical(LGP) :: endOfBlock
1540 integer(I4B) :: idiv
1541 integer,
allocatable,
dimension(:) :: iachk
1542 integer,
allocatable,
dimension(:) :: nboundchk
1548 do n = 1, this%maxbound
1549 ndiversions = ndiversions + this%ndiv(n)
1550 i0 = i0 + this%ndiv(n)
1551 this%iadiv(n + 1) = i0
1555 if (ndiversions > 0)
then
1558 allocate (this%divcprior(ndiversions))
1559 call mem_reallocate(this%divflow, ndiversions,
'DIVFLOW', this%memoryPath)
1560 call mem_reallocate(this%divq, ndiversions,
'DIVQ', this%memoryPath)
1564 do n = 1, ndiversions
1565 this%divflow(n) =
dzero
1566 this%divq(n) =
dzero
1570 call this%parser%GetBlock(
'DIVERSIONS', isfound, ierr, &
1571 supportopenclose=.true., &
1572 blockrequired=.false.)
1576 if (this%idiversions /= 0)
then
1577 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
1582 do n = 1, this%maxbound
1583 ndiv = ndiv + this%ndiv(n)
1585 allocate (iachk(this%maxbound + 1))
1586 allocate (nboundchk(ndiv))
1588 do n = 1, this%maxbound
1589 iachk(n + 1) = iachk(n) + this%ndiv(n)
1597 call this%parser%GetNextLine(endofblock)
1598 if (endofblock)
exit
1601 n = this%parser%GetInteger()
1602 if (n < 1 .or. n > this%maxbound)
then
1603 write (cnum,
'(i0)') n
1604 errmsg =
'Reach number should be between 1 and '// &
1611 if (this%ndiv(n) < 1)
then
1612 write (cnum,
'(i0)') n
1613 errmsg =
'Diversions cannot be specified '// &
1614 'for reach '//trim(cnum)
1620 ival = this%parser%GetInteger()
1621 if (ival < 1 .or. ival > this%ndiv(n))
then
1622 write (cnum,
'(i0)') n
1623 errmsg =
'Reach '//trim(cnum)
1624 write (cnum,
'(i0)') this%ndiv(n)
1625 errmsg = trim(
errmsg)//
' diversion number should be between '// &
1626 '1 and '//trim(cnum)//
'.'
1632 ipos = iachk(n) + ival - 1
1633 nboundchk(ipos) = nboundchk(ipos) + 1
1638 ival = this%parser%GetInteger()
1639 if (ival < 1 .or. ival > this%maxbound)
then
1640 write (cnum,
'(i0)') ival
1641 errmsg =
'Diversion target reach number should be '// &
1642 'between 1 and '//trim(cnum)//
'.'
1647 jpos = this%iadiv(n) + idiv - 1
1648 this%divreach(jpos) = idivreach
1651 call this%parser%GetStringCaps(cval)
1663 errmsg =
'Invalid cprior type '//trim(cval)//
'.'
1668 this%divcprior(jpos) = cval
1671 write (this%iout,
'(1x,a)')
'END OF '//trim(adjustl(this%text))// &
1674 do n = 1, this%maxbound
1675 do j = 1, this%ndiv(n)
1676 ipos = iachk(n) + j - 1
1679 if (nboundchk(ipos) == 0)
then
1680 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1681 'No data specified for reach', n,
'diversion', j
1683 else if (nboundchk(ipos) > 1)
then
1684 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
1685 'Data for reach', n,
'diversion', j, &
1686 'specified', nboundchk(ipos),
'times'
1694 deallocate (nboundchk)
1698 write (
errmsg,
'(a,1x,a)') &
1699 'A diversions block should not be', &
1700 'specified if diversions are not specified.'
1704 if (this%idiversions /= 0)
then
1705 call store_error(
'REQUIRED DIVERSIONS BLOCK NOT FOUND.')
1711 call this%parser%StoreErrorUnit()
1729 class(
sfrtype),
intent(inout) :: this
1731 character(len=LINELENGTH) :: title
1732 character(len=LINELENGTH) :: line
1733 character(len=LINELENGTH) :: crossfile
1734 integer(I4B) :: ierr
1736 integer(I4B) :: ichkustrm
1737 integer(I4B) :: ichkcross
1738 integer(I4B) :: ncrossptstot
1739 logical(LGP) :: isfound
1740 logical(LGP) :: endOfBlock
1743 character(len=*),
parameter :: fmtblkerr = &
1744 &
"('Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
1745 character(len=*),
parameter :: fmtlsp = &
1746 &
"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
1747 character(len=*),
parameter :: fmtnbd = &
1748 "(1X,/1X,'The number of active ',A,'S (',I6, &
1749 &') is greater than maximum (',I6,')')"
1759 this%nbound = this%maxbound
1763 if (this%ionper <
kper)
then
1766 call this%parser%GetBlock(
'PERIOD', isfound, ierr, &
1767 supportopenclose=.true., &
1768 blockrequired=.false.)
1772 call this%read_check_ionper()
1778 this%ionper =
nper + 1
1781 call this%parser%GetCurrentLine(line)
1782 write (
errmsg, fmtblkerr) adjustl(trim(line))
1784 call this%parser%StoreErrorUnit()
1790 if (this%ionper ==
kper)
then
1794 call cross_data%initialize(this%ncrossptstot, this%ncrosspts, &
1796 this%station, this%xsheight, &
1800 if (this%iprpak /= 0)
then
1803 title = trim(adjustl(this%text))//
' PACKAGE ('// &
1804 trim(adjustl(this%packName))//
') DATA FOR PERIOD'
1805 write (title,
'(a,1x,i6)') trim(adjustl(title)),
kper
1806 call table_cr(this%inputtab, this%packName, title)
1807 call this%inputtab%table_df(1, 4, this%iout, finalize=.false.)
1809 call this%inputtab%initialize_column(
text, 10, alignment=
tabcenter)
1811 call this%inputtab%initialize_column(
text, 20, alignment=
tableft)
1813 write (
text,
'(a,1x,i6)')
'VALUE', n
1814 call this%inputtab%initialize_column(
text, 15, alignment=
tabcenter)
1820 call this%parser%GetNextLine(endofblock)
1821 if (endofblock)
exit
1822 n = this%parser%GetInteger()
1823 if (n < 1 .or. n > this%maxbound)
then
1824 write (
errmsg,
'(a,1x,a,1x,i0,a)') &
1825 'Reach number (RNO) must be greater than 0 and', &
1826 'less than or equal to', this%maxbound,
'.'
1832 call this%sfr_set_stressperiod(n, ichkustrm, crossfile)
1835 if (this%iprpak /= 0)
then
1836 call this%parser%GetCurrentLine(line)
1837 call this%inputtab%line_to_columns(line)
1841 if (trim(adjustl(crossfile)) /=
'NONE')
then
1842 call cross_data%read_table(n, this%width(n), &
1843 trim(adjustl(crossfile)))
1848 if (this%iprpak /= 0)
then
1849 call this%inputtab%finalize_table()
1856 ncrossptstot = cross_data%get_ncrossptstot()
1859 if (ncrossptstot /= this%ncrossptstot)
then
1860 this%ncrossptstot = ncrossptstot
1861 call mem_reallocate(this%station, this%ncrossptstot,
'STATION', &
1863 call mem_reallocate(this%xsheight, this%ncrossptstot,
'XSHEIGHT', &
1865 call mem_reallocate(this%xsrough, this%ncrossptstot,
'XSROUGH', &
1870 call cross_data%output(this%width, this%rough, kstp=1,
kper=
kper)
1873 call cross_data%pack(this%ncrossptstot, this%ncrosspts, &
1880 call cross_data%destroy()
1881 deallocate (cross_data)
1882 nullify (cross_data)
1886 write (this%iout, fmtlsp) trim(this%filtyp)
1890 if (ichkustrm /= 0)
then
1891 call this%sfr_check_ustrf()
1896 call this%parser%StoreErrorUnit()
1916 integer(I4B) :: iaux
1924 call this%TsManager%ad()
1929 call this%sfr_check_ustrf()
1935 if (this%naux > 0)
then
1936 do n = 1, this%maxbound
1937 do iaux = 1, this%naux
1938 if (this%noupdateauxvar(iaux) /= 0) cycle
1939 this%auxvar(iaux, n) = this%rauxvar(iaux, n)
1945 do n = 1, this%maxbound
1946 this%usflow(n) =
dzero
1947 if (this%iboundpak(n) < 0)
then
1948 this%stage(n) = this%sstage(n)
1953 if (this%imover == 1)
then
1954 call this%pakmvrobj%ad()
1960 call this%obs%obs_ad()
1977 integer(I4B) :: igwfnode
1980 if (this%nbound == 0)
return
1983 do n = 1, this%nbound
1984 igwfnode = this%igwftopnode(n)
1985 if (igwfnode > 0)
then
1986 if (this%ibound(igwfnode) == 0)
then
1987 call this%dis%highest_active(igwfnode, this%ibound)
1990 this%igwfnode(n) = igwfnode
1991 this%nodelist(n) = igwfnode
2004 subroutine sfr_fc(this, rhs, ia, idxglo, matrix_sln)
2007 real(DP),
dimension(:),
intent(inout) :: rhs
2008 integer(I4B),
dimension(:),
intent(in) :: ia
2009 integer(I4B),
dimension(:),
intent(in) :: idxglo
2015 integer(I4B) :: ipos
2016 integer(I4B) :: node
2027 sfrpicard:
do i = 1, this%maxsfrpicard
2033 if (this%imover == 1)
then
2034 call this%pakmvrobj%fc()
2038 reachsolve:
do j = 1, this%nbound
2039 n = this%isfrorder(j)
2040 node = this%igwfnode(n)
2042 hgwf = this%xnew(node)
2049 this%stage0(n) = this%stage(n)
2050 this%usflow0(n) = this%usflow(n)
2057 if (this%iboundpak(n) /= 0)
then
2058 call this%sfr_solve(n, hgwf, hhcof, rrhs)
2060 this%depth(n) =
dzero
2061 this%stage(n) = this%strtop(n)
2063 call this%sfr_update_flows(n, v, v)
2069 this%hcof(n) = hhcof
2073 ds = s0 - this%stage(n)
2076 if (abs(ds) > abs(dsmax))
then
2083 if (abs(dsmax) <= this%dmaxchg)
then
2090 do n = 1, this%nbound
2091 node = this%nodelist(n)
2093 rhs(node) = rhs(node) + this%rhs(n)
2095 call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(n))
2108 subroutine sfr_fn(this, rhs, ia, idxglo, matrix_sln)
2111 real(DP),
dimension(:),
intent(inout) :: rhs
2112 integer(I4B),
dimension(:),
intent(in) :: ia
2113 integer(I4B),
dimension(:),
intent(in) :: idxglo
2119 integer(I4B) :: ipos
2129 do j = 1, this%nbound
2130 i = this%isfrorder(j)
2132 if (this%iboundpak(i) < 1) cycle
2134 n = this%nodelist(i)
2137 rterm = this%hcof(i) * this%xnew(n)
2139 hgwf = this%xnew(n) +
dem4
2140 call this%sfr_solve(i, hgwf, hcof1, rhs1, update=.false.)
2141 q1 = rhs1 - hcof1 * hgwf
2143 q2 = this%rhs(i) - this%hcof(i) * this%xnew(n)
2145 drterm = (q2 - q1) /
dem4
2148 call matrix_sln%add_value_pos(idxglo(ipos), drterm - this%hcof(i))
2149 rhs(n) = rhs(n) - rterm + drterm * this%xnew(n)
2162 subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
2166 class(
sfrtype),
intent(inout) :: this
2167 integer(I4B),
intent(in) :: innertot
2168 integer(I4B),
intent(in) :: kiter
2169 integer(I4B),
intent(in) :: iend
2170 integer(I4B),
intent(in) :: icnvgmod
2171 character(len=LENPAKLOC),
intent(inout) :: cpak
2172 integer(I4B),
intent(inout) :: ipak
2173 real(DP),
intent(inout) :: dpak
2175 character(len=LENPAKLOC) :: cloc
2176 character(len=LINELENGTH) :: tag
2177 integer(I4B) :: icheck
2178 integer(I4B) :: ipakfail
2179 integer(I4B) :: locdhmax
2180 integer(I4B) :: locrmax
2181 integer(I4B) :: locdqfrommvrmax
2182 integer(I4B) :: ntabrows
2183 integer(I4B) :: ntabcols
2187 real(DP) :: qtolfact
2192 real(DP) :: dqfrommvr
2193 real(DP) :: dqfrommvrmax
2196 icheck = this%iconvchk
2204 dqfrommvrmax =
dzero
2208 if (this%ipakcsv == 0)
then
2209 if (icnvgmod == 0)
then
2217 if (.not.
associated(this%pakcsvtab))
then
2222 if (this%imover == 1)
then
2223 ntabcols = ntabcols + 2
2227 call table_cr(this%pakcsvtab, this%packName,
'')
2228 call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
2229 lineseparator=.false., separator=
',', &
2233 tag =
'total_inner_iterations'
2234 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
2236 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
2238 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
2240 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
2242 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
2244 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
2246 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
2248 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
2249 tag =
'dinflowmax_loc'
2250 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
2251 if (this%imover == 1)
then
2252 tag =
'dqfrommvrmax'
2253 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
2254 tag =
'dqfrommvrmax_loc'
2255 call this%pakcsvtab%initialize_column(tag, 16, alignment=
tableft)
2261 if (icheck /= 0)
then
2262 final_check:
do n = 1, this%maxbound
2263 if (this%iboundpak(n) == 0) cycle
2266 qtolfact =
delt / this%calc_surface_area(n)
2269 dh = this%stage0(n) - this%stage(n)
2272 if (this%gwfiss == 0)
then
2273 r = this%usflow0(n) - this%usflow(n)
2281 if (this%imover == 1)
then
2282 q = this%pakmvrobj%get_qfrommvr(n)
2283 q0 = this%pakmvrobj%get_qfrommvr0(n)
2284 dqfrommvr = qtolfact * (q0 - q)
2293 dqfrommvrmax = dqfrommvr
2296 if (abs(dh) > abs(dhmax))
then
2300 if (abs(r) > abs(rmax))
then
2304 if (abs(dqfrommvr) > abs(dqfrommvrmax))
then
2305 dqfrommvrmax = dqfrommvr
2312 if (abs(dhmax) > abs(dpak))
then
2315 write (cloc,
"(a,'-',a)") trim(this%packName),
'stage'
2318 if (abs(rmax) > abs(dpak))
then
2321 write (cloc,
"(a,'-',a)") trim(this%packName),
'inflow'
2324 if (this%imover == 1)
then
2325 if (abs(dqfrommvrmax) > abs(dpak))
then
2326 ipak = locdqfrommvrmax
2328 write (cloc,
"(a,'-',a)") trim(this%packName),
'qfrommvr'
2334 if (this%ipakcsv /= 0)
then
2337 call this%pakcsvtab%add_term(innertot)
2338 call this%pakcsvtab%add_term(
totim)
2339 call this%pakcsvtab%add_term(
kper)
2340 call this%pakcsvtab%add_term(
kstp)
2341 call this%pakcsvtab%add_term(kiter)
2342 call this%pakcsvtab%add_term(dhmax)
2343 call this%pakcsvtab%add_term(locdhmax)
2344 call this%pakcsvtab%add_term(rmax)
2345 call this%pakcsvtab%add_term(locrmax)
2346 if (this%imover == 1)
then
2347 call this%pakcsvtab%add_term(dqfrommvrmax)
2348 call this%pakcsvtab%add_term(locdqfrommvrmax)
2353 call this%pakcsvtab%finalize_table()
2371 class(
sfrtype),
intent(inout) :: this
2372 real(DP),
dimension(:),
intent(in) :: x
2373 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
2374 integer(I4B),
optional,
intent(in) :: iadv
2381 real(DP) :: qoutflow
2382 real(DP) :: qfrommvr
2387 call this%BndType%bnd_cq(x, flowja, iadv=1)
2390 do n = 1, this%maxbound
2395 if (this%imover == 1)
then
2396 qfrommvr = this%pakmvrobj%get_qfrommvr(n)
2397 qtomvr = this%pakmvrobj%get_qtomvr(n)
2398 if (qtomvr >
dzero)
then
2404 qext = this%dsflow(n)
2406 if (qext >
dzero)
then
2409 do i = this%ia(n) + 1, this%ia(n + 1) - 1
2410 if (this%idir(i) > 0) cycle
2412 if (this%iboundpak(n2) == 0) cycle
2418 if (qext <
dzero)
then
2419 if (qtomvr <
dzero)
then
2420 qext = qext - qtomvr
2423 qoutflow = this%dsflow(n)
2424 if (qoutflow >
dzero)
then
2425 qoutflow = -qoutflow
2431 this%qextoutflow(n) = qext
2432 this%qoutflow(n) = qoutflow
2437 call this%sfr_fill_budobj()
2453 integer(I4B),
intent(in) :: icbcfl
2454 integer(I4B),
intent(in) :: ibudfl
2456 integer(I4B) :: ibinun
2457 character(len=20),
dimension(:),
allocatable :: cellidstr
2459 integer(I4B) :: node
2463 if (this%ibudgetout /= 0)
then
2464 ibinun = this%ibudgetout
2466 if (icbcfl == 0) ibinun = 0
2467 if (ibinun > 0)
then
2468 call this%budobj%save_flows(this%dis, ibinun,
kstp,
kper,
delt, &
2473 if (ibudfl /= 0 .and. this%iprflow /= 0)
then
2479 if (this%ianynone > 0)
then
2480 allocate (cellidstr(this%maxbound))
2481 do n = 1, this%maxbound
2482 node = this%igwfnode(n)
2484 call this%dis%noder_to_string(node, cellidstr(n))
2486 cellidstr(n) =
'NONE'
2489 call this%budobj%write_flowtable(this%dis,
kstp,
kper, cellidstr)
2490 deallocate (cellidstr)
2492 call this%budobj%write_flowtable(this%dis,
kstp,
kper)
2511 integer(I4B),
intent(in) :: idvsave
2512 integer(I4B),
intent(in) :: idvprint
2514 character(len=20) :: cellid
2515 integer(I4B) :: ibinun
2517 integer(I4B) :: node
2530 if (this%istageout /= 0)
then
2531 ibinun = this%istageout
2533 if (idvsave == 0) ibinun = 0
2536 if (ibinun > 0)
then
2537 do n = 1, this%maxbound
2540 if (this%iboundpak(n) == 0)
then
2542 else if (d ==
dzero)
then
2548 this%maxbound, 1, 1, ibinun)
2552 if (idvprint /= 0 .and. this%iprhed /= 0)
then
2555 call this%stagetab%set_kstpkper(
kstp,
kper)
2558 do n = 1, this%maxbound
2559 node = this%igwfnode(n)
2561 call this%dis%noder_to_string(node, cellid)
2562 hgwf = this%xnew(node)
2566 if (this%inamedbound == 1)
then
2567 call this%stagetab%add_term(this%boundname(n))
2569 call this%stagetab%add_term(n)
2570 call this%stagetab%add_term(cellid)
2571 if (this%iboundpak(n) /= 0)
then
2572 depth = this%depth(n)
2573 stage = this%stage(n)
2574 w = this%calc_top_width_wet(n, depth)
2575 call this%sfr_calc_cond(n, depth, cond, stage, hgwf)
2582 if (depth ==
dzero)
then
2583 call this%stagetab%add_term(
dhdry)
2585 call this%stagetab%add_term(stage)
2587 call this%stagetab%add_term(depth)
2588 call this%stagetab%add_term(w)
2590 if (this%iboundpak(n) /= 0)
then
2591 sbot = this%strtop(n) - this%bthick(n)
2592 if (hgwf < sbot)
then
2597 grad = grad / this%bthick(n)
2601 call this%stagetab%add_term(hgwf)
2602 call this%stagetab%add_term(cond)
2603 call this%stagetab%add_term(grad)
2605 call this%stagetab%add_term(
'--')
2606 call this%stagetab%add_term(
'--')
2607 call this%stagetab%add_term(
'--')
2626 integer(I4B),
intent(in) :: kstp
2627 integer(I4B),
intent(in) :: kper
2628 integer(I4B),
intent(in) :: iout
2629 integer(I4B),
intent(in) :: ibudfl
2631 call this%budobj%write_budtable(kstp, kper, iout, ibudfl,
totim,
delt)
2651 deallocate (this%csfrbudget)
2654 deallocate (this%cauxcbc)
2702 if (
associated(this%divcprior))
then
2703 deallocate (this%divcprior)
2717 call this%budobj%budgetobject_da()
2718 deallocate (this%budobj)
2719 nullify (this%budobj)
2722 if (this%iprhed > 0)
then
2723 call this%stagetab%table_da()
2724 deallocate (this%stagetab)
2725 nullify (this%stagetab)
2729 if (this%ipakcsv > 0)
then
2730 if (
associated(this%pakcsvtab))
then
2731 call this%pakcsvtab%table_da()
2732 deallocate (this%pakcsvtab)
2733 nullify (this%pakcsvtab)
2759 nullify (this%gwfiss)
2762 call this%BndType%bnd_da()
2776 class(
sfrtype),
intent(inout) :: this
2779 this%listlabel = trim(this%filtyp)//
' NO.'
2780 if (this%dis%ndim == 3)
then
2781 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
2782 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'ROW'
2783 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'COL'
2784 elseif (this%dis%ndim == 2)
then
2785 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
2786 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'CELL2D'
2788 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'NODE'
2790 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'STRESS RATE'
2791 if (this%inamedbound == 1)
then
2792 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'BOUNDARY NAME'
2830 integer(I4B) :: indx
2834 call this%obs%StoreObsType(
'stage', .false., indx)
2839 call this%obs%StoreObsType(
'inflow', .true., indx)
2844 call this%obs%StoreObsType(
'ext-inflow', .true., indx)
2849 call this%obs%StoreObsType(
'rainfall', .true., indx)
2854 call this%obs%StoreObsType(
'runoff', .true., indx)
2859 call this%obs%StoreObsType(
'evaporation', .true., indx)
2864 call this%obs%StoreObsType(
'outflow', .true., indx)
2869 call this%obs%StoreObsType(
'ext-outflow', .true., indx)
2874 call this%obs%StoreObsType(
'to-mvr', .true., indx)
2879 call this%obs%StoreObsType(
'from-mvr', .true., indx)
2884 call this%obs%StoreObsType(
'sfr', .true., indx)
2889 call this%obs%StoreObsType(
'upstream-flow', .true., indx)
2894 call this%obs%StoreObsType(
'downstream-flow', .true., indx)
2899 call this%obs%StoreObsType(
'depth', .false., indx)
2904 call this%obs%StoreObsType(
'wet-perimeter', .false., indx)
2909 call this%obs%StoreObsType(
'wet-area', .false., indx)
2914 call this%obs%StoreObsType(
'wet-width', .false., indx)
2934 character(len=100) :: msg
2938 if (this%obs%npakobs > 0)
then
2939 call this%obs%obs_bd_clear()
2940 do i = 1, this%obs%npakobs
2941 obsrv => this%obs%pakobs(i)%obsrv
2942 do j = 1, obsrv%indxbnds_count
2943 n = obsrv%indxbnds(j)
2945 select case (obsrv%ObsTypeId)
2950 if (this%imover == 1)
then
2951 v = this%pakmvrobj%get_qtomvr(n)
2958 if (this%imover == 1)
then
2959 v = this%pakmvrobj%get_qfrommvr(n)
2966 v = this%qoutflow(n)
2967 case (
'EXT-OUTFLOW')
2968 v = this%qextoutflow(n)
2970 if (this%iboundpak(n) /= 0)
then
2976 v = this%simrunoff(n)
2977 case (
'EVAPORATION')
2981 case (
'UPSTREAM-FLOW')
2983 if (this%imover == 1)
then
2984 v = v + this%pakmvrobj%get_qfrommvr(n)
2986 case (
'DOWNSTREAM-FLOW')
2993 case (
'WET-PERIMETER')
2994 v = this%calc_perimeter_wet(n, this%depth(n))
2996 v = this%calc_area_wet(n, this%depth(n))
2998 v = this%calc_top_width_wet(n, this%depth(n))
3000 msg =
'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
3003 call this%obs%SaveOneSimval(obsrv, v)
3009 call this%parser%StoreErrorUnit()
3026 class(
sfrtype),
intent(inout) :: this
3031 character(len=LENBOUNDNAME) :: bname
3032 logical(LGP) :: jfound
3035 10
format(
'Boundary "', a,
'" for observation "', a, &
3036 '" is invalid in package "', a,
'"')
3037 30
format(
'Boundary name not provided for observation "', a, &
3038 '" in package "', a,
'"')
3044 do i = 1, this%obs%npakobs
3045 obsrv => this%obs%pakobs(i)%obsrv
3048 nn1 = obsrv%NodeNumber
3050 bname = obsrv%FeatureName
3051 if (bname /=
'')
then
3056 do j = 1, this%maxbound
3057 if (this%boundname(j) == bname)
then
3059 call obsrv%AddObsIndex(j)
3062 if (.not. jfound)
then
3064 trim(bname), trim(obsrv%name), trim(this%packName)
3068 write (
errmsg, 30) trim(obsrv%name), trim(this%packName)
3071 else if (nn1 < 1 .or. nn1 > this%maxbound)
then
3072 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
3073 trim(adjustl(obsrv%ObsTypeId)), &
3074 'reach must be greater than 0 and less than or equal to', &
3075 this%maxbound,
'(specified value is ', nn1,
')'
3078 if (obsrv%indxbnds_count == 0)
then
3079 call obsrv%AddObsIndex(nn1)
3081 errmsg =
'Programming error in sfr_rp_obs'
3088 if (obsrv%ObsTypeId ==
'STAGE' .or. &
3089 obsrv%ObsTypeId ==
'DEPTH' .or. &
3090 obsrv%ObsTypeId ==
'WET-PERIMETER' .or. &
3091 obsrv%ObsTypeId ==
'WET-AREA' .or. &
3092 obsrv%ObsTypeId ==
'WET-WIDTH')
then
3093 nn1 = obsrv%NodeNumber
3095 if (obsrv%indxbnds_count > 1)
then
3096 write (
errmsg,
'(a,3(1x,a))') &
3097 trim(adjustl(obsrv%ObsTypeId)), &
3098 'for observation', trim(adjustl(obsrv%Name)), &
3099 ' must be assigned to a reach with a unique boundname.'
3106 do j = 1, obsrv%indxbnds_count
3107 nn1 = obsrv%indxbnds(j)
3108 if (nn1 < 1 .or. nn1 > this%maxbound)
then
3109 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
3110 trim(adjustl(obsrv%ObsTypeId)), &
3111 'reach must be greater than 0 and less than or equal to', &
3112 this%maxbound,
'(specified value is ', nn1,
')'
3120 call this%parser%StoreErrorUnit()
3140 integer(I4B),
intent(in) :: inunitobs
3141 integer(I4B),
intent(in) :: iout
3144 integer(I4B) :: icol
3145 integer(I4B) :: istart
3146 integer(I4B) :: istop
3147 character(len=LINELENGTH) :: string
3148 character(len=LENBOUNDNAME) :: bndname
3151 string = obsrv%IDstring
3161 obsrv%FeatureName = bndname
3165 obsrv%NodeNumber = nn1
3184 class(
sfrtype),
intent(inout) :: this
3185 integer(I4B),
intent(in) :: n
3186 integer(I4B),
intent(inout) :: ichkustrm
3187 character(len=LINELENGTH),
intent(inout) :: crossfile
3189 character(len=10) :: cnum
3190 character(len=LINELENGTH) :: text
3191 character(len=LINELENGTH) :: caux
3192 character(len=LINELENGTH) :: keyword
3193 integer(I4B) :: ival
3196 integer(I4B) :: idiv
3197 integer(I4B) :: ixserror
3198 character(len=10) :: cp
3200 real(DP),
pointer :: bndElem => null()
3206 call this%parser%GetStringCaps(keyword)
3207 select case (keyword)
3210 call this%parser%GetStringCaps(text)
3211 if (text ==
'INACTIVE')
then
3212 this%iboundpak(n) = 0
3213 else if (text ==
'ACTIVE')
then
3214 this%iboundpak(n) = 1
3215 else if (text ==
'SIMPLE')
then
3216 this%iboundpak(n) = -1
3219 'Unknown '//trim(this%text)//
' sfr status keyword: ', trim(text)
3223 call this%parser%GetString(text)
3225 bndelem => this%hk(n)
3227 this%packName,
'BND', &
3228 this%tsManager, this%iprpak, &
3231 call this%parser%GetString(text)
3233 bndelem => this%rough(n)
3235 this%packName,
'BND', &
3236 this%tsManager, this%iprpak, &
3239 call this%parser%GetString(text)
3241 bndelem => this%sstage(n)
3243 this%packName,
'BND', &
3244 this%tsManager, this%iprpak,
'STAGE')
3246 call this%parser%GetString(text)
3248 bndelem => this%rain(n)
3250 this%packName,
'BND', &
3251 this%tsManager, this%iprpak,
'RAIN')
3252 case (
'EVAPORATION')
3253 call this%parser%GetString(text)
3255 bndelem => this%evap(n)
3257 this%packName,
'BND', &
3258 this%tsManager, this%iprpak, &
3261 call this%parser%GetString(text)
3263 bndelem => this%runoff(n)
3265 this%packName,
'BND', &
3266 this%tsManager, this%iprpak, &
3269 call this%parser%GetString(text)
3271 bndelem => this%inflow(n)
3273 this%packName,
'BND', &
3274 this%tsManager, this%iprpak, &
3279 if (this%ndiv(n) < 1)
then
3280 write (cnum,
'(i0)') n
3281 errmsg =
'diversions cannot be specified for reach '//trim(cnum)
3286 ival = this%parser%GetInteger()
3287 if (ival < 1 .or. ival > this%ndiv(n))
then
3288 write (cnum,
'(i0)') n
3289 errmsg =
'Reach '//trim(cnum)
3290 write (cnum,
'(i0)') this%ndiv(n)
3291 errmsg = trim(
errmsg)//
' diversion number should be between 1 '// &
3292 'and '//trim(cnum)//
'.'
3298 call this%parser%GetString(text)
3299 ii = this%iadiv(n) + idiv - 1
3301 bndelem => this%divflow(ii)
3303 this%packName,
'BND', &
3304 this%tsManager, this%iprpak, &
3308 cp = this%divcprior(ii)
3309 divq = this%divflow(ii)
3310 if (cp ==
'FRACTION' .and. (divq <
dzero .or. divq >
done))
then
3311 write (
errmsg,
'(a,1x,i0,a)') &
3312 'cprior is type FRACTION for diversion no.', ii, &
3313 ', but divflow not within the range 0.0 to 1.0'
3316 case (
'UPSTREAM_FRACTION')
3318 call this%parser%GetString(text)
3320 bndelem => this%ustrf(n)
3322 this%packName,
'BND', &
3323 this%tsManager, this%iprpak,
'USTRF')
3325 case (
'CROSS_SECTION')
3329 call this%parser%GetStringCaps(keyword)
3330 select case (keyword)
3332 call this%parser%GetStringCaps(keyword)
3333 if (trim(adjustl(keyword)) /=
'FILEIN')
then
3334 errmsg =
'TAB6 keyword must be followed by "FILEIN" '// &
3339 if (ixserror == 0)
then
3340 call this%parser%GetString(crossfile)
3343 write (
errmsg,
'(a,1x,i4,1x,a)') &
3344 'CROSS-SECTION TABLE ENTRY for REACH ', n, &
3345 'MUST INCLUDE TAB6 KEYWORD'
3350 call this%parser%GetStringCaps(caux)
3351 do jj = 1, this%naux
3352 if (trim(adjustl(caux)) /= trim(adjustl(this%auxname(jj)))) cycle
3353 call this%parser%GetString(text)
3355 bndelem => this%rauxvar(jj, ii)
3357 this%packName,
'AUX', &
3358 this%tsManager, this%iprpak, &
3364 write (
errmsg,
'(a,a)') &
3365 'Unknown '//trim(this%text)//
' sfr data keyword: ', &
3382 integer(I4B),
intent(in) :: n
3383 real(DP),
intent(in) :: h
3384 real(DP),
intent(inout) :: hcof
3385 real(DP),
intent(inout) :: rhs
3386 logical(LGP),
intent(in),
optional :: update
3388 logical(LGP) :: lupdate
3401 real(DP) :: qfrommvr
3414 if (
present(update))
then
3424 if (this%iboundpak(n) == 0)
then
3425 this%depth(n) =
dzero
3427 this%usflow(n) =
dzero
3428 this%simevap(n) =
dzero
3429 this%simrunoff(n) =
dzero
3430 this%dsflow(n) =
dzero
3431 this%gwflow(n) =
dzero
3442 do i = this%ia(n) + 1, this%ia(n + 1) - 1
3443 if (this%idir(i) < 0) cycle
3445 do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
3446 if (this%idir(ii) > 0) cycle
3447 if (this%ja(ii) /= n) cycle
3448 qu = qu + this%qconn(ii)
3454 sa = this%calc_surface_area(n)
3455 sa_wet = this%calc_surface_area_wet(n, this%depth(n))
3457 qr = this%rain(n) * sa
3458 qe = this%evap(n) * sa_wet
3459 qro = this%runoff(n)
3463 if (this%imover == 1)
then
3464 qfrommvr = this%pakmvrobj%get_qfrommvr(n)
3468 qsrc = qu + qi + qr - qe + qro + qfrommvr
3471 call this%sfr_adjust_ro_ev(qsrc, qu, qi, qr, qro, qe, qfrommvr)
3474 this%simevap(n) = qe
3475 this%simrunoff(n) = qro
3478 if (this%iboundpak(n) < 0)
then
3479 call this%sfr_calc_constant(n, d1, hgwf, qgwf, qd)
3481 call this%sfr_calc_steady(n, d1, hgwf, qu, qi, &
3482 qfrommvr, qr, qe, qro, &
3488 bt = tp - this%bthick(n)
3495 this%stage(n) = hsfr
3497 call this%sfr_update_flows(n, qd, qgwf)
3503 if (this%gwfiss == 0)
then
3514 call this%sfr_calc_qgwf(n, d1, hgwf, qgwf, gwfhcof, gwfrhs)
3517 if (abs(sumleak) >
dzero)
then
3522 else if ((sumleak - qsrc) < -
dem30)
then
3523 if (this%gwfiss == 0)
then
3524 rhs = rhs + gwfrhs - sumrch
3531 if (this%gwfiss == 0)
then
3532 rhs = rhs - sumleak - sumrch
3539 else if (hgwf < bt)
then
3553 class(
sfrtype),
intent(inout) :: this
3554 integer(I4B),
intent(in) :: n
3555 real(DP),
intent(inout) :: qd
3556 real(DP),
intent(in) :: qgwf
3560 integer(I4B) :: idiv
3561 integer(I4B) :: jpos
3571 this%gwflow(n) = qgwf
3574 if (qd >
dzero)
then
3577 do i = this%ia(n) + 1, this%ia(n + 1) - 1
3578 if (this%idir(i) > 0) cycle
3580 if (idiv == 0) cycle
3581 jpos = this%iadiv(n) + idiv - 1
3582 call this%sfr_calc_div(n, idiv, qd, qdiv)
3583 this%qconn(i) = qdiv
3584 this%divq(jpos) = qdiv
3590 if (this%imover == 1)
then
3591 call this%pakmvrobj%accumulate_qformvr(n, qd)
3592 qd = max(qd - this%pakmvrobj%get_qtomvr(n),
dzero)
3596 do i = this%ia(n) + 1, this%ia(n + 1) - 1
3597 if (this%idir(i) > 0) cycle
3598 if (this%idiv(i) > 0) cycle
3600 if (this%iboundpak(n2) == 0) cycle
3601 f = this%ustrf(n2) / this%ftotnd(n)
3602 this%qconn(i) = qd * f
3605 do i = this%ia(n) + 1, this%ia(n + 1) - 1
3606 if (this%idir(i) > 0) cycle
3607 this%qconn(i) =
dzero
3609 if (idiv == 0) cycle
3610 jpos = this%iadiv(n) + idiv - 1
3611 this%divq(jpos) =
dzero
3628 real(DP),
intent(inout) :: qc
3629 real(DP),
intent(in) :: qu
3630 real(DP),
intent(in) :: qi
3631 real(DP),
intent(in) :: qr
3632 real(DP),
intent(inout) :: qro
3633 real(DP),
intent(inout) :: qe
3634 real(DP),
intent(in) :: qfrommvr
3639 if (qc <
dzero)
then
3642 qt = qu + qi + qr + qro + qfrommvr
3645 if (qt <
dzero)
then
3646 if (qro <
dzero)
then
3647 qro = -(qu + qi + qr + qfrommvr)
3653 if (qe >
dzero)
then
3654 qe = qu + qi + qr + qro + qfrommvr
3657 qc = qu + qi + qr - qe + qro + qfrommvr
3669 integer(I4B),
intent(in) :: n
3670 real(DP),
intent(in) :: depth
3671 real(DP),
intent(in) :: hgwf
3672 real(DP),
intent(inout) :: qgwf
3673 real(DP),
intent(inout) :: qd
3681 call this%sfr_calc_qsource(n, depth, qsrc)
3684 call this%sfr_calc_qgwf(n, depth, hgwf, qgwf)
3685 if (-qgwf > qsrc) qgwf = -qsrc
3706 integer(I4B),
intent(in) :: n
3707 real(DP),
intent(in) :: depth
3708 real(DP),
intent(inout) :: qsrc
3715 real(DP) :: qfrommvr
3725 qro = this%runoff(n)
3728 a = this%calc_surface_area(n)
3729 ae = this%calc_surface_area_wet(n, depth)
3730 qr = this%rain(n) * a
3731 qe = this%evap(n) * ae
3735 if (this%imover == 1)
then
3736 qfrommvr = this%pakmvrobj%get_qfrommvr(n)
3740 qsrc = qu + qi + qr - qe + qro + qfrommvr
3743 call this%sfr_adjust_ro_ev(qsrc, qu, qi, qr, qro, qe, qfrommvr)
3758 integer(I4B),
intent(in) :: n
3759 real(DP),
intent(in) :: depth
3760 real(DP),
intent(inout) :: qman
3762 integer(I4B) :: npts
3777 if (depth >
dzero)
then
3778 npts = this%ncrosspts(n)
3789 i0 = this%iacross(n)
3790 i1 = this%iacross(n + 1) - 1
3795 this%station(i0:i1), &
3796 this%xsheight(i0:i1), &
3797 this%xsrough(i0:i1), &
3804 aw = this%calc_area_wet(n, depth)
3805 wp = this%calc_perimeter_wet(n, depth)
3806 if (wp >
dzero)
then
3811 qman = this%unitconv * aw * (rh**
dtwothirds) * sqrt(s) / r
3832 integer(I4B),
intent(in) :: n
3833 real(DP),
intent(in) :: depth
3834 real(DP),
intent(in) :: hgwf
3835 real(DP),
intent(inout) :: qgwf
3836 real(DP),
intent(inout),
optional :: gwfhcof
3837 real(DP),
intent(inout),
optional :: gwfrhs
3839 integer(I4B) :: node
3847 real(DP) :: gwfhcof0
3854 node = this%igwfnode(n)
3855 if (node < 1)
return
3858 if (this%ibound(node) == 0)
return
3865 bt = tp - this%bthick(n)
3868 if (h_temp < bt)
then
3873 call this%sfr_calc_cond(n, depth, cond, hsfr, h_temp)
3876 qgwf = sat * cond * (h_temp - hsfr)
3877 gwfrhs0 = -sat * cond * hsfr
3878 gwfhcof0 = -sat * cond
3881 if (this%idense /= 0)
then
3882 call this%sfr_calculate_density_exchange(n, hsfr, hgwf, cond, tp, &
3883 qgwf, gwfhcof0, gwfrhs0)
3887 if (
present(gwfhcof)) gwfhcof = gwfhcof0
3888 if (
present(gwfrhs)) gwfrhs = gwfrhs0
3905 integer(I4B),
intent(in) :: n
3907 integer(I4B) :: node
3910 node = this%igwfnode(n)
3911 if (node > 0 .and. this%hk(n) >
dzero)
then
3924 integer(I4B),
intent(in) :: n
3925 real(DP),
intent(in) :: depth
3926 real(DP),
intent(inout) :: cond
3927 real(DP),
intent(in),
optional :: hsfr
3928 real(DP),
intent(in),
optional :: h_temp
3930 integer(I4B) :: node
3932 real(DP) :: vscratio
3942 node = this%igwfnode(n)
3944 if (this%ibound(node) > 0)
then
3947 if (this%ivsc == 1)
then
3948 if (hsfr > h_temp)
then
3950 vscratio = this%viscratios(1, n)
3952 vscratio = this%viscratios(2, n)
3955 wp = this%calc_perimeter_wet(n, depth)
3956 cond = this%hk(n) * vscratio * this%length(n) * wp / this%bthick(n)
3967 integer(I4B),
intent(in) :: n
3968 real(DP),
intent(inout) :: d1
3969 real(DP),
intent(in) :: hgwf
3970 real(DP),
intent(inout) :: qgwf
3971 real(DP),
intent(inout) :: qd
3975 this%stage(n) = this%sstage(n)
3976 d1 = max(
dzero, this%stage(n) - this%strtop(n))
3978 call this%sfr_calc_qsource(n, d1, qsrc)
3980 if (this%sfr_gwf_conn(n) == 1)
then
3981 call this%sfr_calc_qgwf(n, d1, hgwf, qgwf)
3988 if (qgwf > qsrc)
then
3990 call this%sfr_calc_qsource(n, d1, qsrc)
4000 qu, qi, qfrommvr, qr, qe, qro, &
4004 integer(I4B),
intent(in) :: n
4005 real(DP),
intent(inout) :: d1
4006 real(DP),
intent(in) :: hgwf
4007 real(DP),
intent(in) :: qu
4008 real(DP),
intent(in) :: qi
4009 real(DP),
intent(in) :: qfrommvr
4010 real(DP),
intent(in) :: qr
4011 real(DP),
intent(in) :: qe
4012 real(DP),
intent(in) :: qro
4013 real(DP),
intent(inout) :: qgwf
4014 real(DP),
intent(inout) :: qd
4017 integer(I4B) :: isolve
4019 integer(I4B) :: iic2
4020 integer(I4B) :: iic3
4021 integer(I4B) :: iic4
4022 integer(I4B) :: ibflg
4063 qc = qu + qi + qr - qe + qro + qfrommvr
4067 qsrcmp = qu + qi + qfrommvr +
dhalf * (qr - qe + qro)
4071 call this%sfr_calc_reach_depth(n, qmp, d1)
4075 call this%sfr_calc_qsource(n, d1, qsrc)
4080 bt = tp - this%bthick(n)
4082 qd = max(qsrc,
dzero)
4087 if (hsfr <= tp .and. hgwf <= tp) isolve = 0
4088 if (hgwf <= tp .and. qc <
dem30) isolve = 0
4089 if (this%sfr_gwf_conn(n) == 0) isolve = 0
4092 calc_solution:
if (isolve /= 0)
then
4096 if (d1 >
dem30)
then
4097 if ((tp - hgwf) >
dem30)
then
4100 en2 =
d1p1 * d1 - (tp - hgwf)
4102 else if ((tp - hgwf) >
dem30)
then
4105 en2 =
dp99 * (hgwf - tp)
4111 call this%sfr_calc_qgwf(n,
dzero, hgwf, qgwf1)
4113 qen1 = qmp -
dhalf * qgwf1
4119 call this%sfr_calc_qgwf(n, en2, hgwf, qgwf2)
4122 call this%sfr_calc_qgwf(n, en2, bt, qgwf2)
4125 if (qgwf2 > qsrc) qgwf2 = qsrc
4127 call this%sfr_calc_reach_depth(n, (qsrcmp -
dhalf * qgwf1), d1)
4128 call this%sfr_calc_reach_depth(n, (qsrcmp -
dhalf * qgwf2), d2)
4130 if (d1 >
dem30)
then
4136 if (d2 >
dem30)
then
4138 if (f2 <
dem30) en2 = d2
4145 dpp =
dhalf * (en1 + en2)
4154 do i = 1, this%maxsfrit
4157 d2 = d1 +
dtwo * this%deps
4159 call this%sfr_calc_qman(n, d1, q1)
4160 call this%sfr_calc_qman(n, d2, q2)
4162 call this%sfr_calc_qgwf(n, d1, hgwf, qgwf1)
4164 call this%sfr_calc_qgwf(n, d2, hgwf, qgwf2)
4167 if (qgwf1 >= qsrc)
then
4169 dpp =
dhalf * (en1 + en2)
4170 call this%sfr_calc_qgwf(n, dpp, hgwf, qgwfp)
4172 if (qgwfp > qsrc) qgwfp = qsrc
4173 call this%sfr_calc_reach_depth(n, (qsrcmp -
dhalf * qgwfp), dx)
4176 fhstr1 = (qsrcmp -
dhalf * qgwf1) - q1
4177 fhstr2 = (qsrcmp -
dhalf * qgwf2) - q2
4180 if (ibflg == 0)
then
4182 if (abs(d1 - d2) >
dzero)
then
4183 derv = (fhstr1 - fhstr2) / (d1 - d2)
4185 if (abs(derv) >
dem30)
then
4186 dlh = -fhstr1 / derv
4193 if ((dpp >= en2) .or. (dpp <= en1))
then
4194 if (abs(dlh) > abs(dlhold) .or. dpp <
dem30)
then
4196 dpp =
dhalf * (en1 + en2)
4203 if (qgwf1 * qgwfold <
dem30)
then
4208 if (qgwf1 <
dem30)
then
4213 if (dlh * dlhold <
dem30 .or. abs(dlh) > abs(dlhold))
then
4217 if (iic3 > 7 .and. iic > 12)
then
4223 if (iic2 > 7 .or. iic > 12 .or. iic4 == 1)
then
4225 dpp =
dhalf * (en1 + en2)
4229 call this%sfr_calc_qgwf(n, dpp, hgwf, qgwfp)
4231 if (qgwfp > qsrc)
then
4233 if (abs(en1 - en2) < this%dmaxchg *
dem6)
then
4234 call this%sfr_calc_reach_depth(n, (qsrcmp -
dhalf * qgwfp), dpp)
4237 call this%sfr_calc_reach_depth(n, (qsrcmp -
dhalf * qgwfp), dx)
4242 if (ibflg == 1)
then
4246 if (f1 * fp <
dzero)
then
4254 err = min(abs(fp), abs(en2 - en1))
4260 if (err < this%dmaxchg)
then
4270 if (ibflg == 1)
then
4280 call this%sfr_calc_qgwf(n, d1, hgwf, qgwf)
4284 if (qgwf > qsrc)
then
4286 call this%sfr_calc_qsource(n, d1, qsrc)
4291 end if calc_solution
4306 integer(I4B),
intent(in) :: n
4307 integer(I4B),
intent(in) :: i
4308 real(DP),
intent(inout) :: qd
4309 real(DP),
intent(inout) :: qdiv
4311 character(len=10) :: cp
4312 integer(I4B) :: jpos
4317 jpos = this%iadiv(n) + i - 1
4318 n2 = this%divreach(jpos)
4319 cp = this%divcprior(jpos)
4320 v = this%divflow(jpos)
4364 integer(I4B),
intent(in) :: n
4365 real(DP),
intent(in) :: q1
4366 real(DP),
intent(inout) :: d1
4378 if (q1 >
dzero)
then
4379 if (this%ncrosspts(n) > 1)
then
4380 call this%sfr_calc_xs_depth(n, q1, d1)
4382 w = this%station(this%iacross(n))
4383 qconst = this%unitconv * w * sqrt(s) / r
4384 d1 = (q1 / qconst)**
dp6
4403 integer(I4B),
intent(in) :: n
4404 real(DP),
intent(in) :: qrch
4405 real(DP),
intent(inout) :: d
4407 integer(I4B) :: iter
4408 real(DP) :: perturbation
4414 real(DP) :: residual
4417 perturbation = this%deps *
dtwo
4420 residual = q0 - qrch
4423 nriter:
do iter = 1, this%maxsfrit
4424 call this%sfr_calc_qman(n, d + perturbation, q1)
4426 if (dq /=
dzero)
then
4427 derv = perturbation / (q1 - q0)
4431 dd = derv * residual
4433 call this%sfr_calc_qman(n, d, q0)
4434 residual = q0 - qrch
4437 if (abs(dd) < this%dmaxchg)
then
4457 character(len=*),
parameter :: fmtunitconv_error = &
4458 &
"('SFR (',a,') UNIT_CONVERSION SPECIFIED VALUE (',g0,') AND', &
4459 &1x,'LENGTH_CONVERSION OR TIME_CONVERSION SPECIFIED.')"
4460 character(len=*),
parameter :: fmtunitconv = &
4461 &
"(1x,'SFR PACKAGE (',a,') CONVERSION DATA',&
4462 &/4x,'UNIT CONVERSION VALUE (',g0,').',/)"
4465 if (this%lengthconv /=
dnodata .or. this%timeconv /=
dnodata)
then
4466 if (this%unitconv /=
done)
then
4467 write (
errmsg, fmtunitconv_error) &
4468 trim(adjustl(this%packName)), this%unitconv
4471 if (this%lengthconv /=
dnodata)
then
4472 this%unitconv = this%unitconv * this%lengthconv**
donethird
4474 if (this%timeconv /=
dnodata)
then
4475 this%unitconv = this%unitconv * this%timeconv
4477 write (this%iout, fmtunitconv) &
4478 trim(adjustl(this%packName)), this%unitconv
4497 character(len=5) :: crch
4498 character(len=10) :: cval
4499 character(len=30) :: nodestr
4500 character(len=LINELENGTH) :: title
4501 character(len=LINELENGTH) :: text
4508 if (this%iprpak /= 0)
then
4509 title = trim(adjustl(this%text))//
' PACKAGE ('// &
4510 trim(adjustl(this%packName))//
') STATIC REACH DATA'
4511 call table_cr(this%inputtab, this%packName, title)
4512 call this%inputtab%table_df(this%maxbound, 10, this%iout)
4514 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4516 call this%inputtab%initialize_column(text, 20, alignment=
tableft)
4518 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4520 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4522 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4524 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4526 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4528 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4530 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4531 text =
'UPSTREAM FRACTION'
4532 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
4538 do n = 1, this%maxbound
4539 write (crch,
'(i5)') n
4540 nn = this%igwfnode(n)
4542 btgwf = this%dis%bot(nn)
4543 call this%dis%noder_to_string(nn, nodestr)
4548 if (this%length(n) <=
dzero)
then
4549 errmsg =
'Reach '//crch//
' length must be greater than 0.0.'
4553 if (this%width(n) <=
dzero)
then
4554 errmsg =
'Reach '//crch//
' width must be greater than 0.0.'
4558 if (this%slope(n) <=
dzero)
then
4559 errmsg =
'Reach '//crch//
' slope must be greater than 0.0.'
4564 bt = this%strtop(n) - this%bthick(n)
4565 if (bt <= btgwf .and. this%icheck /= 0)
then
4566 write (cval,
'(f10.4)') bt
4567 errmsg =
'Reach '//crch//
' bed bottom (rtp-rbth ='// &
4568 cval//
') must be greater than the bottom of cell ('// &
4570 write (cval,
'(f10.4)') btgwf
4574 if (this%hk(n) <
dzero)
then
4575 errmsg =
'Reach '//crch//
' hk must be greater than or equal to 0.0.'
4580 if (this%rough(n) <=
dzero)
then
4581 errmsg =
'Reach '//crch//
" Manning's roughness "// &
4582 'coefficient must be greater than 0.0.'
4586 if (this%ustrf(n) <
dzero)
then
4587 errmsg =
'Reach '//crch//
' upstream fraction must be greater '// &
4588 'than or equal to 0.0.'
4592 if (this%iprpak /= 0)
then
4593 call this%inputtab%add_term(n)
4594 call this%inputtab%add_term(nodestr)
4595 call this%inputtab%add_term(this%length(n))
4596 call this%inputtab%add_term(this%width(n))
4597 call this%inputtab%add_term(this%slope(n))
4598 call this%inputtab%add_term(this%strtop(n))
4599 call this%inputtab%add_term(this%bthick(n))
4600 call this%inputtab%add_term(this%hk(n))
4601 call this%inputtab%add_term(this%rough(n))
4602 call this%inputtab%add_term(this%ustrf(n))
4621 logical(LGP) :: lreorder
4622 character(len=5) :: crch
4623 character(len=5) :: crch2
4624 character(len=LINELENGTH) :: text
4625 character(len=LINELENGTH) :: title
4632 integer(I4B) :: ifound
4633 integer(I4B) :: ierr
4634 integer(I4B) :: maxconn
4635 integer(I4B) :: ntabcol
4639 do j = 1, this%MAXBOUND
4640 n = this%isfrorder(j)
4649 write (this%iout,
'(/,1x,a)') &
4650 trim(adjustl(this%text))//
' PACKAGE ('// &
4651 trim(adjustl(this%packName))//
') REACH SOLUTION HAS BEEN '// &
4652 'REORDERED USING A DAG'
4655 if (this%iprpak /= 0)
then
4659 title = trim(adjustl(this%text))//
' PACKAGE ('// &
4660 trim(adjustl(this%packName))//
') REACH SOLUTION ORDER'
4661 call table_cr(this%inputtab, this%packName, title)
4662 call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4664 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4666 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4669 do j = 1, this%maxbound
4670 n = this%isfrorder(j)
4671 call this%inputtab%add_term(j)
4672 call this%inputtab%add_term(n)
4678 if (this%iprpak /= 0)
then
4682 do n = 1, this%maxbound
4683 maxconn = max(maxconn, this%nconnreach(n))
4685 ntabcol = 1 + maxconn
4688 title = trim(adjustl(this%text))//
' PACKAGE ('// &
4689 trim(adjustl(this%packName))//
') STATIC REACH CONNECTION DATA'
4690 call table_cr(this%inputtab, this%packName, title)
4691 call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4693 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4695 write (text,
'(a,1x,i6)')
'CONN', n
4696 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4702 do n = 1, this%MAXBOUND
4703 write (crch,
'(i5)') n
4704 eachconn:
do i = this%ia(n) + 1, this%ia(n + 1) - 1
4706 write (crch2,
'(i5)') nn
4708 connreach:
do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4715 if (ifound /= 1)
then
4716 errmsg =
'Reach '//crch//
' is connected to '// &
4717 'reach '//crch2//
' but reach '//crch2// &
4718 ' is not connected to reach '//crch//
'.'
4724 if (this%iprpak /= 0)
then
4725 call this%inputtab%add_term(n)
4726 do i = this%ia(n) + 1, this%ia(n + 1) - 1
4727 call this%inputtab%add_term(this%ja(i))
4729 nn = maxconn - this%nconnreach(n)
4731 call this%inputtab%add_term(
' ')
4740 do n = 1, this%maxbound
4741 write (crch,
'(i5)') n
4742 eachconnv:
do i = this%ia(n) + 1, this%ia(n + 1) - 1
4745 if (this%idir(i) < 0) cycle eachconnv
4747 write (crch2,
'(i5)') nn
4748 connreachv:
do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4750 if (this%idir(ii) < 0) cycle connreachv
4757 errmsg =
'Reach '//crch//
' is connected to '// &
4758 'reach '//crch2//
' but streamflow from reach '// &
4759 crch//
' to reach '//crch2//
' is not permitted.'
4769 call this%parser%StoreErrorUnit()
4774 do n = 1, this%maxbound
4775 write (crch,
'(i5)') n
4776 eachconnds:
do i = this%ia(n) + 1, this%ia(n + 1) - 1
4778 if (this%idir(i) > 0) cycle eachconnds
4779 write (crch2,
'(i5)') nn
4781 connreachds:
do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4784 if (this%idir(i) /= this%idir(ii))
then
4790 if (ifound /= 1)
then
4791 errmsg =
'Reach '//crch//
' downstream connected reach '// &
4792 'is reach '//crch2//
' but reach '//crch//
' is not'// &
4793 ' the upstream connected reach for reach '//crch2//
'.'
4800 if (this%iprpak /= 0)
then
4804 do n = 1, this%maxbound
4806 do i = this%ia(n) + 1, this%ia(n + 1) - 1
4807 if (this%idir(i) > 0)
then
4811 maxconn = max(maxconn, ii)
4813 ntabcol = 1 + maxconn
4816 title = trim(adjustl(this%text))//
' PACKAGE ('// &
4817 trim(adjustl(this%packName))//
') STATIC UPSTREAM REACH '// &
4819 call table_cr(this%inputtab, this%packName, title)
4820 call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4822 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4824 write (text,
'(a,1x,i6)')
'UPSTREAM CONN', n
4825 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4829 do n = 1, this%maxbound
4830 call this%inputtab%add_term(n)
4832 do i = this%ia(n) + 1, this%ia(n + 1) - 1
4833 if (this%idir(i) > 0)
then
4834 call this%inputtab%add_term(this%ja(i))
4840 call this%inputtab%add_term(
' ')
4846 do n = 1, this%maxbound
4848 do i = this%ia(n) + 1, this%ia(n + 1) - 1
4849 if (this%idir(i) < 0)
then
4853 maxconn = max(maxconn, ii)
4855 ntabcol = 1 + maxconn
4858 title = trim(adjustl(this%text))//
' PACKAGE ('// &
4859 trim(adjustl(this%packName))//
') STATIC DOWNSTREAM '// &
4860 'REACH CONNECTION DATA'
4861 call table_cr(this%inputtab, this%packName, title)
4862 call this%inputtab%table_df(this%maxbound, ntabcol, this%iout)
4864 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4866 write (text,
'(a,1x,i6)')
'DOWNSTREAM CONN', n
4867 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4871 do n = 1, this%maxbound
4872 call this%inputtab%add_term(n)
4874 do i = this%ia(n) + 1, this%ia(n + 1) - 1
4875 if (this%idir(i) < 0)
then
4876 call this%inputtab%add_term(this%ja(i))
4882 call this%inputtab%add_term(
' ')
4902 character(len=LINELENGTH) :: title
4903 character(len=LINELENGTH) :: text
4904 character(len=5) :: crch
4905 character(len=5) :: cdiv
4906 character(len=5) :: crch2
4907 character(len=10) :: cprior
4908 integer(I4B) :: maxdiv
4913 integer(I4B) :: idiv
4914 integer(I4B) :: ifound
4915 integer(I4B) :: jpos
4917 10
format(
'Diversion ', i0,
' of reach ', i0, &
4918 ' is invalid or has not been defined.')
4921 if (this%iprpak /= 0)
then
4925 do n = 1, this%maxbound
4926 maxdiv = maxdiv + this%ndiv(n)
4930 title = trim(adjustl(this%text))//
' PACKAGE ('// &
4931 trim(adjustl(this%packName))//
') REACH DIVERSION DATA'
4932 call table_cr(this%inputtab, this%packName, title)
4933 call this%inputtab%table_df(maxdiv, 4, this%iout)
4935 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4937 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4939 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4941 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
4945 do n = 1, this%maxbound
4946 if (this%ndiv(n) < 1) cycle
4947 write (crch,
'(i5)') n
4949 do idiv = 1, this%ndiv(n)
4952 jpos = this%iadiv(n) + idiv - 1
4955 write (cdiv,
'(i5)') idiv
4958 nn = this%divreach(jpos)
4959 write (crch2,
'(i5)') nn
4963 if (nn < 1 .or. nn > this%maxbound)
then
4964 write (
errmsg, 10) idiv, n
4968 connreach:
do ii = this%ia(nn) + 1, this%ia(nn + 1) - 1
4971 if (this%idir(ii) > 0)
then
4977 if (ifound /= 1)
then
4978 errmsg =
'Reach '//crch//
' is not a upstream reach for '// &
4979 'reach '//crch2//
' as a result diversion '//cdiv// &
4980 ' from reach '//crch//
' to reach '//crch2// &
4981 ' is not possible. Check reach connectivity.'
4985 cprior = this%divcprior(jpos)
4988 if (this%iprpak /= 0)
then
4989 call this%inputtab%add_term(n)
4990 call this%inputtab%add_term(idiv)
4991 call this%inputtab%add_term(nn)
4992 call this%inputtab%add_term(cprior)
5012 character(len=LINELENGTH) :: title
5013 character(len=LINELENGTH) :: text
5014 logical(LGP) :: lcycle
5015 logical(LGP) :: ladd
5016 character(len=5) :: crch
5017 character(len=5) :: crch2
5018 character(len=10) :: cval
5019 integer(I4B) :: maxcols
5020 integer(I4B) :: npairs
5021 integer(I4B) :: ipair
5025 integer(I4B) :: idiv
5028 integer(I4B) :: jpos
5034 if (this%iprpak /= 0)
then
5038 do n = 1, this%maxbound
5040 ec:
do i = this%ia(n) + 1, this%ia(n + 1) - 1
5043 if (this%idir(i) > 0) cycle ec
5047 if (this%iboundpak(n2) == 0) cycle ec
5051 npairs = max(npairs, ipair)
5054 maxcols = 1 + npairs * 2
5057 title = trim(adjustl(this%text))//
' PACKAGE ('// &
5058 trim(adjustl(this%packName))//
') CONNECTED REACH UPSTREAM '// &
5060 call table_cr(this%inputtab, this%packName, title)
5061 call this%inputtab%table_df(this%maxbound, maxcols, this%iout)
5063 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
5065 write (cval,
'(i10)') i
5066 text =
'DOWNSTREAM REACH '//trim(adjustl(cval))
5067 call this%inputtab%initialize_column(text, 10, alignment=
tabcenter)
5068 text =
'FRACTION '//trim(adjustl(cval))
5069 call this%inputtab%initialize_column(text, 12, alignment=
tabcenter)
5074 do n = 1, this%maxbound
5075 do idiv = 1, this%ndiv(n)
5077 i1 = this%iadiv(n + 1) - 1
5079 do i = this%ia(n) + 1, this%ia(n + 1) - 1
5081 if (this%divreach(jpos) == n2)
then
5082 this%idiv(i) = jpos - i0 + 1
5092 do n = 1, this%maxbound
5096 do i = this%ia(n) + 1, this%ia(n + 1) - 1
5097 if (this%idir(i) < 0)
then
5103 do idiv = 1, this%ndiv(n)
5104 jpos = this%iadiv(n) + idiv - 1
5105 n2 = this%divreach(jpos)
5107 if (f /=
dzero)
then
5108 write (
errmsg,
'(a,2(1x,i0,1x,a),1x,a,g0,a,2(1x,a))') &
5109 'Reach', n,
'is connected to reach', n2,
'by a diversion', &
5110 'but the upstream fraction is not equal to zero (', f,
'). Check', &
5111 trim(this%packName),
'package diversion and package data.'
5115 write (
warnmsg,
'(a,3(1x,a))') &
5117 'A warning instead of an error is issued because', &
5118 'the reach is only connected to the diversion reach in the ', &
5119 'downstream direction.'
5129 do n = 1, this%maxbound
5133 write (crch,
'(i5)') n
5134 if (this%iprpak /= 0)
then
5135 call this%inputtab%add_term(n)
5138 eachconn:
do i = this%ia(n) + 1, this%ia(n + 1) - 1
5142 this%qconn(i) =
dzero
5145 if (this%idir(i) > 0)
then
5151 if (this%iboundpak(n2) == 0)
then
5158 write (crch2,
'(i5)') n2
5161 f = f + this%ustrf(n2)
5162 write (cval,
'(f10.4)') this%ustrf(n2)
5165 if (this%iprpak /= 0)
then
5166 call this%inputtab%add_term(n2)
5167 call this%inputtab%add_term(this%ustrf(n2))
5169 eachdiv:
do idiv = 1, this%ndiv(n)
5170 jpos = this%iadiv(n) + idiv - 1
5171 if (this%divreach(jpos) == n2)
then
5177 rval = rval + this%ustrf(n2)
5180 this%ftotnd(n) = rval
5183 if (this%iprpak /= 0)
then
5185 do i = ipair, npairs
5186 call this%inputtab%add_term(
' ')
5187 call this%inputtab%add_term(
' ')
5195 write (
errmsg,
'(a,1x,i0,1x,a,g0,a,3(1x,a))') &
5196 'Upstream fractions for reach ', n,
'is not equal to one (', f, &
5197 '). Check', trim(this%packName),
'package reach connectivity and', &
5219 integer(I4B) :: nbudterm
5224 integer(I4B) :: maxlist
5225 integer(I4B) :: naux
5228 character(len=LENBUDTXT) :: text
5229 character(len=LENBUDTXT),
dimension(1) :: auxtxt
5236 if (this%imover == 1) nbudterm = nbudterm + 2
5237 if (this%naux > 0) nbudterm = nbudterm + 1
5241 call this%budobj%budgetobject_df(this%maxbound, nbudterm, 0, 0, &
5242 ibudcsv=this%ibudcsv)
5246 text =
' FLOW-JA-FACE'
5248 maxlist = this%nconn
5250 auxtxt(1) =
' FLOW-AREA'
5251 call this%budobj%budterm(idx)%initialize(text, &
5256 maxlist, .false., .false., &
5260 call this%budobj%budterm(idx)%reset(this%nconn)
5262 do n = 1, this%maxbound
5264 do i = this%ia(n) + 1, this%ia(n + 1) - 1
5266 call this%budobj%budterm(idx)%update_term(n1, n2, q)
5273 maxlist = this%maxbound - this%ianynone
5275 auxtxt(1) =
' FLOW-AREA'
5276 call this%budobj%budterm(idx)%initialize(text, &
5281 maxlist, .false., .true., &
5283 call this%budobj%budterm(idx)%reset(this%maxbound)
5285 do n = 1, this%maxbound
5286 n2 = this%igwfnode(n)
5288 call this%budobj%budterm(idx)%update_term(n, n2, q)
5295 maxlist = this%maxbound
5297 call this%budobj%budterm(idx)%initialize(text, &
5302 maxlist, .false., .false., &
5306 text =
' EVAPORATION'
5308 maxlist = this%maxbound
5310 call this%budobj%budterm(idx)%initialize(text, &
5315 maxlist, .false., .false., &
5321 maxlist = this%maxbound
5323 call this%budobj%budterm(idx)%initialize(text, &
5328 maxlist, .false., .false., &
5332 text =
' EXT-INFLOW'
5334 maxlist = this%maxbound
5336 call this%budobj%budterm(idx)%initialize(text, &
5341 maxlist, .false., .false., &
5345 text =
' EXT-OUTFLOW'
5347 maxlist = this%maxbound
5349 call this%budobj%budterm(idx)%initialize(text, &
5354 maxlist, .false., .false., &
5360 maxlist = this%maxbound
5362 auxtxt(1) =
' VOLUME'
5363 call this%budobj%budterm(idx)%initialize(text, &
5368 maxlist, .false., .false., &
5372 if (this%imover == 1)
then
5377 maxlist = this%maxbound
5379 call this%budobj%budterm(idx)%initialize(text, &
5384 maxlist, .false., .false., &
5390 maxlist = this%maxbound
5392 call this%budobj%budterm(idx)%initialize(text, &
5397 maxlist, .false., .false., &
5408 maxlist = this%maxbound
5409 call this%budobj%budterm(idx)%initialize(text, &
5414 maxlist, .false., .false., &
5419 if (this%iprflow /= 0)
then
5420 call this%budobj%flowtable_df(this%iout, cellids=
'GWF')
5438 integer(I4B) :: naux
5445 integer(I4B) :: idiv
5446 integer(I4B) :: jpos
5460 call this%budobj%budterm(idx)%reset(this%nconn)
5461 do n = 1, this%maxbound
5465 do i = this%ia(n) + 1, this%ia(n + 1) - 1
5467 if (this%iboundpak(n) /= 0)
then
5469 if (this%idir(i) < 0)
then
5475 do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
5476 if (this%idir(ii) > 0) cycle
5477 if (this%ja(ii) /= n) cycle
5483 call this%sfr_calc_reach_depth(n, qt, d)
5484 ca = this%calc_area_wet(n, d)
5489 this%qauxcbc(1) = ca
5490 call this%budobj%budterm(idx)%update_term(n1, n2, q, this%qauxcbc)
5496 call this%budobj%budterm(idx)%reset(this%maxbound - this%ianynone)
5497 do n = 1, this%maxbound
5498 n2 = this%igwfnode(n)
5500 if (this%iboundpak(n) /= 0)
then
5502 if (this%depth(n) >
dzero)
then
5503 wp = this%calc_perimeter_wet(n, this%depth(n))
5512 this%qauxcbc(1) =
dzero
5515 call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
5521 call this%budobj%budterm(idx)%reset(this%maxbound)
5522 do n = 1, this%maxbound
5523 if (this%iboundpak(n) /= 0)
then
5524 a = this%calc_surface_area(n)
5525 q = this%rain(n) * a
5529 call this%budobj%budterm(idx)%update_term(n, n, q)
5534 call this%budobj%budterm(idx)%reset(this%maxbound)
5535 do n = 1, this%maxbound
5536 if (this%iboundpak(n) /= 0)
then
5537 q = -this%simevap(n)
5541 call this%budobj%budterm(idx)%update_term(n, n, q)
5546 call this%budobj%budterm(idx)%reset(this%maxbound)
5547 do n = 1, this%maxbound
5548 if (this%iboundpak(n) /= 0)
then
5549 q = this%simrunoff(n)
5553 call this%budobj%budterm(idx)%update_term(n, n, q)
5558 call this%budobj%budterm(idx)%reset(this%maxbound)
5559 do n = 1, this%maxbound
5560 if (this%iboundpak(n) /= 0)
then
5565 call this%budobj%budterm(idx)%update_term(n, n, q)
5570 call this%budobj%budterm(idx)%reset(this%maxbound)
5571 do n = 1, this%maxbound
5573 if (this%iboundpak(n) /= 0)
then
5574 do i = this%ia(n) + 1, this%ia(n + 1) - 1
5575 if (this%idir(i) > 0) cycle
5578 jpos = this%iadiv(n) + idiv - 1
5579 q = q + this%divq(jpos)
5581 q = q + this%qconn(i)
5584 q = q - this%dsflow(n)
5585 if (this%imover == 1)
then
5586 q = q + this%pakmvrobj%get_qtomvr(n)
5589 if (this%imover == 1)
then
5590 q = this%pakmvrobj%get_qfrommvr(n)
5593 call this%budobj%budterm(idx)%update_term(n, n, q)
5598 call this%budobj%budterm(idx)%reset(this%maxbound)
5599 do n = 1, this%maxbound
5601 if (this%iboundpak(n) /= 0)
then
5603 a = this%calc_surface_area_wet(n, d)
5604 this%qauxcbc(1) = a * d
5607 this%qauxcbc(1) =
dzero
5609 call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc)
5613 if (this%imover == 1)
then
5617 call this%budobj%budterm(idx)%reset(this%maxbound)
5618 do n = 1, this%maxbound
5620 if (this%iboundpak(n) /= 0)
then
5621 q = this%pakmvrobj%get_qfrommvr(n)
5623 call this%budobj%budterm(idx)%update_term(n, n, q)
5628 call this%budobj%budterm(idx)%reset(this%maxbound)
5629 do n = 1, this%maxbound
5630 if (this%iboundpak(n) /= 0)
then
5631 q = this%pakmvrobj%get_qtomvr(n)
5638 call this%budobj%budterm(idx)%update_term(n, n, q)
5646 call this%budobj%budterm(idx)%reset(this%maxbound)
5647 do n = 1, this%maxbound
5649 call this%budobj%budterm(idx)%update_term(n, n, q, this%auxvar(:, n))
5654 call this%budobj%accumulate_terms()
5671 integer(I4B) :: nterms
5672 character(len=LINELENGTH) :: title
5673 character(len=LINELENGTH) :: text
5676 if (this%iprhed > 0)
then
5683 if (this%inamedbound == 1)
then
5688 title = trim(adjustl(this%text))//
' PACKAGE ('// &
5689 trim(adjustl(this%packName))//
') STAGES FOR EACH CONTROL VOLUME'
5692 call table_cr(this%stagetab, this%packName, title)
5693 call this%stagetab%table_df(this%maxbound, nterms, this%iout, &
5697 if (this%inamedbound == 1)
then
5699 call this%stagetab%initialize_column(text,
lenboundname, &
5705 call this%stagetab%initialize_column(text, 10, alignment=
tabcenter)
5709 call this%stagetab%initialize_column(text, 20, alignment=
tableft)
5713 call this%stagetab%initialize_column(text, 12, alignment=
tabcenter)
5717 call this%stagetab%initialize_column(text, 12, alignment=
tabcenter)
5721 call this%stagetab%initialize_column(text, 12, alignment=
tabcenter)
5725 call this%stagetab%initialize_column(text, 12, alignment=
tabcenter)
5728 text =
'STREAMBED CONDUCTANCE'
5729 call this%stagetab%initialize_column(text, 12, alignment=
tabcenter)
5732 text =
'STREAMBED GRADIENT'
5733 call this%stagetab%initialize_column(text, 12, alignment=
tabcenter)
5752 integer(I4B),
intent(in) :: n
5753 real(dp),
intent(in) :: depth
5755 integer(I4B) :: npts
5760 npts = this%ncrosspts(n)
5761 i0 = this%iacross(n)
5762 i1 = this%iacross(n + 1) - 1
5765 this%xsheight(i0:i1), depth)
5784 integer(I4B),
intent(in) :: n
5785 real(dp),
intent(in) :: depth
5787 integer(I4B) :: npts
5792 npts = this%ncrosspts(n)
5793 i0 = this%iacross(n)
5794 i1 = this%iacross(n + 1) - 1
5797 this%xsheight(i0:i1), depth)
5816 integer(I4B),
intent(in) :: n
5818 integer(I4B) :: npts
5821 real(dp) :: top_width
5824 npts = this%ncrosspts(n)
5825 i0 = this%iacross(n)
5826 i1 = this%iacross(n + 1) - 1
5830 top_width = this%station(i0)
5848 integer(I4B),
intent(in) :: n
5849 real(dp),
intent(in) :: depth
5851 real(dp) :: top_width
5854 top_width = this%calc_top_width_wet(n, depth)
5871 integer(I4B),
intent(in) :: n
5872 real(dp),
intent(in) :: depth
5874 integer(I4B) :: npts
5880 npts = this%ncrosspts(n)
5881 i0 = this%iacross(n)
5882 i1 = this%iacross(n + 1) - 1
5886 this%station(i0:i1), &
5887 this%xsheight(i0:i1), &
5906 class(
sfrtype),
intent(inout) :: this
5913 call mem_reallocate(this%denseterms, 3, this%MAXBOUND,
'DENSETERMS', &
5915 do i = 1, this%maxbound
5917 this%denseterms(j, i) =
dzero
5920 write (this%iout,
'(/1x,a)')
'DENSITY TERMS HAVE BEEN ACTIVATED FOR SFR &
5921 &PACKAGE: '//trim(adjustl(this%packName))
5937 class(
sfrtype),
intent(inout) :: this
5944 call mem_reallocate(this%viscratios, 2, this%MAXBOUND,
'VISCRATIOS', &
5946 do i = 1, this%maxbound
5948 this%viscratios(j, i) =
done
5951 write (this%iout,
'(/1x,a)')
'VISCOSITY HAS BEEN ACTIVATED FOR SFR &
5952 &PACKAGE: '//trim(adjustl(this%packName))
5971 bots, flow, gwfhcof, gwfrhs)
5973 class(
sfrtype),
intent(inout) :: this
5974 integer(I4B),
intent(in) :: n
5975 real(DP),
intent(in) :: stage
5976 real(DP),
intent(in) :: head
5977 real(DP),
intent(in) :: cond
5978 real(DP),
intent(in) :: bots
5979 real(DP),
intent(inout) :: flow
5980 real(DP),
intent(inout) :: gwfhcof
5981 real(DP),
intent(inout) :: gwfrhs
5986 real(DP) :: rdensesfr
5987 real(DP) :: rdensegwf
5988 real(DP) :: rdenseavg
5994 logical(LGP) :: stage_below_bot
5995 logical(LGP) :: head_below_bot
5998 if (stage >= bots)
then
6000 stage_below_bot = .false.
6001 rdensesfr = this%denseterms(1, n)
6004 stage_below_bot = .true.
6005 rdensesfr = this%denseterms(2, n)
6009 if (head >= bots)
then
6011 head_below_bot = .false.
6012 rdensegwf = this%denseterms(2, n)
6015 head_below_bot = .true.
6016 rdensegwf = this%denseterms(1, n)
6020 if (rdensegwf ==
dzero)
return
6023 if (stage_below_bot .and. head_below_bot)
then
6030 rdenseavg =
dhalf * (rdensesfr + rdensegwf)
6034 d1 = cond * (rdenseavg -
done)
6035 gwfhcof = gwfhcof - d1
6036 gwfrhs = gwfrhs - d1 * ss
6041 if (.not. stage_below_bot .and. .not. head_below_bot)
then
6045 elevgwf = this%denseterms(3, n)
6047 elevavg =
dhalf * (elevsfr + elevgwf)
6048 havg =
dhalf * (hh + ss)
6049 d2 = cond * (havg - elevavg) * (rdensegwf - rdensesfr)
6050 gwfrhs = gwfrhs + d2
This module contains the base boundary package.
This module contains the BudgetModule.
subroutine, public budgetobject_cr(this, name)
Create a new budget object.
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
real(dp), parameter dhdry
real dry cell constant
@ tabcenter
centered table column
@ tabright
right justified table column
@ tableft
left justified table column
@ mnormal
normal output mode
real(dp), parameter dtwothirds
real constant 2/3
integer(i4b), parameter lenpackagename
maximum length of the package name
real(dp), parameter dp9
real constant 9/10
real(dp), parameter deight
real constant 8
real(dp), parameter dfivethirds
real constant 5/3
real(dp), parameter dp999
real constant 999/1000
integer(i4b), parameter namedboundflag
named bound flag
real(dp), parameter donethird
real constant 1/3
real(dp), parameter dnodata
real no data constant
real(dp), parameter d1p1
real constant 1.1
real(dp), parameter dhnoflo
real no flow constant
real(dp), parameter dhundred
real constant 100
integer(i4b), parameter lenpakloc
maximum length of a package location
integer(i4b), parameter lentimeseriesname
maximum length of a time series name
real(dp), parameter dep20
real constant 1e20
real(dp), parameter dp6
real constant 3/5
integer(i4b), parameter maxadpit
maximum advanced package Newton-Raphson iterations
real(dp), parameter dhalf
real constant 1/2
integer(i4b), parameter lenftype
maximum length of a package type (DIS, WEL, OC, etc.)
real(dp), parameter dpi
real constant
real(dp), parameter dp99
real constant 99/100
integer(i4b), parameter lenboundname
maximum length of a bound name
real(dp), parameter dem4
real constant 1e-4
real(dp), parameter dem30
real constant 1e-30
real(dp), parameter dem6
real constant 1e-6
real(dp), parameter dzero
real constant zero
real(dp), parameter dem5
real constant 1e-5
real(dp), parameter dprec
real constant machine precision
integer(i4b), parameter maxcharlen
maximum length of char string
real(dp), parameter dp7
real constant 7/10
real(dp), parameter dem2
real constant 1e-2
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 stateless sfr subroutines and functions.
real(dp) function, public get_wetted_topwidth(npts, stations, heights, d)
Calculate the wetted top width for a reach.
real(dp) function, public get_wetted_perimeter(npts, stations, heights, d)
Calculate the wetted perimeter for a reach.
real(dp) function, public get_cross_section_area(npts, stations, heights, d)
Calculate the cross-sectional area for a reach.
real(dp) function, public get_saturated_topwidth(npts, stations)
Calculate the saturated top width for a reach.
real(dp) function, public get_mannings_section(npts, stations, heights, roughfracs, roughness, conv_fact, slope, d)
Calculate the manning's discharge for a reach.
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 derived types ObserveType and ObsDataType.
subroutine, public cross_section_cr(this, iout, iprpak, nreaches)
Create a cross-section object.
This module contains the SFR package methods.
real(dp) function calc_top_width_wet(this, n, depth)
Calculate wetted top width.
subroutine sfr_setup_tableobj(this)
Setup stage table object for package.
subroutine sfr_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
@ brief Convergence check for package.
subroutine sfr_cq(this, x, flowja, iadv)
@ brief Calculate package flows.
subroutine sfr_ot_package_flows(this, icbcfl, ibudfl)
@ brief Output package flow terms.
subroutine sfr_activate_viscosity(this)
Activate viscosity terms.
subroutine sfr_calc_constant(this, n, d1, hgwf, qgwf, qd)
subroutine sfr_da(this)
@ brief Deallocate package memory
subroutine sfr_read_diversions(this)
@ brief Read diversions for the package
subroutine sfr_calc_reach_depth(this, n, q1, d1)
Calculate the depth at the midpoint.
real(dp) function calc_surface_area(this, n)
Calculate maximum surface area.
subroutine sfr_check_ustrf(this)
Check upstream fraction data.
subroutine sfr_read_connectiondata(this)
@ brief Read connectiondata for the package
subroutine sfr_calc_xs_depth(this, n, qrch, d)
Calculate the depth at the midpoint of a irregular cross-section.
subroutine sfr_set_stressperiod(this, n, ichkustrm, crossfile)
Set period data.
subroutine sfr_check_diversions(this)
Check diversions data.
subroutine sfr_ot_dv(this, idvsave, idvprint)
@ brief Output package dependent-variable terms.
subroutine sfr_calculate_density_exchange(this, n, stage, head, cond, bots, flow, gwfhcof, gwfrhs)
Calculate density terms.
subroutine sfr_fc(this, rhs, ia, idxglo, matrix_sln)
@ brief Copy hcof and rhs terms into solution.
subroutine, public sfr_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
@ brief Create a new package object
subroutine sfr_fn(this, rhs, ia, idxglo, matrix_sln)
@ brief Add Newton-Raphson terms for package into solution.
subroutine sfr_process_obsid(obsrv, dis, inunitobs, iout)
Process observation IDs for a package.
subroutine define_listlabel(this)
@ brief Define the list label for the package
subroutine sfr_df_obs(this)
Define the observation types available in the package.
subroutine sfr_options(this, option, found)
@ brief Read additional options for package
subroutine sfr_calc_cond(this, n, depth, cond, hsfr, h_temp)
Calculate reach-aquifer conductance.
subroutine sfr_check_connections(this)
Check connection data.
subroutine sfr_allocate_arrays(this)
@ brief Allocate arrays
subroutine sfr_bd_obs(this)
Save observations for the package.
subroutine sfr_setup_budobj(this)
Setup budget object for package.
subroutine sfr_rp(this)
@ brief Read and prepare period data for package
subroutine sfr_ar(this)
@ brief Allocate and read method for package
subroutine sfr_calc_qman(this, n, depth, qman)
Calculate streamflow.
subroutine sfr_check_reaches(this)
Check reach data.
real(dp) function calc_perimeter_wet(this, n, depth)
Calculate wetted perimeter.
subroutine sfr_calc_steady(this, n, d1, hgwf, qu, qi, qfrommvr, qr, qe, qro, qgwf, qd)
subroutine sfr_read_packagedata(this)
@ brief Read packagedata for the package
subroutine sfr_update_flows(this, n, qd, qgwf)
Update flow terms.
subroutine sfr_calc_qd(this, n, depth, hgwf, qgwf, qd)
Calculate downstream flow term.
subroutine sfr_fill_budobj(this)
Copy flow terms into budget object for package.
subroutine sfr_cf(this)
@ brief Formulate the package hcof and rhs terms.
subroutine sfr_ot_bdsummary(this, kstp, kper, iout, ibudfl)
@ brief Output advanced package budget summary.
subroutine sfr_calc_div(this, n, i, qd, qdiv)
Calculate diversion flow.
subroutine sfr_calc_qgwf(this, n, depth, hgwf, qgwf, gwfhcof, gwfrhs)
Calculate reach-aquifer exchange.
character(len=lenftype) ftype
package ftype string
character(len=lenpackagename) text
package budget string
subroutine sfr_calc_qsource(this, n, depth, qsrc)
Calculate sum of sources.
subroutine sfr_allocate_scalars(this)
@ brief Allocate scalars
logical function sfr_obs_supported(this)
Determine if observations are supported.
subroutine sfr_solve(this, n, h, hcof, rhs, update)
Solve reach continuity equation.
subroutine sfr_read_dimensions(this)
@ brief Read dimensions for package
real(dp) function calc_area_wet(this, n, depth)
Calculate wetted area.
subroutine sfr_read_crossection(this)
@ brief Read crosssection block for the package
subroutine sfr_rp_obs(this)
Read and prepare observations for a package.
subroutine sfr_activate_density(this)
Activate density terms.
subroutine sfr_adjust_ro_ev(this, qc, qu, qi, qr, qro, qe, qfrommvr)
Adjust runoff and evaporation.
subroutine sfr_check_conversion(this)
Check unit conversion data.
real(dp) function calc_surface_area_wet(this, n, depth)
Calculate wetted surface area.
integer(i4b) function sfr_gwf_conn(this, n)
Determine if a reach is connected to a gwf cell.
subroutine sfr_ad(this)
@ brief Advance the package
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.
subroutine, public deprecation_warning(cblock, cvar, cver, endmsg, iunit)
Store deprecation warning message.
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
This module contains simulation variables.
character(len=maxcharlen) errmsg
error message string
character(len=maxcharlen) warnmsg
warning message string
real(dp) function squadraticsaturation(top, bot, x, eps)
@ brief sQuadraticSaturation
real(dp) function scubicsaturation(top, bot, x, eps)
@ brief sCubicSaturation
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function
real(dp) function sqsaturationderivative(top, bot, x, c1, c2)
@ brief sQSaturationDerivative
subroutine schsmooth(d, smooth, dwdh)
@ brief sChSmooth
real(dp) function sqsaturation(top, bot, x, c1, c2)
@ brief sQSaturation
subroutine, public table_cr(this, name, title)
real(dp), pointer, public pertim
time relative to start of stress period
real(dp), pointer, public totim
time relative to start of simulation
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
subroutine, public read_value_or_time_series_adv(textInput, ii, jj, bndElem, pkgName, auxOrBnd, tsManager, iprpak, varName)
Call this subroutine from advanced packages to define timeseries link for a variable (varName).
logical function, public var_timeseries(tsManager, pkgName, varName, auxOrBnd)
Determine if a timeseries link with varName is defined.
Derived type for the Budget object.