41 character(len=LENFTYPE) ::
ftype =
'LAK'
42 character(len=LENPACKAGENAME) ::
text =
' LAK'
45 real(dp),
dimension(:),
pointer,
contiguous :: tabstage => null()
46 real(dp),
dimension(:),
pointer,
contiguous :: tabvolume => null()
47 real(dp),
dimension(:),
pointer,
contiguous :: tabsarea => null()
48 real(dp),
dimension(:),
pointer,
contiguous :: tabwarea => null()
54 character(len=16),
dimension(:),
pointer,
contiguous :: clakbudget => null()
55 character(len=16),
dimension(:),
pointer,
contiguous :: cauxcbc => null()
58 integer(I4B),
pointer :: iprhed => null()
59 integer(I4B),
pointer :: istageout => null()
60 integer(I4B),
pointer :: ibudgetout => null()
61 integer(I4B),
pointer :: ibudcsv => null()
62 integer(I4B),
pointer :: ipakcsv => null()
63 character(len=:),
allocatable :: pakcsvfile
64 integer(I4B),
pointer :: cbcauxitems => null()
65 integer(I4B),
pointer :: nlakes => null()
66 integer(I4B),
pointer :: noutlets => null()
67 integer(I4B),
pointer :: ntables => null()
68 real(dp),
pointer :: convlength => null()
69 real(dp),
pointer :: convtime => null()
70 real(dp),
pointer :: outdmax => null()
71 integer(I4B),
pointer :: igwhcopt => null()
72 integer(I4B),
pointer :: iconvchk => null()
73 integer(I4B),
pointer :: maxlakit => null()
74 real(dp),
pointer :: surfdep => null()
75 real(dp),
pointer :: dmaxchg => null()
76 real(dp),
pointer :: delh => null()
77 integer(I4B),
pointer :: check_attr => null()
80 integer(I4B),
pointer :: iimplicit => null()
81 integer(I4B),
pointer :: iforceleg => null()
82 integer(I4B),
pointer :: iforceleglak => null()
84 integer(I4B),
pointer :: bditems => null()
87 integer(I4B),
dimension(:),
pointer,
contiguous :: nlakeconn => null()
88 integer(I4B),
dimension(:),
pointer,
contiguous :: idxlakeconn => null()
89 integer(I4B),
dimension(:),
pointer,
contiguous :: ntabrow => null()
90 real(dp),
dimension(:),
pointer,
contiguous :: strt => null()
91 real(dp),
dimension(:),
pointer,
contiguous :: laketop => null()
92 real(dp),
dimension(:),
pointer,
contiguous :: lakebot => null()
93 real(dp),
dimension(:),
pointer,
contiguous :: sareamax => null()
94 character(len=LENBOUNDNAME),
dimension(:),
pointer, &
95 contiguous :: lakename => null()
96 character(len=8),
dimension(:),
pointer,
contiguous :: status => null()
97 real(dp),
dimension(:),
pointer,
contiguous :: avail => null()
98 real(dp),
dimension(:),
pointer,
contiguous :: lkgwsink => null()
99 real(dp),
dimension(:),
pointer,
contiguous :: stage => null()
100 real(dp),
dimension(:),
pointer,
contiguous :: rainfall => null()
101 real(dp),
dimension(:),
pointer,
contiguous :: evaporation => null()
102 real(dp),
dimension(:),
pointer,
contiguous :: runoff => null()
103 real(dp),
dimension(:),
pointer,
contiguous :: inflow => null()
104 real(dp),
dimension(:),
pointer,
contiguous :: withdrawal => null()
105 real(dp),
dimension(:, :),
pointer,
contiguous :: lauxvar => null()
108 integer(I4B),
dimension(:),
pointer,
contiguous :: ialaktab => null()
109 real(dp),
dimension(:),
pointer,
contiguous :: tabstage => null()
110 real(dp),
dimension(:),
pointer,
contiguous :: tabvolume => null()
111 real(dp),
dimension(:),
pointer,
contiguous :: tabsarea => null()
112 real(dp),
dimension(:),
pointer,
contiguous :: tabwarea => null()
115 integer(I4B),
dimension(:),
pointer,
contiguous :: ncncvr => null()
119 integer(I4B),
dimension(:),
pointer,
contiguous :: ilegacy => null()
120 integer(I4B),
dimension(:),
pointer,
contiguous :: nstuck => null()
121 real(dp),
dimension(:),
pointer,
contiguous :: surfin => null()
122 real(dp),
dimension(:),
pointer,
contiguous :: surfout => null()
123 real(dp),
dimension(:),
pointer,
contiguous :: surfout1 => null()
124 real(dp),
dimension(:),
pointer,
contiguous :: precip => null()
125 real(dp),
dimension(:),
pointer,
contiguous :: precip1 => null()
126 real(dp),
dimension(:),
pointer,
contiguous :: evap => null()
127 real(dp),
dimension(:),
pointer,
contiguous :: evap1 => null()
128 real(dp),
dimension(:),
pointer,
contiguous :: evapo => null()
129 real(dp),
dimension(:),
pointer,
contiguous :: withr => null()
130 real(dp),
dimension(:),
pointer,
contiguous :: withr1 => null()
131 real(dp),
dimension(:),
pointer,
contiguous :: flwin => null()
132 real(dp),
dimension(:),
pointer,
contiguous :: flwiter => null()
133 real(dp),
dimension(:),
pointer,
contiguous :: flwiter1 => null()
134 real(dp),
dimension(:),
pointer,
contiguous :: seep => null()
135 real(dp),
dimension(:),
pointer,
contiguous :: seep1 => null()
136 real(dp),
dimension(:),
pointer,
contiguous :: seep0 => null()
137 real(dp),
dimension(:),
pointer,
contiguous :: stageiter => null()
138 real(dp),
dimension(:),
pointer,
contiguous :: chterm => null()
141 integer(I4B),
dimension(:),
pointer,
contiguous :: iseepc => null()
142 integer(I4B),
dimension(:),
pointer,
contiguous :: idhc => null()
143 real(dp),
dimension(:),
pointer,
contiguous :: en1 => null()
144 real(dp),
dimension(:),
pointer,
contiguous :: en2 => null()
145 real(dp),
dimension(:),
pointer,
contiguous :: r1 => null()
146 real(dp),
dimension(:),
pointer,
contiguous :: r2 => null()
147 real(dp),
dimension(:),
pointer,
contiguous :: dh0 => null()
148 real(dp),
dimension(:),
pointer,
contiguous :: s0 => null()
149 real(dp),
dimension(:),
pointer,
contiguous :: qgwf0 => null()
151 integer(I4B),
dimension(:),
pointer,
contiguous :: idxlocnode => null()
152 integer(I4B),
dimension(:),
pointer,
contiguous :: idxdiag => null()
153 integer(I4B),
dimension(:),
pointer,
contiguous :: idxoffdglo => null()
154 integer(I4B),
dimension(:),
pointer,
contiguous :: idxsymdglo => null()
155 integer(I4B),
dimension(:),
pointer,
contiguous :: idxsymoffdglo => null()
158 integer(I4B),
dimension(:),
pointer,
contiguous :: imap => null()
159 integer(I4B),
dimension(:),
pointer,
contiguous :: cellid => null()
160 integer(I4B),
dimension(:),
pointer,
contiguous :: nodesontop => null()
161 integer(I4B),
dimension(:),
pointer,
contiguous :: ictype => null()
162 real(dp),
dimension(:),
pointer,
contiguous :: bedleak => null()
163 real(dp),
dimension(:),
pointer,
contiguous :: belev => null()
164 real(dp),
dimension(:),
pointer,
contiguous :: telev => null()
165 real(dp),
dimension(:),
pointer,
contiguous :: connlength => null()
166 real(dp),
dimension(:),
pointer,
contiguous :: connwidth => null()
167 real(dp),
dimension(:),
pointer,
contiguous :: sarea => null()
168 real(dp),
dimension(:),
pointer,
contiguous :: warea => null()
169 real(dp),
dimension(:),
pointer,
contiguous :: satcond => null()
170 real(dp),
dimension(:),
pointer,
contiguous :: simcond => null()
171 real(dp),
dimension(:),
pointer,
contiguous :: simlakgw => null()
174 integer(I4B),
dimension(:),
pointer,
contiguous :: lakein => null()
175 integer(I4B),
dimension(:),
pointer,
contiguous :: lakeout => null()
176 integer(I4B),
dimension(:),
pointer,
contiguous :: iouttype => null()
177 real(dp),
dimension(:),
pointer,
contiguous :: outrate => null()
178 real(dp),
dimension(:),
pointer,
contiguous :: outinvert => null()
179 real(dp),
dimension(:),
pointer,
contiguous :: outwidth => null()
180 real(dp),
dimension(:),
pointer,
contiguous :: outrough => null()
181 real(dp),
dimension(:),
pointer,
contiguous :: outslope => null()
182 real(dp),
dimension(:),
pointer,
contiguous :: simoutrate => null()
185 real(dp),
dimension(:),
pointer,
contiguous :: qauxcbc => null()
186 real(dp),
dimension(:),
pointer,
contiguous :: dbuff => null()
187 real(dp),
dimension(:),
pointer,
contiguous :: qleak => null()
191 real(dp),
dimension(:),
pointer,
contiguous :: holdconn => null()
192 real(dp),
dimension(:),
pointer,
contiguous :: qsto => null()
195 integer(I4B),
pointer :: gwfiss => null()
196 real(dp),
dimension(:),
pointer,
contiguous :: gwfk11 => null()
197 real(dp),
dimension(:),
pointer,
contiguous :: gwfk33 => null()
198 real(dp),
dimension(:),
pointer,
contiguous :: gwfsat => null()
199 integer(I4B),
pointer :: gwfik33 => null()
202 integer(I4B),
dimension(:),
pointer,
contiguous :: iboundpak => null()
203 real(dp),
dimension(:),
pointer,
contiguous :: xnewpak => null()
204 real(dp),
dimension(:),
pointer,
contiguous :: xoldpak => null()
214 integer(I4B),
pointer :: idense
215 real(dp),
dimension(:, :),
pointer,
contiguous :: denseterms => null()
218 real(dp),
dimension(:, :),
pointer,
contiguous :: viscratios => null()
240 procedure,
private :: lak_set_legacy
241 procedure,
private :: lak_check_disconnected
295 procedure,
private :: lak_fc_implicit
296 procedure,
private :: lak_budget_nogwf
317 module subroutine lak_budget_nogwf(this, n, stage, b)
318 class(laktype),
intent(inout) :: this
319 integer(I4B),
intent(in) :: n
320 real(dp),
intent(in) :: stage
321 real(dp),
intent(inout) :: b
326 module subroutine lak_fc_implicit(this, rhs, matrix_sln)
327 class(laktype) :: this
328 real(dp),
dimension(:),
intent(inout) :: rhs
334 module subroutine lak_set_legacy(this, kiter, icnvgmod)
335 class(laktype),
intent(inout) :: this
336 integer(I4B),
intent(in) :: kiter
337 integer(I4B),
intent(in) :: icnvgmod
342 module subroutine lak_check_disconnected(this)
343 class(laktype),
intent(inout) :: this
351 subroutine lak_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
353 class(
bndtype),
pointer :: packobj
354 integer(I4B),
intent(in) :: id
355 integer(I4B),
intent(in) :: ibcnum
356 integer(I4B),
intent(in) :: inunit
357 integer(I4B),
intent(in) :: iout
358 character(len=*),
intent(in) :: namemodel
359 character(len=*),
intent(in) :: pakname
361 type(laktype),
pointer :: lakobj
368 call packobj%set_names(ibcnum, namemodel, pakname, ftype)
372 call lakobj%lak_allocate_scalars()
375 call packobj%pack_initialize()
377 packobj%inunit = inunit
380 packobj%ibcnum = ibcnum
385 end subroutine lak_create
389 subroutine lak_allocate_scalars(this)
391 class(laktype),
intent(inout) :: this
394 call this%BndType%allocate_scalars()
397 call mem_allocate(this%iprhed,
'IPRHED', this%memoryPath)
398 call mem_allocate(this%istageout,
'ISTAGEOUT', this%memoryPath)
399 call mem_allocate(this%ibudgetout,
'IBUDGETOUT', this%memoryPath)
400 call mem_allocate(this%ibudcsv,
'IBUDCSV', this%memoryPath)
401 call mem_allocate(this%ipakcsv,
'IPAKCSV', this%memoryPath)
402 call mem_allocate(this%nlakes,
'NLAKES', this%memoryPath)
403 call mem_allocate(this%noutlets,
'NOUTLETS', this%memoryPath)
404 call mem_allocate(this%ntables,
'NTABLES', this%memoryPath)
405 call mem_allocate(this%convlength,
'CONVLENGTH', this%memoryPath)
406 call mem_allocate(this%convtime,
'CONVTIME', this%memoryPath)
407 call mem_allocate(this%outdmax,
'OUTDMAX', this%memoryPath)
408 call mem_allocate(this%igwhcopt,
'IGWHCOPT', this%memoryPath)
409 call mem_allocate(this%iconvchk,
'ICONVCHK', this%memoryPath)
410 call mem_allocate(this%maxlakit,
'MAXLAKIT', this%memoryPath)
411 call mem_allocate(this%surfdep,
'SURFDEP', this%memoryPath)
412 call mem_allocate(this%dmaxchg,
'DMAXCHG', this%memoryPath)
414 call mem_allocate(this%check_attr,
'CHECK_ATTR', this%memoryPath)
415 call mem_allocate(this%iimplicit,
'IIMPLICIT', this%memoryPath)
416 call mem_allocate(this%iforceleg,
'IFORCELEG', this%memoryPath)
417 call mem_allocate(this%iforceleglak,
'IFORCELEGLAK', this%memoryPath)
418 call mem_allocate(this%bditems,
'BDITEMS', this%memoryPath)
419 call mem_allocate(this%cbcauxitems,
'CBCAUXITEMS', this%memoryPath)
420 call mem_allocate(this%idense,
'IDENSE', this%memoryPath)
431 this%convlength =
done
439 this%delh =
dp999 * this%dmaxchg
442 this%iforceleglak = 0
447 end subroutine lak_allocate_scalars
451 subroutine lak_allocate_arrays(this)
454 class(laktype),
intent(inout) :: this
459 call this%BndType%allocate_arrays()
462 allocate (this%clakbudget(this%bditems))
465 this%clakbudget(1) =
' GWF'
466 this%clakbudget(2) =
' RAINFALL'
467 this%clakbudget(3) =
' EVAPORATION'
468 this%clakbudget(4) =
' RUNOFF'
469 this%clakbudget(5) =
' EXT-INFLOW'
470 this%clakbudget(6) =
' WITHDRAWAL'
471 this%clakbudget(7) =
' EXT-OUTFLOW'
472 this%clakbudget(8) =
' STORAGE'
473 this%clakbudget(9) =
' CONSTANT'
474 this%clakbudget(10) =
' FROM-MVR'
475 this%clakbudget(11) =
' TO-MVR'
478 if (this%istageout > 0)
then
479 call mem_allocate(this%dbuff, this%nlakes,
'DBUFF', this%memoryPath)
480 do i = 1, this%nlakes
481 this%dbuff(i) =
dzero
484 call mem_allocate(this%dbuff, 0,
'DBUFF', this%memoryPath)
488 allocate (this%cauxcbc(this%cbcauxitems))
491 call mem_allocate(this%qauxcbc, this%cbcauxitems,
'QAUXCBC', this%memoryPath)
492 do i = 1, this%cbcauxitems
493 this%qauxcbc(i) =
dzero
497 call mem_allocate(this%qleak, this%maxbound,
'QLEAK', this%memoryPath)
498 do i = 1, this%maxbound
499 this%qleak(i) =
dzero
503 if (this%iimplicit /= 0)
then
504 call mem_allocate(this%holdconn, this%maxbound,
'HOLDCONN', this%memoryPath)
505 do i = 1, this%maxbound
506 this%holdconn(i) =
dzero
509 call mem_allocate(this%holdconn, 0,
'HOLDCONN', this%memoryPath)
511 call mem_allocate(this%qsto, this%nlakes,
'QSTO', this%memoryPath)
512 do i = 1, this%nlakes
517 call mem_allocate(this%denseterms, 3, 0,
'DENSETERMS', this%memoryPath)
520 call mem_allocate(this%viscratios, 2, 0,
'VISCRATIOS', this%memoryPath)
521 end subroutine lak_allocate_arrays
525 subroutine lak_read_lakes(this)
531 class(laktype),
intent(inout) :: this
533 character(len=LINELENGTH) :: text
534 character(len=LENBOUNDNAME) :: bndName, bndNameTemp
535 character(len=9) :: cno
536 character(len=50),
dimension(:),
allocatable :: caux
537 integer(I4B) :: ierr, ival
538 logical(LGP) :: isfound, endOfBlock
540 integer(I4B) :: ii, jj
544 integer(I4B) :: nconn
545 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
546 real(DP),
pointer :: bndElem => null()
552 call mem_allocate(this%nlakeconn, this%nlakes,
'NLAKECONN', this%memoryPath)
553 call mem_allocate(this%idxlakeconn, this%nlakes + 1,
'IDXLAKECONN', &
555 call mem_allocate(this%ntabrow, this%nlakes,
'NTABROW', this%memoryPath)
556 call mem_allocate(this%strt, this%nlakes,
'STRT', this%memoryPath)
557 call mem_allocate(this%laketop, this%nlakes,
'LAKETOP', this%memoryPath)
558 call mem_allocate(this%lakebot, this%nlakes,
'LAKEBOT', this%memoryPath)
559 call mem_allocate(this%sareamax, this%nlakes,
'SAREAMAX', this%memoryPath)
560 call mem_allocate(this%stage, this%nlakes,
'STAGE', this%memoryPath)
561 call mem_allocate(this%rainfall, this%nlakes,
'RAINFALL', this%memoryPath)
562 call mem_allocate(this%evaporation, this%nlakes,
'EVAPORATION', &
564 call mem_allocate(this%runoff, this%nlakes,
'RUNOFF', this%memoryPath)
565 call mem_allocate(this%inflow, this%nlakes,
'INFLOW', this%memoryPath)
566 call mem_allocate(this%withdrawal, this%nlakes,
'WITHDRAWAL', this%memoryPath)
567 call mem_allocate(this%lauxvar, this%naux, this%nlakes,
'LAUXVAR', &
569 call mem_allocate(this%avail, this%nlakes,
'AVAIL', this%memoryPath)
570 call mem_allocate(this%lkgwsink, this%nlakes,
'LKGWSINK', this%memoryPath)
571 call mem_allocate(this%ncncvr, this%nlakes,
'NCNCVR', this%memoryPath)
572 call mem_allocate(this%ilegacy, this%nlakes,
'ILEGACY', this%memoryPath)
573 call mem_allocate(this%nstuck, this%nlakes,
'NSTUCK', this%memoryPath)
574 call mem_allocate(this%surfin, this%nlakes,
'SURFIN', this%memoryPath)
575 call mem_allocate(this%surfout, this%nlakes,
'SURFOUT', this%memoryPath)
576 call mem_allocate(this%surfout1, this%nlakes,
'SURFOUT1', this%memoryPath)
577 call mem_allocate(this%precip, this%nlakes,
'PRECIP', this%memoryPath)
578 call mem_allocate(this%precip1, this%nlakes,
'PRECIP1', this%memoryPath)
579 call mem_allocate(this%evap, this%nlakes,
'EVAP', this%memoryPath)
580 call mem_allocate(this%evap1, this%nlakes,
'EVAP1', this%memoryPath)
581 call mem_allocate(this%evapo, this%nlakes,
'EVAPO', this%memoryPath)
582 call mem_allocate(this%withr, this%nlakes,
'WITHR', this%memoryPath)
583 call mem_allocate(this%withr1, this%nlakes,
'WITHR1', this%memoryPath)
584 call mem_allocate(this%flwin, this%nlakes,
'FLWIN', this%memoryPath)
585 call mem_allocate(this%flwiter, this%nlakes,
'FLWITER', this%memoryPath)
586 call mem_allocate(this%flwiter1, this%nlakes,
'FLWITER1', this%memoryPath)
587 call mem_allocate(this%seep, this%nlakes,
'SEEP', this%memoryPath)
588 call mem_allocate(this%seep1, this%nlakes,
'SEEP1', this%memoryPath)
589 call mem_allocate(this%seep0, this%nlakes,
'SEEP0', this%memoryPath)
590 call mem_allocate(this%stageiter, this%nlakes,
'STAGEITER', this%memoryPath)
591 call mem_allocate(this%chterm, this%nlakes,
'CHTERM', this%memoryPath)
597 if (this%iimplicit == 0)
then
598 call mem_allocate(this%iboundpak, this%nlakes,
'IBOUND', this%memoryPath)
599 call mem_allocate(this%xnewpak, this%nlakes,
'XNEWPAK', this%memoryPath)
601 call mem_allocate(this%xoldpak, this%nlakes,
'XOLDPAK', this%memoryPath)
604 call mem_allocate(this%iseepc, this%nlakes,
'ISEEPC', this%memoryPath)
605 call mem_allocate(this%idhc, this%nlakes,
'IDHC', this%memoryPath)
606 call mem_allocate(this%en1, this%nlakes,
'EN1', this%memoryPath)
607 call mem_allocate(this%en2, this%nlakes,
'EN2', this%memoryPath)
608 call mem_allocate(this%r1, this%nlakes,
'R1', this%memoryPath)
609 call mem_allocate(this%r2, this%nlakes,
'R2', this%memoryPath)
610 call mem_allocate(this%dh0, this%nlakes,
'DH0', this%memoryPath)
611 call mem_allocate(this%s0, this%nlakes,
'S0', this%memoryPath)
612 call mem_allocate(this%qgwf0, this%nlakes,
'QGWF0', this%memoryPath)
615 allocate (this%lakename(this%nlakes))
616 allocate (this%status(this%nlakes))
618 do n = 1, this%nlakes
620 this%status(n) =
'ACTIVE'
621 this%laketop(n) = -dep20
622 this%lakebot(n) = dep20
623 this%sareamax(n) = dzero
628 if (this%iimplicit == 0)
then
629 this%iboundpak(n) = 1
630 this%xnewpak(n) = dep20
632 this%xoldpak(n) = dep20
635 this%rainfall(n) = dzero
636 this%evaporation(n) = dzero
637 this%runoff(n) = dzero
638 this%inflow(n) = dzero
639 this%withdrawal(n) = dzero
645 if (this%naux > 0)
then
646 allocate (caux(this%naux))
650 allocate (nboundchk(this%nlakes))
651 do n = 1, this%nlakes
657 call this%parser%GetBlock(
'PACKAGEDATA', isfound, ierr, &
658 supportopenclose=.true.)
662 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
667 call this%parser%GetNextLine(endofblock)
669 n = this%parser%GetInteger()
671 if (n < 1 .or. n > this%nlakes)
then
672 write (
errmsg,
'(a,1x,i0)')
'lakeno MUST BE > 0 and <= ', this%nlakes
678 nboundchk(n) = nboundchk(n) + 1
681 this%strt(n) = this%parser%GetDouble()
684 ival = this%parser%GetInteger()
687 write (
errmsg,
'(a,1x,i0)')
'nlakeconn MUST BE >= 0 for lake ', n
694 if (this%iimplicit /= 0 .and. ival == 0)
then
695 write (
errmsg,
'(a,1x,i0,1x,a)') &
696 'lake', n,
'has no connections; the IMPLICIT option requires &
697 &each lake to have at least one GWF connection.'
702 this%nlakeconn(n) = ival
705 do iaux = 1, this%naux
706 call this%parser%GetString(caux(iaux))
710 write (cno,
'(i9.9)') n
711 bndname =
'Lake'//cno
714 if (this%inamedbound /= 0)
then
715 call this%parser%GetStringCaps(bndnametemp)
716 if (bndnametemp /=
'')
then
717 bndname = bndnametemp
720 this%lakename(n) = bndname
727 bndelem => this%lauxvar(jj, ii)
729 this%packName,
'AUX', &
730 this%tsManager, this%iprpak, &
738 do n = 1, this%nlakes
739 if (nboundchk(n) == 0)
then
740 write (
errmsg,
'(a,1x,i0)')
'NO DATA SPECIFIED FOR LAKE', n
742 else if (nboundchk(n) > 1)
then
743 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
744 'DATA FOR LAKE', n,
'SPECIFIED', nboundchk(n),
'TIMES'
749 write (this%iout,
'(1x,a)')
'END OF '//trim(adjustl(this%text))// &
752 call store_error(
'REQUIRED PACKAGEDATA BLOCK NOT FOUND.')
757 call this%parser%StoreErrorUnit()
761 this%MAXBOUND = nconn
762 write (this%iout,
'(//4x,a,i7)')
'MAXBOUND = ', this%maxbound
765 this%idxlakeconn(1) = 1
766 do n = 1, this%nlakes
767 this%idxlakeconn(n + 1) = this%idxlakeconn(n) + this%nlakeconn(n)
771 if (this%naux > 0)
then
776 deallocate (nboundchk)
777 end subroutine lak_read_lakes
781 subroutine lak_read_lake_connections(this)
785 class(laktype),
intent(inout) :: this
787 character(len=LINELENGTH) :: keyword, cellid
788 integer(I4B) :: ierr, ival
789 logical(LGP) :: isfound, endOfBlock
790 logical(LGP) :: is_lake_bed
794 integer(I4B) :: ipos, ipos0
795 integer(I4B) :: icellid, icellid0
798 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
799 character(len=LENVARNAME) :: ctypenm
802 allocate (nboundchk(this%MAXBOUND))
803 do n = 1, this%MAXBOUND
808 call this%parser%GetBlock(
'CONNECTIONDATA', isfound, ierr, &
809 supportopenclose=.true.)
814 call mem_allocate(this%imap, this%MAXBOUND,
'IMAP', this%memoryPath)
815 call mem_allocate(this%cellid, this%MAXBOUND,
'CELLID', this%memoryPath)
816 call mem_allocate(this%nodesontop, this%MAXBOUND,
'NODESONTOP', &
818 call mem_allocate(this%ictype, this%MAXBOUND,
'ICTYPE', this%memoryPath)
819 call mem_allocate(this%bedleak, this%MAXBOUND,
'BEDLEAK', this%memoryPath)
820 call mem_allocate(this%belev, this%MAXBOUND,
'BELEV', this%memoryPath)
821 call mem_allocate(this%telev, this%MAXBOUND,
'TELEV', this%memoryPath)
822 call mem_allocate(this%connlength, this%MAXBOUND,
'CONNLENGTH', &
824 call mem_allocate(this%connwidth, this%MAXBOUND,
'CONNWIDTH', &
826 call mem_allocate(this%sarea, this%MAXBOUND,
'SAREA', this%memoryPath)
827 call mem_allocate(this%warea, this%MAXBOUND,
'WAREA', this%memoryPath)
828 call mem_allocate(this%satcond, this%MAXBOUND,
'SATCOND', this%memoryPath)
829 call mem_allocate(this%simcond, this%MAXBOUND,
'SIMCOND', this%memoryPath)
830 call mem_allocate(this%simlakgw, this%MAXBOUND,
'SIMLAKGW', this%memoryPath)
833 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
836 call this%parser%GetNextLine(endofblock)
838 n = this%parser%GetInteger()
840 if (n < 1 .or. n > this%nlakes)
then
841 write (
errmsg,
'(a,1x,i0)')
'lakeno MUST BE > 0 and <= ', this%nlakes
847 ival = this%parser%GetInteger()
848 if (ival < 1 .or. ival > this%nlakeconn(n))
then
849 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
850 'iconn FOR LAKE ', n,
'MUST BE > 1 and <= ', this%nlakeconn(n)
856 ipos = this%idxlakeconn(n) + ival - 1
863 nboundchk(ipos) = nboundchk(ipos) + 1
866 call this%parser%GetCellid(this%dis%ndim, cellid)
867 nn = this%dis%noder_from_cellid(cellid, &
868 this%parser%iuactive, this%iout)
872 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
873 'INVALID cellid FOR LAKE ', n,
'connection', j
878 this%cellid(ipos) = nn
879 this%nodesontop(ipos) = nn
882 call this%parser%GetStringCaps(keyword)
883 select case (keyword)
885 this%ictype(ipos) = 0
887 this%ictype(ipos) = 1
889 this%ictype(ipos) = 2
891 this%ictype(ipos) = 3
893 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,a,a)') &
894 'UNKNOWN ctype FOR LAKE ', n,
'connection', j, &
895 '(', trim(keyword),
')'
898 write (ctypenm,
'(a16)') keyword
902 call this%parser%GetStringCaps(keyword)
903 select case (keyword)
905 is_lake_bed = .false.
906 this%bedleak(ipos) = dnodata
909 write (
warnmsg,
'(2(a,1x,i0,1x),a,1pe8.1,a)') &
910 'BEDLEAK for connection', j,
'in lake', n,
'is specified to '// &
911 'be NONE. Lake connections where the lake-GWF connection '// &
912 'conductance is solely a function of aquifer properties '// &
913 'in the connected GWF cell should be specified with a '// &
914 'DNODATA (', dnodata,
') value.'
917 call deprecation_warning(
'CONNECTIONDATA',
'bedleak=NONE',
'6.4.3', &
918 warnmsg, this%parser%GetUnit())
920 read (keyword, *) rval
922 is_lake_bed = .false.
926 this%bedleak(ipos) = rval
929 if (is_lake_bed .and. this%bedleak(ipos) < dzero)
then
930 write (
errmsg,
'(a,1x,i0,1x,a)')
'bedleak FOR LAKE ', n,
'MUST BE >= 0'
935 this%belev(ipos) = this%parser%GetDouble()
938 this%telev(ipos) = this%parser%GetDouble()
941 rval = this%parser%GetDouble()
942 if (rval <= dzero)
then
943 if (this%ictype(ipos) == 1 .or. this%ictype(ipos) == 2 .or. &
944 this%ictype(ipos) == 3)
then
945 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,a,1x,a)') &
946 'connection length (connlen) FOR LAKE ', n, &
947 ', CONNECTION NO.', j,
', MUST BE > 0 FOR SPECIFIED ', &
948 'connection type (ctype)', ctypenm
954 this%connlength(ipos) = rval
957 rval = this%parser%GetDouble()
958 if (rval < dzero)
then
959 if (this%ictype(ipos) == 1)
then
960 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
961 'cell width (connwidth) FOR LAKE ', n, &
962 ' HORIZONTAL CONNECTION ', j,
'MUST BE >= 0'
968 this%connwidth(ipos) = rval
970 write (this%iout,
'(1x,a)') &
971 'END OF '//trim(adjustl(this%text))//
' CONNECTIONDATA'
973 call store_error(
'REQUIRED CONNECTIONDATA BLOCK NOT FOUND.')
978 call this%parser%StoreErrorUnit()
982 do n = 1, this%nlakes
984 do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
985 if (this%ictype(ipos) /= 2 .and. this%ictype(ipos) /= 3) cycle
988 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
989 'nlakeconn FOR LAKE', n,
'EMBEDDED CONNECTION', j,
' EXCEEDS 1.'
996 do n = 1, this%nlakes
997 ipos0 = this%idxlakeconn(n)
998 icellid0 = this%cellid(ipos0)
999 if (this%ictype(ipos0) /= 2 .and. this%ictype(ipos0) /= 3) cycle
1000 do nn = 1, this%nlakes
1003 do ipos = this%idxlakeconn(nn), this%idxlakeconn(nn + 1) - 1
1005 icellid = this%cellid(ipos)
1006 if (icellid == icellid0)
then
1007 if (this%ictype(ipos) == 0)
then
1008 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
1009 'EMBEDDED LAKE', n, &
1010 'CANNOT COINCIDE WITH VERTICAL CONNECTION', j, &
1020 do n = 1, this%nlakes
1022 do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1024 nn = this%cellid(ipos)
1025 top = this%dis%top(nn)
1026 bot = this%dis%bot(nn)
1028 if (this%ictype(ipos) == 0)
then
1029 this%telev(ipos) = top + this%surfdep
1030 this%belev(ipos) = top
1031 this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
1033 else if (this%ictype(ipos) == 1)
then
1034 if (this%belev(ipos) == this%telev(ipos))
then
1035 this%telev(ipos) = top
1036 this%belev(ipos) = bot
1038 if (this%belev(ipos) >= this%telev(ipos))
then
1039 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1040 'telev FOR LAKE ', n,
' HORIZONTAL CONNECTION ', j, &
1043 else if (this%belev(ipos) < bot)
then
1044 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
1045 'belev FOR LAKE ', n,
' HORIZONTAL CONNECTION ', j, &
1046 'MUST BE >= cell bottom (', bot,
')'
1048 else if (this%telev(ipos) > top)
then
1049 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
1050 'telev FOR LAKE ', n,
' HORIZONTAL CONNECTION ', j, &
1051 'MUST BE <= cell top (', top,
')'
1055 this%laketop(n) = max(this%telev(ipos), this%laketop(n))
1056 this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
1058 else if (this%ictype(ipos) == 2 .or. this%ictype(ipos) == 3)
then
1059 this%telev(ipos) = top
1060 this%belev(ipos) = bot
1061 this%lakebot(n) = bot
1065 if (nboundchk(ipos) == 0)
then
1066 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1067 'NO DATA SPECIFIED FOR LAKE', n,
'CONNECTION', j
1069 else if (nboundchk(ipos) > 1)
then
1070 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
1071 'DATA FOR LAKE', n,
'CONNECTION', j, &
1072 'SPECIFIED', nboundchk(ipos),
'TIMES'
1078 if (this%laketop(n) == -dep20)
then
1079 this%laketop(n) = this%lakebot(n) + 100.
1084 deallocate (nboundchk)
1088 call this%parser%StoreErrorUnit()
1090 end subroutine lak_read_lake_connections
1094 subroutine lak_read_tables(this)
1098 class(laktype),
intent(inout) :: this
1100 type(laktabtype),
dimension(:),
allocatable :: laketables
1101 character(len=LINELENGTH) :: line
1102 character(len=LINELENGTH) :: keyword
1103 integer(I4B) :: ierr
1104 logical(LGP) :: isfound, endOfBlock
1106 integer(I4B) :: iconn
1107 integer(I4B) :: ntabs
1108 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
1111 if (this%ntables < 1)
return
1114 allocate (nboundchk(this%nlakes))
1115 do n = 1, this%nlakes
1120 allocate (laketables(this%nlakes))
1123 call this%parser%GetBlock(
'TABLES', isfound, ierr, &
1124 supportopenclose=.true.)
1130 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
1133 call this%parser%GetNextLine(endofblock)
1134 if (endofblock)
exit
1135 n = this%parser%GetInteger()
1137 if (n < 1 .or. n > this%nlakes)
then
1138 write (
errmsg,
'(a,1x,i0)')
'lakeno MUST BE > 0 and <= ', this%nlakes
1145 nboundchk(n) = nboundchk(n) + 1
1148 call this%parser%GetStringCaps(keyword)
1149 select case (keyword)
1151 call this%parser%GetStringCaps(keyword)
1152 if (trim(adjustl(keyword)) /=
'FILEIN')
then
1153 errmsg =
'TAB6 keyword must be followed by "FILEIN" '// &
1158 call this%parser%GetString(line)
1159 call this%lak_read_table(n, line, laketables(n))
1161 write (
errmsg,
'(a,1x,i0,1x,a)') &
1162 'LAKE TABLE ENTRY for LAKE ', n,
'MUST INCLUDE TAB6 KEYWORD'
1168 write (this%iout,
'(1x,a)') &
1169 'END OF '//trim(adjustl(this%text))//
' LAKE_TABLES'
1172 if (ntabs < this%ntables)
then
1173 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1174 'TABLE DATA ARE SPECIFIED', ntabs, &
1175 'TIMES BUT NTABLES IS SET TO', this%ntables
1178 do n = 1, this%nlakes
1179 if (this%ntabrow(n) > 0 .and. nboundchk(n) > 1)
then
1180 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1181 'TABLE DATA FOR LAKE', n,
'SPECIFIED', nboundchk(n),
'TIMES'
1186 call store_error(
'REQUIRED TABLES BLOCK NOT FOUND.')
1190 deallocate (nboundchk)
1194 call this%parser%StoreErrorUnit()
1198 call this%laktables_to_vectors(laketables)
1201 do n = 1, this%nlakes
1202 if (this%ntabrow(n) > 0)
then
1203 deallocate (laketables(n)%tabstage)
1204 deallocate (laketables(n)%tabvolume)
1205 deallocate (laketables(n)%tabsarea)
1206 iconn = this%idxlakeconn(n)
1207 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1208 deallocate (laketables(n)%tabwarea)
1212 deallocate (laketables)
1213 end subroutine lak_read_tables
1218 subroutine laktables_to_vectors(this, laketables)
1219 class(laktype),
intent(inout) :: this
1220 type(laktabtype),
intent(in),
dimension(:),
contiguous :: laketables
1222 integer(I4B) :: ntabrows
1224 integer(I4B) :: ipos
1225 integer(I4B) :: iconn
1228 call mem_allocate(this%ialaktab, this%nlakes + 1,
'IALAKTAB', this%memoryPath)
1231 this%ialaktab(1) = 1
1232 do n = 1, this%nlakes
1234 this%ialaktab(n + 1) = this%ialaktab(n) + this%ntabrow(n)
1238 ntabrows = this%ialaktab(this%nlakes + 1) - 1
1239 call mem_allocate(this%tabstage, ntabrows,
'TABSTAGE', this%memoryPath)
1240 call mem_allocate(this%tabvolume, ntabrows,
'TABVOLUME', this%memoryPath)
1241 call mem_allocate(this%tabsarea, ntabrows,
'TABSAREA', this%memoryPath)
1242 call mem_allocate(this%tabwarea, ntabrows,
'TABWAREA', this%memoryPath)
1245 do n = 1, this%nlakes
1247 do ipos = this%ialaktab(n), this%ialaktab(n + 1) - 1
1248 this%tabstage(ipos) = laketables(n)%tabstage(j)
1249 this%tabvolume(ipos) = laketables(n)%tabvolume(j)
1250 this%tabsarea(ipos) = laketables(n)%tabsarea(j)
1251 iconn = this%idxlakeconn(n)
1252 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1255 this%tabwarea(ipos) = laketables(n)%tabwarea(j)
1257 this%tabwarea(ipos) =
dzero
1262 end subroutine laktables_to_vectors
1266 subroutine lak_read_table(this, ilak, filename, laketable)
1271 class(laktype),
intent(inout) :: this
1272 integer(I4B),
intent(in) :: ilak
1273 character(len=*),
intent(in) :: filename
1274 type(laktabtype),
intent(inout) :: laketable
1276 character(len=LINELENGTH) :: keyword
1277 integer(I4B) :: ierr
1278 logical(LGP) :: isfound, endOfBlock
1281 integer(I4B) :: ipos
1283 integer(I4B) :: jmin
1284 integer(I4B) :: iconn
1292 character(len=*),
parameter :: fmttaberr = &
1293 &
'(a,1x,i0,1x,a,1x,g15.6,1x,a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.6,1x,a)'
1301 call openfile(iu, this%iout, filename,
'LAKE TABLE')
1302 call parser%Initialize(iu, this%iout)
1305 call parser%GetBlock(
'DIMENSIONS', isfound, ierr, supportopenclose=.true.)
1310 if (this%iprpak /= 0)
then
1311 write (this%iout,
'(/1x,a)') &
1312 'PROCESSING '//trim(adjustl(this%text))//
' DIMENSIONS'
1315 call parser%GetNextLine(endofblock)
1316 if (endofblock)
exit
1317 call parser%GetStringCaps(keyword)
1318 select case (keyword)
1320 n = parser%GetInteger()
1323 write (
errmsg,
'(a)')
'LAKE TABLE NROW MUST BE > 0'
1327 j = parser%GetInteger()
1329 if (this%ictype(ilak) == 2 .or. this%ictype(ilak) == 3)
then
1335 write (
errmsg,
'(a,1x,i0)')
'LAKE TABLE NCOL MUST BE >= ', jmin
1340 write (
errmsg,
'(a,a)') &
1341 'UNKNOWN '//trim(this%text)//
' DIMENSIONS KEYWORD: ', trim(keyword)
1345 if (this%iprpak /= 0)
then
1346 write (this%iout,
'(1x,a)') &
1347 'END OF '//trim(adjustl(this%text))//
' DIMENSIONS'
1350 call store_error(
'REQUIRED DIMENSIONS BLOCK NOT FOUND.')
1356 'NROW NOT SPECIFIED IN THE LAKE TABLE DIMENSIONS BLOCK'
1361 'NCOL NOT SPECIFIED IN THE LAKE TABLE DIMENSIONS BLOCK'
1370 this%ntabrow(ilak) = n
1371 allocate (laketable%tabstage(n))
1372 allocate (laketable%tabvolume(n))
1373 allocate (laketable%tabsarea(n))
1374 ipos = this%idxlakeconn(ilak)
1375 if (this%ictype(ipos) == 2 .or. this%ictype(ipos) == 3)
then
1376 allocate (laketable%tabwarea(n))
1380 call parser%GetBlock(
'TABLE', isfound, ierr, supportopenclose=.true.)
1386 if (this%iprpak /= 0)
then
1387 write (this%iout,
'(/1x,a)') &
1388 'PROCESSING '//trim(adjustl(this%text))//
' TABLE'
1390 iconn = this%idxlakeconn(ilak)
1393 call parser%GetNextLine(endofblock)
1394 if (endofblock)
exit
1396 if (ipos > this%ntabrow(ilak))
then
1399 laketable%tabstage(ipos) = parser%GetDouble()
1400 laketable%tabvolume(ipos) = parser%GetDouble()
1401 laketable%tabsarea(ipos) = parser%GetDouble()
1402 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1403 laketable%tabwarea(ipos) = parser%GetDouble()
1405 end do readtabledata
1407 if (this%iprpak /= 0)
then
1408 write (this%iout,
'(1x,a)') &
1409 'END OF '//trim(adjustl(this%text))//
' TABLE'
1412 call store_error(
'REQUIRED TABLE BLOCK NOT FOUND.')
1416 if (ipos /= this%ntabrow(ilak))
then
1417 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1418 'NROW SET TO', this%ntabrow(ilak),
'BUT', ipos,
'ROWS WERE READ'
1423 iconn = this%idxlakeconn(ilak)
1424 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1425 do n = 1, this%ntabrow(ilak)
1426 vol = laketable%tabvolume(n)
1427 sa = laketable%tabsarea(n)
1428 wa = laketable%tabwarea(n)
1431 if (vol > dzero)
exit
1433 this%lakebot(ilak) = laketable%tabstage(n)
1434 this%belev(ilak) = laketable%tabstage(n)
1437 n = this%ntabrow(ilak)
1438 this%sareamax(ilak) = laketable%tabsarea(n)
1442 do n = 2, this%ntabrow(ilak)
1443 v = laketable%tabstage(n)
1444 v0 = laketable%tabstage(n - 1)
1446 write (
errmsg, fmttaberr) &
1447 'TABLE STAGE ENTRY', n,
'(', laketable%tabstage(n),
') FOR LAKE ', &
1448 ilak,
'MUST BE GREATER THAN THE PREVIOUS STAGE ENTRY', &
1449 n - 1,
'(', laketable%tabstage(n - 1),
')'
1452 v = laketable%tabvolume(n)
1453 v0 = laketable%tabvolume(n - 1)
1455 write (
errmsg, fmttaberr) &
1456 'TABLE VOLUME ENTRY', n,
'(', laketable%tabvolume(n), &
1458 ilak,
'MUST BE GREATER THAN THE PREVIOUS VOLUME ENTRY', &
1459 n - 1,
'(', laketable%tabvolume(n - 1),
')'
1462 v = laketable%tabsarea(n)
1463 v0 = laketable%tabsarea(n - 1)
1465 write (
errmsg, fmttaberr) &
1466 'TABLE SURFACE AREA ENTRY', n,
'(', &
1467 laketable%tabsarea(n),
') FOR LAKE ', ilak, &
1468 'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS SURFACE AREA ENTRY', &
1469 n - 1,
'(', laketable%tabsarea(n - 1),
')'
1472 iconn = this%idxlakeconn(ilak)
1473 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1474 v = laketable%tabwarea(n)
1475 v0 = laketable%tabwarea(n - 1)
1477 write (
errmsg, fmttaberr) &
1478 'TABLE EXCHANGE AREA ENTRY', n,
'(', &
1479 laketable%tabwarea(n),
') FOR LAKE ', ilak, &
1480 'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS EXCHANGE AREA '// &
1481 'ENTRY', n - 1,
'(', laketable%tabwarea(n - 1),
')'
1490 call parser%StoreErrorUnit()
1495 end subroutine lak_read_table
1499 subroutine lak_read_outlets(this)
1504 class(laktype),
intent(inout) :: this
1506 character(len=LINELENGTH) :: text, keyword
1507 character(len=LENBOUNDNAME) :: bndName
1508 character(len=9) :: citem
1509 integer(I4B) :: ierr, ival
1510 logical(LGP) :: isfound, endOfBlock
1513 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
1514 real(DP),
pointer :: bndElem => null()
1517 call this%parser%GetBlock(
'OUTLETS', isfound, ierr, &
1518 supportopenclose=.true., blockrequired=.false.)
1522 if (this%noutlets > 0)
then
1525 allocate (nboundchk(this%noutlets))
1526 do n = 1, this%noutlets
1531 call mem_allocate(this%lakein, this%NOUTLETS,
'LAKEIN', this%memoryPath)
1532 call mem_allocate(this%lakeout, this%NOUTLETS,
'LAKEOUT', this%memoryPath)
1533 call mem_allocate(this%iouttype, this%NOUTLETS,
'IOUTTYPE', &
1535 call mem_allocate(this%outrate, this%NOUTLETS,
'OUTRATE', this%memoryPath)
1536 call mem_allocate(this%outinvert, this%NOUTLETS,
'OUTINVERT', &
1538 call mem_allocate(this%outwidth, this%NOUTLETS,
'OUTWIDTH', &
1540 call mem_allocate(this%outrough, this%NOUTLETS,
'OUTROUGH', &
1542 call mem_allocate(this%outslope, this%NOUTLETS,
'OUTSLOPE', &
1544 call mem_allocate(this%simoutrate, this%NOUTLETS,
'SIMOUTRATE', &
1548 do n = 1, this%noutlets
1549 this%outrate(n) = dzero
1553 write (this%iout,
'(/1x,a)') &
1554 'PROCESSING '//trim(adjustl(this%text))//
' OUTLETS'
1556 call this%parser%GetNextLine(endofblock)
1557 if (endofblock)
exit
1558 n = this%parser%GetInteger()
1560 if (n < 1 .or. n > this%noutlets)
then
1561 write (
errmsg,
'(a,1x,i0)') &
1562 'outletno MUST BE > 0 and <= ', this%noutlets
1568 nboundchk(n) = nboundchk(n) + 1
1571 ival = this%parser%GetInteger()
1572 if (ival < 1 .or. ival > this%nlakes)
then
1573 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1574 'lakein FOR OUTLET ', n,
'MUST BE > 0 and <= ', this%nlakes
1578 this%lakein(n) = ival
1581 ival = this%parser%GetInteger()
1582 if (ival < 0 .or. ival > this%nlakes)
then
1583 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1584 'lakeout FOR OUTLET ', n,
'MUST BE >= 0 and <= ', this%nlakes
1588 this%lakeout(n) = ival
1591 call this%parser%GetStringCaps(keyword)
1592 select case (keyword)
1594 this%iouttype(n) = 0
1596 this%iouttype(n) = 1
1598 this%iouttype(n) = 2
1600 write (
errmsg,
'(a,1x,i0,1x,a,a,a)') &
1601 'UNKNOWN couttype FOR OUTLET ', n,
'(', trim(keyword),
')'
1607 write (citem,
'(i9.9)') n
1608 bndname =
'OUTLET'//citem
1614 call this%parser%GetString(text)
1615 bndelem => this%outinvert(n)
1617 this%packName,
'BND', &
1618 this%tsManager, this%iprpak, &
1622 call this%parser%GetString(text)
1623 bndelem => this%outwidth(n)
1625 this%packName,
'BND', &
1626 this%tsManager, this%iprpak,
'WIDTH')
1629 call this%parser%GetString(text)
1630 bndelem => this%outrough(n)
1632 this%packName,
'BND', &
1633 this%tsManager, this%iprpak,
'ROUGH')
1636 call this%parser%GetString(text)
1637 bndelem => this%outslope(n)
1639 this%packName,
'BND', &
1640 this%tsManager, this%iprpak,
'SLOPE')
1642 write (this%iout,
'(1x,a)')
'END OF '//trim(adjustl(this%text))// &
1646 do n = 1, this%noutlets
1647 if (nboundchk(n) == 0)
then
1648 write (
errmsg,
'(a,1x,i0)')
'NO DATA SPECIFIED FOR OUTLET', n
1650 else if (nboundchk(n) > 1)
then
1651 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1652 'DATA FOR OUTLET', n,
'SPECIFIED', nboundchk(n),
'TIMES'
1658 deallocate (nboundchk)
1660 write (
errmsg,
'(a,1x,a)') &
1661 'AN OUTLETS BLOCK SHOULD NOT BE SPECIFIED IF NOUTLETS IS NOT', &
1662 'SPECIFIED OR IS SPECIFIED TO BE 0.'
1667 if (this%noutlets > 0)
then
1668 call store_error(
'REQUIRED OUTLETS BLOCK NOT FOUND.')
1675 call this%parser%StoreErrorUnit()
1677 end subroutine lak_read_outlets
1681 subroutine lak_read_dimensions(this)
1685 class(laktype),
intent(inout) :: this
1687 character(len=LINELENGTH) :: keyword
1688 integer(I4B) :: ierr
1689 logical(LGP) :: isfound, endOfBlock
1696 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
1697 supportopenclose=.true.)
1701 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
1704 call this%parser%GetNextLine(endofblock)
1705 if (endofblock)
exit
1706 call this%parser%GetStringCaps(keyword)
1707 select case (keyword)
1709 this%nlakes = this%parser%GetInteger()
1710 write (this%iout,
'(4x,a,i7)')
'NLAKES = ', this%nlakes
1712 this%noutlets = this%parser%GetInteger()
1713 write (this%iout,
'(4x,a,i7)')
'NOUTLETS = ', this%noutlets
1715 this%ntables = this%parser%GetInteger()
1716 write (this%iout,
'(4x,a,i7)')
'NTABLES = ', this%ntables
1718 write (
errmsg,
'(a,a)') &
1719 'UNKNOWN '//trim(this%text)//
' DIMENSION: ', trim(keyword)
1723 write (this%iout,
'(1x,a)') &
1724 'END OF '//trim(adjustl(this%text))//
' DIMENSIONS'
1726 call store_error(
'REQUIRED DIMENSIONS BLOCK NOT FOUND.')
1729 if (this%nlakes < 0)
then
1731 'NLAKES WAS NOT SPECIFIED OR WAS SPECIFIED INCORRECTLY.'
1735 if (this%iforceleglak /= 0)
then
1736 if (this%iforceleglak < 1 .or. this%iforceleglak > this%nlakes)
then
1737 write (
errmsg,
'(a,i0,a,i0,a)') &
1738 'DEV_FORCE_LEGACY_LAKE (', this%iforceleglak, &
1739 ') MUST BE BETWEEN 1 AND NLAKES (', this%nlakes,
').'
1746 call this%parser%StoreErrorUnit()
1751 if (this%iimplicit /= 0)
then
1752 this%npakeq = this%nlakes
1761 call this%lak_read_lakes()
1764 call this%lak_read_lake_connections()
1767 call this%lak_read_tables()
1770 call this%lak_read_outlets()
1774 call this%define_listlabel()
1777 call this%lak_setup_budobj()
1780 call this%lak_setup_tableobj()
1781 end subroutine lak_read_dimensions
1785 subroutine lak_read_initial_attr(this)
1791 class(laktype),
intent(inout) :: this
1793 character(len=LINELENGTH) :: text
1794 integer(I4B) :: j, jj, n
1811 real(DP),
allocatable,
dimension(:) :: clb, caq
1812 character(len=14) :: cbedleak
1813 character(len=14) :: cbedcond
1814 character(len=10),
dimension(0:3) :: ctype
1815 character(len=15) :: nodestr
1816 real(DP),
pointer :: bndElem => null()
1818 data ctype(0)/
'VERTICAL '/
1819 data ctype(1)/
'HORIZONTAL'/
1820 data ctype(2)/
'EMBEDDEDH '/
1821 data ctype(3)/
'EMBEDDEDV '/
1824 do n = 1, this%nlakes
1825 this%xnewpak(n) = this%strt(n)
1826 write (text,
'(g15.7)') this%strt(n)
1828 bndelem => this%stage(n)
1830 'BND', this%tsManager, this%iprpak, &
1835 do n = 1, this%nlakes
1836 if (this%status(n) ==
'CONSTANT')
then
1837 this%iboundpak(n) = -1
1838 else if (this%status(n) ==
'INACTIVE')
then
1839 this%iboundpak(n) = 0
1840 else if (this%status(n) ==
'ACTIVE ')
then
1841 this%iboundpak(n) = 1
1846 if (this%inamedbound /= 0)
then
1847 do n = 1, this%nlakes
1848 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1849 this%boundname(j) = this%lakename(n)
1855 call this%copy_boundname()
1865 allocate (clb(this%MAXBOUND))
1866 allocate (caq(this%MAXBOUND))
1869 do n = 1, this%nlakes
1870 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1872 top = this%dis%top(nn)
1873 bot = this%dis%bot(nn)
1875 if (this%ictype(j) == 0)
then
1876 area = this%dis%area(nn)
1877 this%sarea(j) = area
1878 this%warea(j) = area
1879 this%sareamax(n) = this%sareamax(n) + area
1880 if (this%gwfik33 == 0)
then
1885 length = dhalf * (top - bot)
1887 else if (this%ictype(j) == 1)
then
1888 area = (this%telev(j) - this%belev(j)) * this%connwidth(j)
1891 if (top == this%telev(j) .and. bot == this%belev(j))
then
1892 if (this%icelltype(nn) == 0)
then
1893 area = this%gwfsat(nn) * (top - bot) * this%connwidth(j)
1896 this%sarea(j) = dzero
1897 this%warea(j) = area
1898 this%sareamax(n) = this%sareamax(n) + dzero
1900 length = this%connlength(j)
1902 else if (this%ictype(j) == 2)
then
1904 this%sarea(j) = dzero
1905 this%warea(j) = area
1906 this%sareamax(n) = this%sareamax(n) + dzero
1908 length = this%connlength(j)
1910 else if (this%ictype(j) == 3)
then
1912 this%sarea(j) = dzero
1913 this%warea(j) = area
1914 this%sareamax(n) = this%sareamax(n) + dzero
1915 if (this%gwfik33 == 0)
then
1920 length = this%connlength(j)
1922 if (
is_close(this%bedleak(j), dnodata))
then
1924 else if (this%bedleak(j) > dzero)
then
1925 clb(j) = done / this%bedleak(j)
1934 if (
is_close(this%bedleak(j), dnodata))
then
1935 this%satcond(j) = area / caq(j)
1936 else if (clb(j) * caq(j) > dzero)
then
1937 this%satcond(j) = area / (clb(j) + caq(j))
1939 this%satcond(j) = dzero
1945 if (this%iprpak > 0)
then
1946 write (this%iout,
'(//,29x,a,/)') &
1947 'INTERFACE CONDUCTANCE BETWEEN LAKE AND AQUIFER CELLS'
1948 write (this%iout,
'(1x,a)') &
1949 &
' LAKE CONNECTION CONNECTION LAKEBED'// &
1950 &
' C O N D U C T A N C E S '
1951 write (this%iout,
'(1x,a)') &
1952 &
' NUMBER NUMBER CELLID DIRECTION LEAKANCE'// &
1953 &
' LAKEBED AQUIFER COMBINED'
1954 write (this%iout,
"(1x,108('-'))")
1955 do n = 1, this%nlakes
1957 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1960 if (this%ictype(j) == 1)
then
1961 fact = this%telev(j) - this%belev(j)
1962 if (abs(fact) > dzero)
then
1967 area = this%warea(j)
1969 if (
is_close(clb(j), dnodata))
then
1972 else if (clb(j) > dzero)
then
1973 c1 = area * fact / clb(j)
1974 write (cbedleak,
'(g14.5)') this%bedleak(j)
1975 write (cbedcond,
'(g14.5)') c1
1977 write (cbedleak,
'(g14.5)') c1
1978 write (cbedcond,
'(g14.5)') c1
1981 if (caq(j) > dzero)
then
1982 c2 = area * fact / caq(j)
1984 call this%dis%noder_to_string(nn, nodestr)
1986 '(1x,i10,1x,i10,1x,a15,1x,a10,2(1x,a14),2(1x,g14.5))') &
1987 n, idx, nodestr, ctype(this%ictype(j)), cbedleak, &
1988 cbedcond, c2, this%satcond(j) * fact
1991 write (this%iout,
"(1x,108('-'))")
1992 write (this%iout,
'(1x,a)') &
1993 'IF VERTICAL CONNECTION, CONDUCTANCE (L^2/T) IS &
1994 &BETWEEN AQUIFER CELL AND OVERLYING LAKE CELL.'
1995 write (this%iout,
'(1x,a)') &
1996 'IF HORIZONTAL CONNECTION, CONDUCTANCES ARE PER &
1997 &UNIT SATURATED THICKNESS (L/T).'
1998 write (this%iout,
'(1x,a)') &
1999 'IF EMBEDDED CONNECTION, CONDUCTANCES ARE PER &
2000 &UNIT EXCHANGE AREA (1/T).'
2005 do n = 1, this%nlakes
2006 write (this%iout,
'(//1x,a,1x,i10)')
'STAGE/VOLUME RELATION FOR LAKE ', n
2007 write (this%iout,
'(/1x,5(a14))')
' STAGE',
' SURFACE AREA', &
2008 &
' WETTED AREA',
' CONDUCTANCE', &
2010 write (this%iout,
"(1x,70('-'))")
2011 dx = (this%laketop(n) - this%lakebot(n)) / 150.
2014 call this%lak_calculate_conductance(n, s, c)
2015 call this%lak_calculate_sarea(n, s, sa)
2016 call this%lak_calculate_warea(n, s, wa, s)
2017 call this%lak_calculate_vol(n, s, v)
2018 write (this%iout,
'(1x,5(E14.5))') s, sa, wa, c, v
2021 write (this%iout,
"(1x,70('-'))")
2023 write (this%iout,
'(//1x,a,1x,i10)')
'STAGE/VOLUME RELATION FOR LAKE ', n
2024 write (this%iout,
'(/1x,4(a14))')
' ',
' ', &
2025 &
' CALCULATED',
' STAGE'
2026 write (this%iout,
'(1x,4(a14))')
' STAGE',
' VOLUME', &
2027 &
' STAGE',
' DIFFERENCE'
2028 write (this%iout,
"(1x,56('-'))")
2029 s = this%lakebot(n) - dx
2031 call this%lak_calculate_vol(n, s, v)
2032 call this%lak_vol2stage(n, v, c)
2033 write (this%iout,
'(1x,4(E14.5))') s, v, c, s - c
2036 write (this%iout,
"(1x,56('-'))")
2041 this%gwfk11 => null()
2042 this%gwfk33 => null()
2043 this%gwfsat => null()
2044 this%gwfik33 => null()
2049 end subroutine lak_read_initial_attr
2055 subroutine lak_linear_interpolation(this, n, x, y, z, v)
2057 class(laktype),
intent(inout) :: this
2058 integer(I4B),
intent(in) :: n
2059 real(DP),
dimension(n),
intent(in) :: x
2060 real(DP),
dimension(n),
intent(in) :: y
2061 real(DP),
intent(in) :: z
2062 real(DP),
intent(inout) :: v
2065 real(DP) :: dx, dydx
2073 else if (z > x(n))
then
2074 dx = x(n) - x(n - 1)
2076 if (abs(dx) >
dzero)
then
2077 dydx = (y(n) - y(n - 1)) / dx
2080 v = y(n) + dydx * dx
2084 dx = x(i) - x(i - 1)
2086 if (z >= x(i - 1) .and. z <= x(i))
then
2087 if (abs(dx) >
dzero)
then
2088 dydx = (y(i) - y(i - 1)) / dx
2091 v = y(i - 1) + dydx * dx
2096 end subroutine lak_linear_interpolation
2100 subroutine lak_calculate_sarea(this, ilak, stage, sarea)
2102 class(laktype),
intent(inout) :: this
2103 integer(I4B),
intent(in) :: ilak
2104 real(DP),
intent(in) :: stage
2105 real(DP),
intent(inout) :: sarea
2108 integer(I4B) :: ifirst
2109 integer(I4B) :: ilast
2116 i = this%ntabrow(ilak)
2118 ifirst = this%ialaktab(ilak)
2119 ilast = this%ialaktab(ilak + 1) - 1
2120 if (stage <= this%tabstage(ifirst))
then
2121 sarea = this%tabsarea(ifirst)
2122 else if (stage >= this%tabstage(ilast))
then
2123 sarea = this%tabsarea(ilast)
2125 call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2126 this%tabsarea(ifirst:ilast), &
2130 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2131 topl = this%telev(i)
2132 botl = this%belev(i)
2134 sa = sat * this%sarea(i)
2138 end subroutine lak_calculate_sarea
2142 subroutine lak_calculate_warea(this, ilak, stage, warea, hin)
2144 class(laktype),
intent(inout) :: this
2145 integer(I4B),
intent(in) :: ilak
2146 real(DP),
intent(in) :: stage
2147 real(DP),
intent(inout) :: warea
2148 real(DP),
optional,
intent(inout) :: hin
2151 integer(I4B) :: igwfnode
2156 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2157 if (
present(hin))
then
2160 igwfnode = this%cellid(i)
2161 head = this%xnew(igwfnode)
2163 call this%lak_calculate_conn_warea(ilak, i, stage, head, wa)
2166 end subroutine lak_calculate_warea
2170 subroutine lak_calculate_conn_warea(this, ilak, iconn, stage, head, wa)
2172 class(laktype),
intent(inout) :: this
2173 integer(I4B),
intent(in) :: ilak
2174 integer(I4B),
intent(in) :: iconn
2175 real(DP),
intent(in) :: stage
2176 real(DP),
intent(in) :: head
2177 real(DP),
intent(inout) :: wa
2180 integer(I4B) :: ifirst
2181 integer(I4B) :: ilast
2182 integer(I4B) :: node
2189 topl = this%telev(iconn)
2190 botl = this%belev(iconn)
2191 call this%lak_calculate_cond_head(iconn, stage, head, vv)
2192 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
2193 if (vv > topl) vv = topl
2194 i = this%ntabrow(ilak)
2195 ifirst = this%ialaktab(ilak)
2196 ilast = this%ialaktab(ilak + 1) - 1
2197 if (vv <= this%tabstage(ifirst))
then
2198 wa = this%tabwarea(ifirst)
2199 else if (vv >= this%tabstage(ilast))
then
2200 wa = this%tabwarea(ilast)
2202 call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2203 this%tabwarea(ifirst:ilast), &
2207 node = this%cellid(iconn)
2209 if (this%icelltype(node) == 0)
then
2215 wa = sat * this%warea(iconn)
2217 end subroutine lak_calculate_conn_warea
2221 subroutine lak_calculate_vol(this, ilak, stage, volume)
2223 class(laktype),
intent(inout) :: this
2224 integer(I4B),
intent(in) :: ilak
2225 real(DP),
intent(in) :: stage
2226 real(DP),
intent(inout) :: volume
2229 integer(I4B) :: ifirst
2230 integer(I4B) :: ilast
2239 i = this%ntabrow(ilak)
2241 ifirst = this%ialaktab(ilak)
2242 ilast = this%ialaktab(ilak + 1) - 1
2243 if (stage <= this%tabstage(ifirst))
then
2244 volume = this%tabvolume(ifirst)
2245 else if (stage >= this%tabstage(ilast))
then
2246 ds = stage - this%tabstage(ilast)
2247 sa = this%tabsarea(ilast)
2248 volume = this%tabvolume(ilast) + ds * sa
2250 call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2251 this%tabvolume(ifirst:ilast), &
2255 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2256 topl = this%telev(i)
2257 botl = this%belev(i)
2259 sa = sat * this%sarea(i)
2260 if (stage < botl)
then
2262 else if (stage > botl .and. stage < topl)
then
2263 v = sa * (stage - botl)
2265 v = sa * (topl - botl) + sa * (stage - topl)
2270 end subroutine lak_calculate_vol
2274 subroutine lak_calculate_conductance(this, ilak, stage, conductance)
2276 class(laktype),
intent(inout) :: this
2277 integer(I4B),
intent(in) :: ilak
2278 real(DP),
intent(in) :: stage
2279 real(DP),
intent(inout) :: conductance
2285 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2286 call this%lak_calculate_conn_conductance(ilak, i, stage, stage, c)
2287 conductance = conductance + c
2289 end subroutine lak_calculate_conductance
2295 subroutine lak_calculate_cond_head(this, iconn, stage, head, vv)
2297 class(laktype),
intent(inout) :: this
2298 integer(I4B),
intent(in) :: iconn
2299 real(DP),
intent(in) :: stage
2300 real(DP),
intent(in) :: head
2301 real(DP),
intent(inout) :: vv
2308 topl = this%telev(iconn)
2309 botl = this%belev(iconn)
2310 ss = min(stage, topl)
2311 hh = min(head, topl)
2312 if (this%igwhcopt > 0)
then
2314 else if (this%inewton > 0)
then
2317 vv =
dhalf * (ss + hh)
2319 end subroutine lak_calculate_cond_head
2324 subroutine lak_calculate_conn_conductance(this, ilak, iconn, stage, head, cond)
2326 class(laktype),
intent(inout) :: this
2327 integer(I4B),
intent(in) :: ilak
2328 integer(I4B),
intent(in) :: iconn
2329 real(DP),
intent(in) :: stage
2330 real(DP),
intent(in) :: head
2331 real(DP),
intent(inout) :: cond
2333 integer(I4B) :: node
2341 real(DP) :: vscratio
2345 topl = this%telev(iconn)
2346 botl = this%belev(iconn)
2347 call this%lak_calculate_cond_head(iconn, stage, head, vv)
2352 if (this%ictype(iconn) == 0)
then
2353 if (abs(topl - botl) <
dprec)
then
2358 else if (this%ictype(iconn) == 1)
then
2359 node = this%cellid(iconn)
2360 if (this%icelltype(node) == 0)
then
2364 else if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
2365 node = this%cellid(iconn)
2366 if (this%icelltype(node) == 0)
then
2367 vv = this%telev(iconn)
2368 call this%lak_calculate_conn_warea(ilak, iconn, vv, vv, wa)
2370 call this%lak_calculate_conn_warea(ilak, iconn, stage, head, wa)
2376 if (this%ivsc == 1)
then
2378 if (stage > head)
then
2379 vscratio = this%viscratios(1, iconn)
2382 vscratio = this%viscratios(2, iconn)
2385 cond = sat * this%satcond(iconn) * vscratio
2386 end subroutine lak_calculate_conn_conductance
2390 subroutine lak_calculate_exchange(this, ilak, stage, totflow)
2392 class(laktype),
intent(inout) :: this
2393 integer(I4B),
intent(in) :: ilak
2394 real(DP),
intent(in) :: stage
2395 real(DP),
intent(inout) :: totflow
2398 integer(I4B) :: igwfnode
2403 do j = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2404 igwfnode = this%cellid(j)
2405 hgwf = this%xnew(igwfnode)
2406 call this%lak_calculate_conn_exchange(ilak, j, stage, hgwf, flow)
2407 totflow = totflow + flow
2409 end subroutine lak_calculate_exchange
2414 subroutine lak_calculate_conn_exchange(this, ilak, iconn, stage, head, flow, &
2417 class(laktype),
intent(inout) :: this
2418 integer(I4B),
intent(in) :: ilak
2419 integer(I4B),
intent(in) :: iconn
2420 real(DP),
intent(in) :: stage
2421 real(DP),
intent(in) :: head
2422 real(DP),
intent(inout) :: flow
2423 real(DP),
intent(inout),
optional :: gwfhcof
2424 real(DP),
intent(inout),
optional :: gwfrhs
2430 real(DP) :: gwfhcof0
2434 call this%lak_calculate_conn_conductance(ilak, iconn, stage, head, cond)
2435 botl = this%belev(iconn)
2438 if (stage >= botl)
then
2445 if (head >= botl)
then
2452 flow = cond * (hh - ss)
2455 if (head >= botl)
then
2457 gwfrhs0 = -cond * ss
2464 if (this%idense /= 0)
then
2465 call this%lak_calculate_density_exchange(iconn, stage, head, cond, botl, &
2466 flow, gwfhcof0, gwfrhs0)
2470 if (
present(gwfhcof)) gwfhcof = gwfhcof0
2471 if (
present(gwfrhs)) gwfrhs = gwfrhs0
2472 end subroutine lak_calculate_conn_exchange
2483 subroutine lak_calculate_conn_exchange_deriv(this, ilak, iconn, stage, &
2484 head, flow, dqds, dqdh)
2486 class(laktype),
intent(inout) :: this
2487 integer(I4B),
intent(in) :: ilak
2488 integer(I4B),
intent(in) :: iconn
2489 real(DP),
intent(in) :: stage
2490 real(DP),
intent(in) :: head
2491 real(DP),
intent(inout) :: flow
2492 real(DP),
intent(inout),
optional :: dqds
2493 real(DP),
intent(inout),
optional :: dqdh
2495 real(DP) :: cond, botl, dps, dph, ss, hh
2497 call this%lak_calculate_conn_conductance(ilak, iconn, stage, head, cond)
2498 botl = this%belev(iconn)
2499 if (stage >= botl)
then
2506 if (head >= botl)
then
2513 flow = cond * (hh - ss)
2514 if (
present(dqds)) dqds = -cond * dps
2515 if (
present(dqdh)) dqdh = cond * dph
2516 end subroutine lak_calculate_conn_exchange_deriv
2521 subroutine lak_estimate_conn_exchange(this, iflag, ilak, iconn, idry, stage, &
2522 head, flow, source, gwfhcof, gwfrhs)
2524 class(laktype),
intent(inout) :: this
2525 integer(I4B),
intent(in) :: iflag
2526 integer(I4B),
intent(in) :: ilak
2527 integer(I4B),
intent(in) :: iconn
2528 integer(I4B),
intent(inout) :: idry
2529 real(DP),
intent(in) :: stage
2530 real(DP),
intent(in) :: head
2531 real(DP),
intent(inout) :: flow
2532 real(DP),
intent(inout) :: source
2533 real(DP),
intent(inout),
optional :: gwfhcof
2534 real(DP),
intent(inout),
optional :: gwfrhs
2536 real(DP) :: gwfhcof0, gwfrhs0
2540 call this%lak_calculate_conn_exchange(ilak, iconn, stage, head, flow, &
2542 if (iflag == 1)
then
2543 if (flow >
dzero)
then
2544 source = source + flow
2546 else if (iflag == 2)
then
2547 if (-flow > source)
then
2551 else if (flow <
dzero)
then
2552 source = source + flow
2557 if (
present(gwfhcof)) gwfhcof = gwfhcof0
2558 if (
present(gwfrhs)) gwfrhs = gwfrhs0
2559 end subroutine lak_estimate_conn_exchange
2564 subroutine lak_calculate_storagechange(this, ilak, stage, stage0, delt, dvr)
2566 class(laktype),
intent(inout) :: this
2567 integer(I4B),
intent(in) :: ilak
2568 real(DP),
intent(in) :: stage
2569 real(DP),
intent(in) :: stage0
2570 real(DP),
intent(in) :: delt
2571 real(DP),
intent(inout) :: dvr
2577 if (this%gwfiss /= 1)
then
2578 call this%lak_calculate_vol(ilak, stage, v)
2579 call this%lak_calculate_vol(ilak, stage0, v0)
2580 dvr = (v0 - v) / delt
2582 end subroutine lak_calculate_storagechange
2586 subroutine lak_calculate_rainfall(this, ilak, stage, ra)
2588 class(laktype),
intent(inout) :: this
2589 integer(I4B),
intent(in) :: ilak
2590 real(DP),
intent(in) :: stage
2591 real(DP),
intent(inout) :: ra
2593 integer(I4B) :: iconn
2597 iconn = this%idxlakeconn(ilak)
2598 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
2599 sa = this%sareamax(ilak)
2601 call this%lak_calculate_sarea(ilak, stage, sa)
2603 ra = this%rainfall(ilak) * sa
2604 end subroutine lak_calculate_rainfall
2608 subroutine lak_calculate_runoff(this, ilak, ro)
2610 class(laktype),
intent(inout) :: this
2611 integer(I4B),
intent(in) :: ilak
2612 real(DP),
intent(inout) :: ro
2615 ro = this%runoff(ilak)
2616 end subroutine lak_calculate_runoff
2620 subroutine lak_calculate_inflow(this, ilak, qin)
2622 class(laktype),
intent(inout) :: this
2623 integer(I4B),
intent(in) :: ilak
2624 real(DP),
intent(inout) :: qin
2627 qin = this%inflow(ilak)
2628 end subroutine lak_calculate_inflow
2632 subroutine lak_calculate_external(this, ilak, ex)
2634 class(laktype),
intent(inout) :: this
2635 integer(I4B),
intent(in) :: ilak
2636 real(DP),
intent(inout) :: ex
2641 if (this%imover == 1)
then
2642 ex = this%pakmvrobj%get_qfrommvr(ilak)
2644 end subroutine lak_calculate_external
2648 subroutine lak_calculate_withdrawal(this, ilak, avail, wr)
2650 class(laktype),
intent(inout) :: this
2651 integer(I4B),
intent(in) :: ilak
2652 real(DP),
intent(inout) :: avail
2653 real(DP),
intent(inout) :: wr
2656 wr = this%withdrawal(ilak)
2657 if (wr > avail)
then
2660 if (wr >
dzero)
then
2665 end subroutine lak_calculate_withdrawal
2670 subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
2672 class(laktype),
intent(inout) :: this
2673 integer(I4B),
intent(in) :: ilak
2674 real(DP),
intent(in) :: stage
2675 real(DP),
intent(inout) :: avail
2676 real(DP),
intent(inout) :: ev
2681 call this%lak_calculate_sarea(ilak, stage, sa)
2682 ev = sa * this%evaporation(ilak)
2683 if (ev > avail)
then
2693 end subroutine lak_calculate_evaporation
2697 subroutine lak_calculate_outlet_inflow(this, ilak, outinf)
2699 class(laktype),
intent(inout) :: this
2700 integer(I4B),
intent(in) :: ilak
2701 real(DP),
intent(inout) :: outinf
2706 do n = 1, this%noutlets
2707 if (this%lakeout(n) == ilak)
then
2708 outinf = outinf - this%simoutrate(n)
2709 if (this%imover == 1)
then
2710 outinf = outinf - this%pakmvrobj%get_qtomvr(n)
2714 end subroutine lak_calculate_outlet_inflow
2718 subroutine lak_calculate_outlet_outflow(this, ilak, stage, avail, outoutf)
2720 class(laktype),
intent(inout) :: this
2721 integer(I4B),
intent(in) :: ilak
2722 real(DP),
intent(in) :: stage
2723 real(DP),
intent(inout) :: avail
2724 real(DP),
intent(inout) :: outoutf
2734 do n = 1, this%noutlets
2735 if (this%lakein(n) == ilak)
then
2737 d = stage - this%outinvert(n)
2738 if (this%outdmax >
dzero)
then
2739 if (d > this%outdmax) d = this%outdmax
2741 g =
dgravity * this%convlength * this%convtime * this%convtime
2742 select case (this%iouttype(n))
2745 rate = this%outrate(n)
2746 if (-rate > avail)
then
2752 c = (this%convlength**
donethird) * this%convtime
2754 if (this%outrough(n) >
dzero)
then
2755 gsm =
done / this%outrough(n)
2757 rate = -c * gsm * this%outwidth(n) * (d**
dfivethirds) * &
2758 sqrt(this%outslope(n))
2767 this%simoutrate(n) = rate
2768 avail = avail + rate
2769 outoutf = outoutf + rate
2772 end subroutine lak_calculate_outlet_outflow
2782 subroutine lak_outlet_outflow_rate(this, ilak, stage, qout)
2784 class(laktype),
intent(inout) :: this
2785 integer(I4B),
intent(in) :: ilak
2786 real(DP),
intent(in) :: stage
2787 real(DP),
intent(inout) :: qout
2790 real(DP) :: g, d, c, gsm, rate
2793 do n = 1, this%noutlets
2794 if (this%lakein(n) /= ilak) cycle
2796 d = stage - this%outinvert(n)
2797 if (this%outdmax >
dzero .and. d > this%outdmax) d = this%outdmax
2798 g =
dgravity * this%convlength * this%convtime * this%convtime
2799 select case (this%iouttype(n))
2801 rate = this%outrate(n)
2804 c = (this%convlength**
donethird) * this%convtime
2806 if (this%outrough(n) >
dzero) gsm =
done / this%outrough(n)
2807 rate = -c * gsm * this%outwidth(n) * (d**
dfivethirds) * &
2808 sqrt(this%outslope(n))
2817 end subroutine lak_outlet_outflow_rate
2821 subroutine lak_get_internal_inlet(this, ilak, outinf)
2823 class(laktype),
intent(inout) :: this
2824 integer(I4B),
intent(in) :: ilak
2825 real(DP),
intent(inout) :: outinf
2830 do n = 1, this%noutlets
2831 if (this%lakeout(n) == ilak)
then
2832 outinf = outinf - this%simoutrate(n)
2833 if (this%imover == 1)
then
2834 outinf = outinf - this%pakmvrobj%get_qtomvr(n)
2838 end subroutine lak_get_internal_inlet
2842 subroutine lak_get_internal_outlet(this, ilak, outoutf)
2844 class(laktype),
intent(inout) :: this
2845 integer(I4B),
intent(in) :: ilak
2846 real(DP),
intent(inout) :: outoutf
2851 do n = 1, this%noutlets
2852 if (this%lakein(n) == ilak)
then
2853 if (this%lakeout(n) < 1) cycle
2854 outoutf = outoutf + this%simoutrate(n)
2857 end subroutine lak_get_internal_outlet
2861 subroutine lak_get_external_outlet(this, ilak, outoutf)
2863 class(laktype),
intent(inout) :: this
2864 integer(I4B),
intent(in) :: ilak
2865 real(DP),
intent(inout) :: outoutf
2870 do n = 1, this%noutlets
2871 if (this%lakein(n) == ilak)
then
2872 if (this%lakeout(n) > 0) cycle
2873 outoutf = outoutf + this%simoutrate(n)
2876 end subroutine lak_get_external_outlet
2880 subroutine lak_get_external_mover(this, ilak, outoutf)
2882 class(laktype),
intent(inout) :: this
2883 integer(I4B),
intent(in) :: ilak
2884 real(DP),
intent(inout) :: outoutf
2889 if (this%imover == 1)
then
2890 do n = 1, this%noutlets
2891 if (this%lakein(n) == ilak)
then
2892 if (this%lakeout(n) > 0) cycle
2893 outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2897 end subroutine lak_get_external_mover
2901 subroutine lak_get_internal_mover(this, ilak, outoutf)
2903 class(laktype),
intent(inout) :: this
2904 integer(I4B),
intent(in) :: ilak
2905 real(DP),
intent(inout) :: outoutf
2910 if (this%imover == 1)
then
2911 do n = 1, this%noutlets
2912 if (this%lakein(n) == ilak)
then
2913 if (this%lakeout(n) < 1) cycle
2914 outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2918 end subroutine lak_get_internal_mover
2922 subroutine lak_get_outlet_tomover(this, ilak, outoutf)
2924 class(laktype),
intent(inout) :: this
2925 integer(I4B),
intent(in) :: ilak
2926 real(DP),
intent(inout) :: outoutf
2931 if (this%imover == 1)
then
2932 do n = 1, this%noutlets
2933 if (this%lakein(n) == ilak)
then
2934 outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2938 end subroutine lak_get_outlet_tomover
2942 subroutine lak_vol2stage(this, ilak, vol, stage)
2944 class(laktype),
intent(inout) :: this
2945 integer(I4B),
intent(in) :: ilak
2946 real(DP),
intent(in) :: vol
2947 real(DP),
intent(inout) :: stage
2951 real(DP) :: s0, s1, sm
2952 real(DP) :: v0, v1, vm
2953 real(DP) :: f0, f1, fm
2955 real(DP) :: en0, en1
2959 s0 = this%lakebot(ilak)
2960 call this%lak_calculate_vol(ilak, s0, v0)
2961 s1 = this%laketop(ilak)
2962 call this%lak_calculate_vol(ilak, s1, v1)
2967 else if (vol >= v1)
then
2968 call this%lak_calculate_sarea(ilak, s1, sa)
2969 stage = s1 + (vol - v1) / sa
2980 secantbisection:
do i = 1, 150
2982 if (denom /=
dzero)
then
2983 ds = f1 * (s1 - s0) / denom
2991 if (sm < en0 .or. sm > en1) ibs = 13
2995 if (ds * ds0 <
dprec .or. abs(ds) > abs(ds0)) ibs = ibs + 1
2997 ds =
dhalf * (s1 - s0)
3001 if (abs(ds) <
dem6)
then
3002 exit secantbisection
3004 call this%lak_calculate_vol(ilak, sm, vm)
3011 end do secantbisection
3013 if (abs(ds) >=
dem6)
then
3014 write (this%iout,
'(1x,a,1x,i0,4(1x,a,1x,g15.6))') &
3015 &
'LAK_VOL2STAGE failed for lake', ilak,
'volume error =', fm, &
3016 &
'finding stage (', stage,
') for volume =', vol, &
3017 &
'final change in stage =', ds
3020 end subroutine lak_vol2stage
3023 function lak_check_valid(this, itemno)
result(ierr)
3027 integer(I4B) :: ierr
3029 class(laktype),
intent(inout) :: this
3030 integer(I4B),
intent(in) :: itemno
3032 integer(I4B) :: ival
3036 if (itemno > 0)
then
3037 if (ival < 1 .or. ival > this%nlakes)
then
3038 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,a)') &
3039 'LAKENO', itemno,
'must be greater than 0 and less than or equal to', &
3045 if (ival < 1 .or. ival > this%noutlets)
then
3046 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,a)') &
3047 'IOUTLET', itemno,
'must be greater than 0 and less than or equal to', &
3053 end function lak_check_valid
3057 subroutine lak_set_stressperiod(this, itemno)
3062 class(laktype),
intent(inout) :: this
3063 integer(I4B),
intent(in) :: itemno
3065 character(len=LINELENGTH) :: text
3066 character(len=LINELENGTH) :: caux
3067 character(len=LINELENGTH) :: keyword
3068 integer(I4B) :: ierr
3071 real(DP),
pointer :: bndElem => null()
3074 call this%parser%GetStringCaps(keyword)
3075 select case (keyword)
3077 ierr = this%lak_check_valid(itemno)
3081 call this%parser%GetStringCaps(text)
3082 this%status(itemno) = text(1:8)
3083 if (text ==
'CONSTANT')
then
3084 this%iboundpak(itemno) = -1
3085 else if (text ==
'INACTIVE')
then
3086 this%iboundpak(itemno) = 0
3087 else if (text ==
'ACTIVE')
then
3088 this%iboundpak(itemno) = 1
3090 write (
errmsg,
'(a,a)') &
3091 'Unknown '//trim(this%text)//
' lak status keyword: ', text//
'.'
3095 ierr = this%lak_check_valid(itemno)
3099 call this%parser%GetString(text)
3101 bndelem => this%stage(itemno)
3103 this%packName,
'BND', this%tsManager, &
3104 this%iprpak,
'STAGE')
3106 ierr = this%lak_check_valid(itemno)
3110 call this%parser%GetString(text)
3112 bndelem => this%rainfall(itemno)
3114 this%packName,
'BND', this%tsManager, &
3115 this%iprpak,
'RAINFALL')
3116 if (this%rainfall(itemno) <
dzero)
then
3117 write (
errmsg,
'(a,i0,a,G0,a)') &
3118 'Lake ', itemno,
' was assigned a rainfall value of ', &
3119 this%rainfall(itemno),
'. Rainfall must be positive.'
3122 case (
'EVAPORATION')
3123 ierr = this%lak_check_valid(itemno)
3127 call this%parser%GetString(text)
3129 bndelem => this%evaporation(itemno)
3131 this%packName,
'BND', this%tsManager, &
3132 this%iprpak,
'EVAPORATION')
3133 if (this%evaporation(itemno) <
dzero)
then
3134 write (
errmsg,
'(a,i0,a,G0,a)') &
3135 'Lake ', itemno,
' was assigned an evaporation value of ', &
3136 this%evaporation(itemno),
'. Evaporation must be positive.'
3140 ierr = this%lak_check_valid(itemno)
3144 call this%parser%GetString(text)
3146 bndelem => this%runoff(itemno)
3148 this%packName,
'BND', this%tsManager, &
3149 this%iprpak,
'RUNOFF')
3150 if (this%runoff(itemno) <
dzero)
then
3151 write (
errmsg,
'(a,i0,a,G0,a)') &
3152 'Lake ', itemno,
' was assigned a runoff value of ', &
3153 this%runoff(itemno),
'. Runoff must be positive.'
3157 ierr = this%lak_check_valid(itemno)
3161 call this%parser%GetString(text)
3163 bndelem => this%inflow(itemno)
3165 this%packName,
'BND', this%tsManager, &
3166 this%iprpak,
'INFLOW')
3167 if (this%inflow(itemno) <
dzero)
then
3168 write (
errmsg,
'(a,i0,a,G0,a)') &
3169 'Lake ', itemno,
' was assigned an inflow value of ', &
3170 this%inflow(itemno),
'. Inflow must be positive.'
3174 ierr = this%lak_check_valid(itemno)
3178 call this%parser%GetString(text)
3180 bndelem => this%withdrawal(itemno)
3182 this%packName,
'BND', this%tsManager, &
3183 this%iprpak,
'WITHDRAWAL')
3184 if (this%withdrawal(itemno) <
dzero)
then
3185 write (
errmsg,
'(a,i0,a,G0,a)') &
3186 'Lake ', itemno,
' was assigned a withdrawal value of ', &
3187 this%withdrawal(itemno),
'. Withdrawal must be positive.'
3191 ierr = this%lak_check_valid(-itemno)
3195 call this%parser%GetString(text)
3197 bndelem => this%outrate(itemno)
3199 this%packName,
'BND', this%tsManager, &
3200 this%iprpak,
'RATE')
3202 ierr = this%lak_check_valid(-itemno)
3206 call this%parser%GetString(text)
3208 bndelem => this%outinvert(itemno)
3210 this%packName,
'BND', this%tsManager, &
3211 this%iprpak,
'INVERT')
3213 ierr = this%lak_check_valid(-itemno)
3217 call this%parser%GetString(text)
3219 bndelem => this%outwidth(itemno)
3221 this%packName,
'BND', this%tsManager, &
3222 this%iprpak,
'WIDTH')
3224 ierr = this%lak_check_valid(-itemno)
3228 call this%parser%GetString(text)
3230 bndelem => this%outrough(itemno)
3232 this%packName,
'BND', this%tsManager, &
3233 this%iprpak,
'ROUGH')
3235 ierr = this%lak_check_valid(-itemno)
3239 call this%parser%GetString(text)
3241 bndelem => this%outslope(itemno)
3243 this%packName,
'BND', this%tsManager, &
3244 this%iprpak,
'SLOPE')
3246 ierr = this%lak_check_valid(itemno)
3250 call this%parser%GetStringCaps(caux)
3251 do jj = 1, this%naux
3252 if (trim(adjustl(caux)) /= trim(adjustl(this%auxname(jj)))) cycle
3253 call this%parser%GetString(text)
3255 bndelem => this%lauxvar(jj, ii)
3257 this%packName,
'AUX', &
3258 this%tsManager, this%iprpak, &
3264 'Unknown '//trim(this%text)//
' lak data keyword: ', &
3270 end subroutine lak_set_stressperiod
3276 subroutine lak_set_attribute_error(this, ilak, keyword, msg)
3280 class(laktype),
intent(inout) :: this
3281 integer(I4B),
intent(in) :: ilak
3282 character(len=*),
intent(in) :: keyword
3283 character(len=*),
intent(in) :: msg
3285 if (len(msg) == 0)
then
3286 write (
errmsg,
'(a,1x,a,1x,i0,1x,a)') &
3287 keyword,
' for LAKE', ilak,
'has already been set.'
3289 write (
errmsg,
'(a,1x,a,1x,i0,1x,a)') keyword,
' for LAKE', ilak, msg
3292 end subroutine lak_set_attribute_error
3298 subroutine lak_options(this, option, found)
3305 class(laktype),
intent(inout) :: this
3306 character(len=*),
intent(inout) :: option
3307 logical(LGP),
intent(inout) :: found
3309 character(len=MAXCHARLEN) :: fname, keyword
3312 character(len=*),
parameter :: fmtlengthconv = &
3313 &
"(4x, 'LENGTH CONVERSION VALUE (',g15.7,') SPECIFIED.')"
3314 character(len=*),
parameter :: fmttimeconv = &
3315 &
"(4x, 'TIME CONVERSION VALUE (',g15.7,') SPECIFIED.')"
3316 character(len=*),
parameter :: fmtoutdmax = &
3317 &
"(4x, 'MAXIMUM OUTLET WATER DEPTH (',g15.7,') SPECIFIED.')"
3318 character(len=*),
parameter :: fmtlakeopt = &
3319 &
"(4x, 'LAKE ', a, ' VALUE (',g15.7,') SPECIFIED.')"
3320 character(len=*),
parameter :: fmtlakbin = &
3321 "(4x, 'LAK ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', &
3322 &a, /4x, 'OPENED ON UNIT: ', I0)"
3323 character(len=*),
parameter :: fmtiter = &
3324 &
"(4x, 'MAXIMUM LAK ITERATION VALUE (',i0,') SPECIFIED.')"
3325 character(len=*),
parameter :: fmtdmaxchg = &
3326 &
"(4x, 'MAXIMUM STAGE CHANGE VALUE (',g0,') SPECIFIED.')"
3329 select case (option)
3330 case (
'PRINT_STAGE')
3332 write (this%iout,
'(4x,a)') trim(adjustl(this%text))// &
3333 ' STAGES WILL BE PRINTED TO LISTING FILE.'
3335 call this%parser%GetStringCaps(keyword)
3336 if (keyword ==
'FILEOUT')
then
3337 call this%parser%GetString(fname)
3339 call openfile(this%istageout, this%iout, fname,
'DATA(BINARY)', &
3341 write (this%iout, fmtlakbin)
'STAGE', trim(adjustl(fname)), &
3344 call store_error(
'OPTIONAL STAGE KEYWORD MUST BE FOLLOWED BY FILEOUT')
3347 call this%parser%GetStringCaps(keyword)
3348 if (keyword ==
'FILEOUT')
then
3349 call this%parser%GetString(fname)
3350 call assign_iounit(this%ibudgetout, this%inunit,
"BUDGET fileout")
3351 call openfile(this%ibudgetout, this%iout, fname,
'DATA(BINARY)', &
3353 write (this%iout, fmtlakbin)
'BUDGET', trim(adjustl(fname)), &
3356 call store_error(
'OPTIONAL BUDGET KEYWORD MUST BE FOLLOWED BY FILEOUT')
3359 call this%parser%GetStringCaps(keyword)
3360 if (keyword ==
'FILEOUT')
then
3361 call this%parser%GetString(fname)
3362 call assign_iounit(this%ibudcsv, this%inunit,
"BUDGETCSV fileout")
3363 call openfile(this%ibudcsv, this%iout, fname,
'CSV', &
3364 filstat_opt=
'REPLACE')
3365 write (this%iout, fmtlakbin)
'BUDGET CSV', trim(adjustl(fname)), &
3368 call store_error(
'OPTIONAL BUDGETCSV KEYWORD MUST BE FOLLOWED BY &
3371 case (
'PACKAGE_CONVERGENCE')
3372 call this%parser%GetStringCaps(keyword)
3373 if (keyword ==
'FILEOUT')
then
3374 call this%parser%GetString(fname)
3378 this%pakcsvfile = trim(adjustl(fname))
3380 call store_error(
'OPTIONAL PACKAGE_CONVERGENCE KEYWORD MUST BE '// &
3381 'FOLLOWED BY FILEOUT')
3385 write (this%iout,
'(4x,A)')
'MOVER OPTION ENABLED'
3386 case (
'LENGTH_CONVERSION')
3387 this%convlength = this%parser%GetDouble()
3388 write (this%iout, fmtlengthconv) this%convlength
3389 case (
'TIME_CONVERSION')
3390 this%convtime = this%parser%GetDouble()
3391 write (this%iout, fmttimeconv) this%convtime
3393 r = this%parser%GetDouble()
3398 write (this%iout, fmtlakeopt)
'SURFDEP', this%surfdep
3399 case (
'MAXIMUM_ITERATIONS')
3400 this%maxlakit = this%parser%GetInteger()
3401 write (this%iout, fmtiter) this%maxlakit
3402 case (
'MAXIMUM_STAGE_CHANGE')
3403 r = this%parser%GetDouble()
3405 this%delh = dp999 * r
3406 write (this%iout, fmtdmaxchg) this%dmaxchg
3412 case (
'DEV_GROUNDWATER_HEAD_CONDUCTANCE')
3413 call this%parser%DevOpt()
3415 write (this%iout,
'(4x,a)') &
3416 'CONDUCTANCE FOR HORIZONTAL CONNECTIONS WILL BE CALCULATED &
3417 &USING THE GROUNDWATER HEAD'
3418 case (
'DEV_MAXIMUM_OUTLET_DEPTH')
3419 call this%parser%DevOpt()
3420 this%outdmax = this%parser%GetDouble()
3421 write (this%iout, fmtoutdmax) this%outdmax
3424 write (this%iout,
'(4x,a)') &
3425 'LAKE STAGE WILL BE SOLVED AS AN UNKNOWN IN THE GROUNDWATER FLOW '// &
3426 'MATRIX (IMPLICIT FORMULATION)'
3427 case (
'DEV_FORCE_LEGACY')
3428 call this%parser%DevOpt()
3430 write (this%iout,
'(4x,a)') &
3431 'EVERY ACTIVE LAKE WILL BE SOLVED WITH THE LEGACY SUBSTITUTION '// &
3432 'SOLVER UNDER THE IMPLICIT FORMULATION'
3433 if (this%iforceleglak /= 0)
then
3434 call store_error(
'DEV_FORCE_LEGACY and DEV_FORCE_LEGACY_LAKE '// &
3435 'cannot both be specified.')
3437 case (
'DEV_FORCE_LEGACY_LAKE')
3438 call this%parser%DevOpt()
3439 this%iforceleglak = this%parser%GetInteger()
3440 write (this%iout,
'(4x,a,i0,a)')
'LAKE ', this%iforceleglak, &
3441 ' WILL BE SOLVED WITH THE LEGACY SUBSTITUTION SOLVER UNDER THE '// &
3442 'IMPLICIT FORMULATION'
3443 if (this%iforceleg /= 0)
then
3444 call store_error(
'DEV_FORCE_LEGACY and DEV_FORCE_LEGACY_LAKE '// &
3445 'cannot both be specified.')
3447 case (
'DEV_NO_FINAL_CHECK')
3448 call this%parser%DevOpt()
3450 write (this%iout,
'(4x,a)') &
3451 'A FINAL CONVERGENCE CHECK OF THE CHANGE IN LAKE STAGES &
3458 end subroutine lak_options
3464 subroutine lak_ar(this)
3469 class(laktype),
intent(inout) :: this
3471 character(len=*),
parameter :: fmtlakbin = &
3472 "(4x, 'LAK ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', &
3473 &a, /4x, 'OPENED ON UNIT: ', I0)"
3480 if (
allocated(this%pakcsvfile))
then
3481 if (this%iimplicit /= 0)
then
3483 'PACKAGE_CONVERGENCE output file "'//trim(this%pakcsvfile)// &
3484 '" is not written when the IMPLICIT option is active; the lake '// &
3485 'stage is part of the solver (IMS) convergence check.'
3489 call openfile(this%ipakcsv, this%iout, this%pakcsvfile,
'CSV', &
3490 filstat_opt=
'REPLACE', mode_opt=
mnormal)
3491 write (this%iout, fmtlakbin)
'PACKAGE_CONVERGENCE', &
3492 trim(this%pakcsvfile), this%ipakcsv
3494 deallocate (this%pakcsvfile)
3497 call this%obs%obs_ar()
3500 call this%lak_allocate_arrays()
3503 call this%read_initial_attr()
3506 if (this%imover /= 0)
then
3507 allocate (this%pakmvrobj)
3508 call this%pakmvrobj%ar(this%noutlets, this%nlakes, this%memoryPath)
3510 end subroutine lak_ar
3516 subroutine lak_rp(this)
3522 class(laktype),
intent(inout) :: this
3524 character(len=LINELENGTH) :: title
3525 character(len=LINELENGTH) :: line
3526 character(len=LINELENGTH) :: text
3527 logical(LGP) :: isfound
3528 logical(LGP) :: endOfBlock
3529 integer(I4B) :: ierr
3530 integer(I4B) :: node
3532 integer(I4B) :: itemno
3535 character(len=*),
parameter :: fmtblkerr = &
3536 &
"('Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
3537 character(len=*),
parameter :: fmtlsp = &
3538 &
"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
3541 this%nbound = this%maxbound
3545 if (this%inunit == 0)
return
3548 if (this%ionper <
kper)
then
3551 call this%parser%GetBlock(
'PERIOD', isfound, ierr, &
3552 supportopenclose=.true., &
3553 blockrequired=.false.)
3557 call this%read_check_ionper()
3563 this%ionper =
nper + 1
3566 call this%parser%GetCurrentLine(line)
3567 write (
errmsg, fmtblkerr) adjustl(trim(line))
3569 call this%parser%StoreErrorUnit()
3575 if (this%ionper ==
kper)
then
3578 if (this%iprpak /= 0)
then
3581 title = trim(adjustl(this%text))//
' PACKAGE ('// &
3582 trim(adjustl(this%packName))//
') DATA FOR PERIOD'
3583 write (title,
'(a,1x,i6)') trim(adjustl(title)),
kper
3584 call table_cr(this%inputtab, this%packName, title)
3585 call this%inputtab%table_df(1, 4, this%iout, finalize=.false.)
3587 call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
3589 call this%inputtab%initialize_column(text, 20, alignment=tableft)
3591 write (text,
'(a,1x,i6)')
'VALUE', n
3592 call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
3599 call this%parser%GetNextLine(endofblock)
3600 if (endofblock)
exit
3603 itemno = this%parser%GetInteger()
3606 call this%lak_set_stressperiod(itemno)
3609 if (this%iprpak /= 0)
then
3610 call this%parser%GetCurrentLine(line)
3611 call this%inputtab%line_to_columns(line)
3615 if (this%iprpak /= 0)
then
3616 call this%inputtab%finalize_table()
3621 write (this%iout, fmtlsp) trim(this%filtyp)
3626 call this%parser%StoreErrorUnit()
3630 do n = 1, this%nlakes
3631 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3632 node = this%cellid(j)
3633 this%nodelist(j) = node
3634 this%bound(1, j) = this%xnewpak(n)
3635 this%bound(2, j) = this%satcond(j)
3636 this%bound(3, j) = this%belev(j)
3641 if (this%imover == 1)
then
3642 do n = 1, this%noutlets
3643 this%pakmvrobj%iprmap(n) = this%lakein(n)
3646 end subroutine lak_rp
3650 subroutine lak_ad(this)
3654 class(laktype) :: this
3658 integer(I4B) :: iaux
3661 call this%TsManager%ad()
3666 if (this%naux > 0)
then
3667 do n = 1, this%nlakes
3668 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3669 do iaux = 1, this%naux
3670 if (this%noupdateauxvar(iaux) /= 0) cycle
3671 this%auxvar(iaux, j) = this%lauxvar(iaux, n)
3682 do n = 1, this%nlakes
3683 this%xoldpak(n) = this%xnewpak(n)
3684 this%stageiter(n) = this%xnewpak(n)
3685 if (this%iboundpak(n) < 0)
then
3686 this%xnewpak(n) = this%stage(n)
3688 this%seep0(n) =
dzero
3694 do n = 1, this%nlakes
3695 this%xnewpak(n) = this%xoldpak(n)
3696 this%stageiter(n) = this%xnewpak(n)
3697 if (this%iboundpak(n) < 0)
then
3698 this%xnewpak(n) = this%stage(n)
3700 this%seep0(n) =
dzero
3705 if (this%imover == 1)
then
3706 call this%pakmvrobj%ad()
3712 call this%obs%obs_ad()
3713 end subroutine lak_ad
3719 subroutine lak_cf(this)
3721 class(laktype) :: this
3723 integer(I4B) :: j, n
3724 integer(I4B) :: igwfnode
3725 real(DP) :: hlak, bottom_lake
3728 do n = 1, this%nlakes
3729 this%seep0(n) = this%seep(n)
3733 do n = 1, this%nlakes
3734 this%s0(n) = this%xnewpak(n)
3735 call this%lak_calculate_exchange(n, this%s0(n), this%qgwf0(n))
3739 do n = 1, this%nlakes
3740 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3742 if (this%ictype(j) /= 0)
then
3745 igwfnode = this%nodesontop(j)
3746 if (this%ibound(igwfnode) == 0)
then
3747 call this%dis%highest_active(igwfnode, this%ibound)
3749 this%nodelist(j) = igwfnode
3750 this%cellid(j) = igwfnode
3757 do n = 1, this%nlakes
3759 hlak = this%xnewpak(n)
3762 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3765 igwfnode = this%cellid(j)
3768 if (this%ibound(igwfnode) < 1)
then
3773 if (this%ictype(j) /= 0)
then
3778 if (this%ictype(j) == 2 .or. this%ictype(j) == 3)
then
3783 bottom_lake = this%belev(j)
3784 if (hlak > bottom_lake .or. this%iboundpak(n) == 0)
then
3787 this%ibound(igwfnode) = 1
3795 call this%lak_bound_update()
3796 end subroutine lak_cf
3800 subroutine lak_fc(this, rhs, ia, idxglo, matrix_sln)
3802 class(laktype) :: this
3803 real(DP),
dimension(:),
intent(inout) :: rhs
3804 integer(I4B),
dimension(:),
intent(in) :: ia
3805 integer(I4B),
dimension(:),
intent(in) :: idxglo
3808 integer(I4B) :: j, n
3809 integer(I4B) :: igwfnode
3810 integer(I4B) :: ipossymd
3813 if (this%imover == 1)
then
3814 call this%pakmvrobj%fc()
3819 if (this%iimplicit /= 0)
then
3820 call this%lak_fc_implicit(rhs, matrix_sln)
3826 call this%lak_solve()
3827 do n = 1, this%nlakes
3828 if (this%iboundpak(n) == 0) cycle
3829 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3830 igwfnode = this%cellid(j)
3831 if (this%ibound(igwfnode) < 1) cycle
3832 ipossymd = idxglo(ia(igwfnode))
3833 call matrix_sln%add_value_pos(ipossymd, this%hcof(j))
3834 rhs(igwfnode) = rhs(igwfnode) + this%rhs(j)
3837 end subroutine lak_fc
3841 subroutine lak_fn(this, rhs, ia, idxglo, matrix_sln)
3843 class(laktype) :: this
3844 real(DP),
dimension(:),
intent(inout) :: rhs
3845 integer(I4B),
dimension(:),
intent(in) :: ia
3846 integer(I4B),
dimension(:),
intent(in) :: idxglo
3849 integer(I4B) :: j, n
3850 integer(I4B) :: ipos
3851 integer(I4B) :: igwfnode
3852 integer(I4B) :: idry
3869 if (this%iimplicit /= 0)
then
3873 do n = 1, this%nlakes
3874 if (this%iboundpak(n) == 0) cycle
3875 hlak = this%xnewpak(n)
3876 call this%lak_calculate_available(n, hlak, avail, &
3877 ra, ro, qinf, ex, this%delh)
3878 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3879 igwfnode = this%cellid(j)
3881 head = this%xnew(igwfnode)
3882 if (-this%hcof(j) >
dzero)
then
3883 if (this%ibound(igwfnode) > 0)
then
3887 call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, &
3888 head + this%delh, q1, avail)
3891 q = this%hcof(j) * head - this%rhs(j)
3893 rterm = this%hcof(j) * head
3895 drterm = (q1 - q) / this%delh
3898 call matrix_sln%add_value_pos(idxglo(ipos), drterm - this%hcof(j))
3899 rhs(igwfnode) = rhs(igwfnode) - rterm + drterm * head
3904 end subroutine lak_fn
3914 subroutine lak_nur(this, neqpak, x, xtemp, dx, inewtonur, dxmax, locmax)
3916 class(laktype),
intent(inout) :: this
3917 integer(I4B),
intent(in) :: neqpak
3918 real(DP),
dimension(neqpak),
intent(inout) :: x
3919 real(DP),
dimension(neqpak),
intent(in) :: xtemp
3920 real(DP),
dimension(neqpak),
intent(inout) :: dx
3921 integer(I4B),
intent(inout) :: inewtonur
3922 real(DP),
intent(inout) :: dxmax
3923 integer(I4B),
intent(inout) :: locmax
3931 if (this%iimplicit == 0)
return
3934 do n = 1, this%nlakes
3935 if (this%iboundpak(n) < 1) cycle
3936 botl = this%lakebot(n)
3940 if (x(n) < botl)
then
3944 if (abs(dxx) > abs(dxmax))
then
3952 end subroutine lak_nur
3956 subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
3960 class(laktype),
intent(inout) :: this
3961 integer(I4B),
intent(in) :: innertot
3962 integer(I4B),
intent(in) :: kiter
3963 integer(I4B),
intent(in) :: iend
3964 integer(I4B),
intent(in) :: icnvgmod
3965 character(len=LENPAKLOC),
intent(inout) :: cpak
3966 integer(I4B),
intent(inout) :: ipak
3967 real(DP),
intent(inout) :: dpak
3969 character(len=LENPAKLOC) :: cloc
3970 character(len=LINELENGTH) :: tag
3971 integer(I4B) :: icheck
3972 integer(I4B) :: ipakfail
3973 integer(I4B) :: locdhmax
3974 integer(I4B) :: locresidmax
3975 integer(I4B) :: locdgwfmax
3976 integer(I4B) :: locdqoutmax
3977 integer(I4B) :: locdqfrommvrmax
3978 integer(I4B) :: ntabrows
3979 integer(I4B) :: ntabcols
3983 real(DP) :: qtolfact
4001 real(DP) :: residmax
4003 real(DP) :: dqoutmax
4004 real(DP) :: dqfrommvr
4005 real(DP) :: dqfrommvrmax
4007 call this%lak_set_legacy(kiter, icnvgmod)
4012 if (iend /= 0 .and. icnvgmod == 0)
then
4013 call this%lak_check_disconnected()
4019 if (this%iimplicit /= 0)
then
4024 icheck = this%iconvchk
4035 dqfrommvrmax =
dzero
4039 if (this%ipakcsv == 0)
then
4040 if (icnvgmod == 0)
then
4048 if (.not.
associated(this%pakcsvtab))
then
4053 if (this%noutlets > 0)
then
4054 ntabcols = ntabcols + 2
4056 if (this%imover == 1)
then
4057 ntabcols = ntabcols + 2
4061 call table_cr(this%pakcsvtab, this%packName,
'')
4062 call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
4063 lineseparator=.false., separator=
',', &
4067 tag =
'total_inner_iterations'
4068 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4070 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4072 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4074 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4076 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4078 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4080 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4082 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4083 tag =
'residmax_loc'
4084 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4086 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4088 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4089 if (this%noutlets > 0)
then
4091 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4092 tag =
'dqoutmax_loc'
4093 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4095 if (this%imover == 1)
then
4096 tag =
'dqfrommvrmax'
4097 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4098 tag =
'dqfrommvrmax_loc'
4099 call this%pakcsvtab%initialize_column(tag, 16, alignment=
tableft)
4105 if (icheck /= 0)
then
4106 final_check:
do n = 1, this%nlakes
4107 if (this%iboundpak(n) < 1) cycle
4111 hlak = this%xnewpak(n)
4117 call this%lak_calculate_sarea(n, hlak, area)
4120 if (area >
dzero)
then
4121 qtolfact =
delt / area
4127 call this%lak_calculate_residual(n, hlak, resid)
4128 resid = resid * qtolfact
4132 if (area >
dzero)
then
4133 gwf0 = this%qgwf0(n)
4134 call this%lak_calculate_exchange(n, hlak, gwf)
4135 dgwf = (gwf0 - gwf) * qtolfact
4140 if (this%noutlets > 0)
then
4141 if (area >
dzero)
then
4142 call this%lak_calculate_available(n, hlak0, inf, ra, ro, qinf, ex)
4143 call this%lak_calculate_outlet_outflow(n, hlak0, inf, qout0)
4144 call this%lak_calculate_available(n, hlak, inf, ra, ro, qinf, ex)
4145 call this%lak_calculate_outlet_outflow(n, hlak, inf, qout)
4146 dqout = (qout0 - qout) * qtolfact
4152 if (this%imover == 1)
then
4153 q = this%pakmvrobj%get_qfrommvr(n)
4154 q0 = this%pakmvrobj%get_qfrommvr0(n)
4155 dqfrommvr = qtolfact * (q0 - q)
4168 dqfrommvrmax = dqfrommvr
4171 if (abs(dh) > abs(dhmax))
then
4175 if (abs(resid) > abs(residmax))
then
4179 if (abs(dgwf) > abs(dgwfmax))
then
4183 if (abs(dqout) > abs(dqoutmax))
then
4187 if (abs(dqfrommvr) > abs(dqfrommvrmax))
then
4188 dqfrommvrmax = dqfrommvr
4195 if (abs(dhmax) > abs(dpak))
then
4198 write (cloc,
"(a,'-',a)") &
4199 trim(this%packName),
'stage'
4202 if (abs(residmax) > abs(dpak))
then
4205 write (cloc,
"(a,'-',a)") &
4206 trim(this%packName),
'residual'
4209 if (abs(dgwfmax) > abs(dpak))
then
4212 write (cloc,
"(a,'-',a)") &
4213 trim(this%packName),
'gwf'
4216 if (this%noutlets > 0)
then
4217 if (abs(dqoutmax) > abs(dpak))
then
4220 write (cloc,
"(a,'-',a)") &
4221 trim(this%packName),
'outlet'
4225 if (this%imover == 1)
then
4226 if (abs(dqfrommvrmax) > abs(dpak))
then
4227 ipak = locdqfrommvrmax
4229 write (cloc,
"(a,'-',a)") trim(this%packName),
'qfrommvr'
4235 if (this%ipakcsv /= 0)
then
4238 call this%pakcsvtab%add_term(innertot)
4239 call this%pakcsvtab%add_term(
totim)
4240 call this%pakcsvtab%add_term(
kper)
4241 call this%pakcsvtab%add_term(
kstp)
4242 call this%pakcsvtab%add_term(kiter)
4243 call this%pakcsvtab%add_term(dhmax)
4244 call this%pakcsvtab%add_term(locdhmax)
4245 call this%pakcsvtab%add_term(residmax)
4246 call this%pakcsvtab%add_term(locresidmax)
4247 call this%pakcsvtab%add_term(dgwfmax)
4248 call this%pakcsvtab%add_term(locdgwfmax)
4249 if (this%noutlets > 0)
then
4250 call this%pakcsvtab%add_term(dqoutmax)
4251 call this%pakcsvtab%add_term(locdqoutmax)
4253 if (this%imover == 1)
then
4254 call this%pakcsvtab%add_term(dqfrommvrmax)
4255 call this%pakcsvtab%add_term(locdqfrommvrmax)
4260 call this%pakcsvtab%finalize_table()
4264 end subroutine lak_cc
4268 subroutine lak_cq(this, x, flowja, iadv)
4272 class(laktype),
intent(inout) :: this
4273 real(DP),
dimension(:),
intent(in) :: x
4274 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
4275 integer(I4B),
optional,
intent(in) :: iadv
4278 real(DP) :: chratin, chratout
4280 integer(I4B) :: j, n, igwfnode
4281 real(DP) :: hlak, head, flow, dqdh
4282 real(DP) :: v0, v1, sa, sf
4284 call this%lak_solve(update=.false.)
4293 if (this%iimplicit /= 0)
then
4294 do n = 1, this%nlakes
4295 if (this%iboundpak(n) < 1 .or. this%ilegacy(n) /= 0) cycle
4296 hlak = this%xnewpak(n)
4299 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4300 igwfnode = this%cellid(j)
4301 if (this%ibound(igwfnode) < 1) cycle
4302 head = this%xnew(igwfnode)
4303 call this%lak_calculate_conn_exchange_deriv(n, j, hlak, head, &
4305 this%hcof(j) = -dqdh
4306 this%rhs(j) = -dqdh * head + flow
4316 call this%lak_calculate_sarea(n, hlak, sa)
4318 if (this%surfdep >
dzero)
then
4320 this%lakebot(n), hlak)
4322 this%evap(n) = -this%evaporation(n) * sa * sf
4323 this%withr(n) = -this%withdrawal(n) * sf
4329 call this%BndType%bnd_cq(x, flowja, iadv=1)
4334 do n = 1, this%nlakes
4335 this%chterm(n) =
dzero
4336 if (this%iboundpak(n) == 0) cycle
4337 hlak = this%xnewpak(n)
4338 call this%lak_calculate_vol(n, hlak, v1)
4341 if (this%iboundpak(n) /= 0)
then
4344 rrate = this%precip(n)
4345 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4348 rrate = this%evap(n)
4349 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4352 rrate = this%runoff(n)
4353 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4356 rrate = this%inflow(n)
4357 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4360 rrate = this%withr(n)
4361 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4365 if (this%iboundpak(n) > 0)
then
4366 if (this%gwfiss /= 1)
then
4367 call this%lak_calculate_vol(n, this%xoldpak(n), v0)
4368 rrate = -(v1 - v0) /
delt
4369 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4372 this%qsto(n) = rrate
4375 call this%lak_get_external_outlet(n, rrate)
4376 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4379 if (this%imover == 1)
then
4380 if (this%iboundpak(n) /= 0)
then
4381 rrate = this%pakmvrobj%get_qfrommvr(n)
4385 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4391 do n = 1, this%nlakes
4393 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4397 rrate = -this%simvals(j)
4398 this%qleak(j) = rrate
4399 if (this%iboundpak(n) /= 0)
then
4400 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4406 call this%lak_fill_budobj()
4407 end subroutine lak_cq
4411 subroutine lak_ot_package_flows(this, icbcfl, ibudfl)
4413 class(laktype) :: this
4414 integer(I4B),
intent(in) :: icbcfl
4415 integer(I4B),
intent(in) :: ibudfl
4416 integer(I4B) :: ibinun
4420 if (this%ibudgetout /= 0)
then
4421 ibinun = this%ibudgetout
4423 if (icbcfl == 0) ibinun = 0
4424 if (ibinun > 0)
then
4425 call this%budobj%save_flows(this%dis, ibinun,
kstp,
kper,
delt, &
4430 if (ibudfl /= 0 .and. this%iprflow /= 0)
then
4431 call this%budobj%write_flowtable(this%dis,
kstp,
kper)
4433 end subroutine lak_ot_package_flows
4437 subroutine lak_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
4438 class(laktype) :: this
4439 integer(I4B),
intent(in) :: icbcfl
4440 integer(I4B),
intent(in) :: ibudfl
4441 integer(I4B),
intent(in) :: icbcun
4442 integer(I4B),
dimension(:),
optional,
intent(in) :: imap
4445 call this%BndType%bnd_ot_model_flows(icbcfl, ibudfl, icbcun, this%imap)
4446 end subroutine lak_ot_model_flows
4450 subroutine lak_ot_dv(this, idvsave, idvprint)
4454 class(laktype) :: this
4455 integer(I4B),
intent(in) :: idvsave
4456 integer(I4B),
intent(in) :: idvprint
4457 integer(I4B) :: ibinun
4467 if (this%istageout /= 0)
then
4468 ibinun = this%istageout
4470 if (idvsave == 0) ibinun = 0
4473 if (ibinun > 0)
then
4474 do n = 1, this%nlakes
4476 d = v - this%lakebot(n)
4477 if (this%iboundpak(n) == 0)
then
4479 else if (d <= dzero)
then
4485 this%nlakes, 1, 1, ibinun)
4489 if (idvprint /= 0 .and. this%iprhed /= 0)
then
4492 call this%stagetab%set_kstpkper(
kstp,
kper)
4495 do n = 1, this%nlakes
4496 if (this%iboundpak(n) == 0)
then
4502 stage = this%xnewpak(n)
4503 call this%lak_calculate_sarea(n, stage, sa)
4504 call this%lak_calculate_warea(n, stage, wa)
4505 call this%lak_calculate_vol(n, stage, v)
4507 if (this%inamedbound == 1)
then
4508 call this%stagetab%add_term(this%lakename(n))
4510 call this%stagetab%add_term(n)
4511 call this%stagetab%add_term(stage)
4512 call this%stagetab%add_term(sa)
4513 call this%stagetab%add_term(wa)
4514 call this%stagetab%add_term(v)
4517 end subroutine lak_ot_dv
4521 subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl)
4525 class(laktype) :: this
4526 integer(I4B),
intent(in) :: kstp
4527 integer(I4B),
intent(in) :: kper
4528 integer(I4B),
intent(in) :: iout
4529 integer(I4B),
intent(in) :: ibudfl
4531 call this%budobj%write_budtable(kstp, kper, iout, ibudfl,
totim,
delt)
4532 end subroutine lak_ot_bdsummary
4536 subroutine lak_da(this)
4540 class(laktype) :: this
4543 deallocate (this%lakename)
4544 deallocate (this%status)
4545 deallocate (this%clakbudget)
4547 deallocate (this%cauxcbc)
4556 if (this%ntables > 0)
then
4565 call this%budobj%budgetobject_da()
4566 deallocate (this%budobj)
4567 nullify (this%budobj)
4570 if (this%noutlets > 0)
then
4583 if (this%iprhed > 0)
then
4584 call this%stagetab%table_da()
4585 deallocate (this%stagetab)
4586 nullify (this%stagetab)
4590 if (this%ipakcsv > 0)
then
4591 if (
associated(this%pakcsvtab))
then
4592 call this%pakcsvtab%table_da()
4593 deallocate (this%pakcsvtab)
4594 nullify (this%pakcsvtab)
4604 if (
allocated(this%pakcsvfile))
deallocate (this%pakcsvfile)
4664 if (this%iimplicit == 0)
then
4709 nullify (this%gwfiss)
4712 call this%BndType%bnd_da()
4713 end subroutine lak_da
4718 subroutine define_listlabel(this)
4720 class(laktype),
intent(inout) :: this
4723 this%listlabel = trim(this%filtyp)//
' NO.'
4724 if (this%dis%ndim == 3)
then
4725 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
4726 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'ROW'
4727 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'COL'
4728 elseif (this%dis%ndim == 2)
then
4729 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
4730 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'CELL2D'
4732 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'NODE'
4734 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'STRESS RATE'
4735 if (this%inamedbound == 1)
then
4736 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'BOUNDARY NAME'
4738 end subroutine define_listlabel
4743 subroutine lak_set_pointers(this, neq, ibound, xnew, xold, flowja)
4747 class(laktype) :: this
4748 integer(I4B),
pointer :: neq
4749 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
4750 real(DP),
dimension(:),
pointer,
contiguous :: xnew
4751 real(DP),
dimension(:),
pointer,
contiguous :: xold
4752 real(DP),
dimension(:),
pointer,
contiguous :: flowja
4755 integer(I4B) :: istart, iend
4758 call this%BndType%set_pointers(neq, ibound, xnew, xold, flowja)
4764 if (this%iimplicit /= 0)
then
4765 istart = this%dis%nodes + this%ioffset + 1
4766 iend = istart + this%nlakes - 1
4767 this%iboundpak => this%ibound(istart:iend)
4768 this%xnewpak => this%xnew(istart:iend)
4769 call mem_checkin(this%xnewpak,
'XNEWPAK', this%memoryPath,
'X', &
4770 this%memoryPathModel)
4773 do n = 1, this%nlakes
4774 this%xnewpak(n) =
dep20
4777 end subroutine lak_set_pointers
4785 subroutine lak_ac(this, moffset, sparse)
4788 class(laktype),
intent(inout) :: this
4789 integer(I4B),
intent(in) :: moffset
4792 integer(I4B) :: j, n
4794 integer(I4B) :: jglo
4795 integer(I4B) :: nglo
4798 if (this%iimplicit == 0)
return
4801 do n = 1, this%nlakes
4802 nglo = moffset + this%dis%nodes + this%ioffset + n
4803 call sparse%addconnection(nglo, nglo, 1)
4804 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4807 call sparse%addconnection(nglo, jglo, 1)
4808 call sparse%addconnection(jglo, nglo, 1)
4811 end subroutine lak_ac
4819 subroutine lak_mc(this, moffset, matrix_sln)
4822 class(laktype),
intent(inout) :: this
4823 integer(I4B),
intent(in) :: moffset
4828 integer(I4B) :: iglo
4829 integer(I4B) :: jglo
4830 integer(I4B) :: ipos
4835 if (this%iimplicit == 0)
then
4836 call mem_allocate(this%idxlocnode, 0,
'IDXLOCNODE', this%memoryPath)
4837 call mem_allocate(this%idxdiag, 0,
'IDXDIAG', this%memoryPath)
4838 call mem_allocate(this%idxoffdglo, 0,
'IDXOFFDGLO', this%memoryPath)
4839 call mem_allocate(this%idxsymdglo, 0,
'IDXSYMDGLO', this%memoryPath)
4840 call mem_allocate(this%idxsymoffdglo, 0,
'IDXSYMOFFDGLO', this%memoryPath)
4843 call mem_allocate(this%idxlocnode, this%nlakes,
'IDXLOCNODE', &
4845 call mem_allocate(this%idxdiag, this%nlakes,
'IDXDIAG', this%memoryPath)
4846 call mem_allocate(this%idxoffdglo, this%maxbound,
'IDXOFFDGLO', &
4848 call mem_allocate(this%idxsymdglo, this%maxbound,
'IDXSYMDGLO', &
4850 call mem_allocate(this%idxsymoffdglo, this%maxbound,
'IDXSYMOFFDGLO', &
4857 do n = 1, this%nlakes
4858 iglo = moffset + this%dis%nodes + this%ioffset + n
4859 this%idxlocnode(n) = this%dis%nodes + this%ioffset + n
4860 this%idxdiag(n) = matrix_sln%get_position_diag(iglo)
4861 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4862 jglo = this%cellid(j) + moffset
4863 this%idxoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
4870 do n = 1, this%nlakes
4871 jglo = moffset + this%dis%nodes + this%ioffset + n
4872 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4873 iglo = this%cellid(j) + moffset
4874 this%idxsymdglo(ipos) = matrix_sln%get_position_diag(iglo)
4875 this%idxsymoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
4879 end subroutine lak_mc
4886 logical function lak_obs_supported(this)
4888 class(laktype) :: this
4890 lak_obs_supported = .true.
4891 end function lak_obs_supported
4896 subroutine lak_df_obs(this)
4898 class(laktype) :: this
4900 integer(I4B) :: indx
4904 call this%obs%StoreObsType(
'stage', .false., indx)
4905 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4909 call this%obs%StoreObsType(
'ext-inflow', .true., indx)
4910 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4914 call this%obs%StoreObsType(
'outlet-inflow', .true., indx)
4915 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4919 call this%obs%StoreObsType(
'inflow', .true., indx)
4920 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4924 call this%obs%StoreObsType(
'from-mvr', .true., indx)
4925 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4929 call this%obs%StoreObsType(
'rainfall', .true., indx)
4930 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4934 call this%obs%StoreObsType(
'runoff', .true., indx)
4935 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4939 call this%obs%StoreObsType(
'lak', .true., indx)
4940 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4944 call this%obs%StoreObsType(
'evaporation', .true., indx)
4945 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4949 call this%obs%StoreObsType(
'withdrawal', .true., indx)
4950 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4954 call this%obs%StoreObsType(
'ext-outflow', .true., indx)
4955 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4959 call this%obs%StoreObsType(
'to-mvr', .true., indx)
4960 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4964 call this%obs%StoreObsType(
'storage', .true., indx)
4965 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4969 call this%obs%StoreObsType(
'constant', .true., indx)
4970 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4974 call this%obs%StoreObsType(
'outlet', .true., indx)
4975 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4979 call this%obs%StoreObsType(
'volume', .true., indx)
4980 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4984 call this%obs%StoreObsType(
'surface-area', .true., indx)
4985 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4989 call this%obs%StoreObsType(
'wetted-area', .true., indx)
4990 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4994 call this%obs%StoreObsType(
'conductance', .true., indx)
4995 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4996 end subroutine lak_df_obs
5001 subroutine lak_bd_obs(this)
5003 class(laktype) :: this
5006 integer(I4B) :: igwfnode
5017 if (this%obs%npakobs > 0)
then
5018 call this%obs%obs_bd_clear()
5019 do i = 1, this%obs%npakobs
5020 obsrv => this%obs%pakobs(i)%obsrv
5021 do j = 1, obsrv%indxbnds_count
5023 jj = obsrv%indxbnds(j)
5024 select case (obsrv%ObsTypeId)
5026 if (this%iboundpak(jj) /= 0)
then
5027 v = this%xnewpak(jj)
5030 if (this%iboundpak(jj) /= 0)
then
5031 call this%lak_calculate_inflow(jj, v)
5033 case (
'OUTLET-INFLOW')
5034 if (this%iboundpak(jj) /= 0)
then
5035 call this%lak_calculate_outlet_inflow(jj, v)
5038 if (this%iboundpak(jj) /= 0)
then
5039 call this%lak_calculate_inflow(jj, v)
5040 call this%lak_calculate_outlet_inflow(jj, v2)
5044 if (this%iboundpak(jj) /= 0)
then
5045 if (this%imover == 1)
then
5046 v = this%pakmvrobj%get_qfrommvr(jj)
5050 if (this%iboundpak(jj) /= 0)
then
5054 if (this%iboundpak(jj) /= 0)
then
5059 if (this%iboundpak(n) /= 0)
then
5060 igwfnode = this%cellid(jj)
5061 hgwf = this%xnew(igwfnode)
5062 if (this%hcof(jj) /=
dzero)
then
5063 v = -(this%hcof(jj) * (this%xnewpak(n) - hgwf))
5068 case (
'EVAPORATION')
5069 if (this%iboundpak(jj) /= 0)
then
5073 if (this%iboundpak(jj) /= 0)
then
5076 case (
'EXT-OUTFLOW')
5078 if (this%iboundpak(n) /= 0)
then
5079 if (this%lakeout(jj) == 0)
then
5080 v = this%simoutrate(jj)
5082 if (this%imover == 1)
then
5083 v = v + this%pakmvrobj%get_qtomvr(jj)
5090 if (this%iboundpak(n) /= 0)
then
5091 if (this%imover == 1)
then
5092 v = this%pakmvrobj%get_qtomvr(jj)
5099 if (this%iboundpak(jj) /= 0)
then
5103 if (this%iboundpak(jj) /= 0)
then
5108 if (this%iboundpak(n) /= 0)
then
5109 v = this%simoutrate(jj)
5112 if (this%iboundpak(jj) /= 0)
then
5113 call this%lak_calculate_vol(jj, this%xnewpak(jj), v)
5115 case (
'SURFACE-AREA')
5116 if (this%iboundpak(jj) /= 0)
then
5117 hlak = this%xnewpak(jj)
5118 call this%lak_calculate_sarea(jj, hlak, v)
5120 case (
'WETTED-AREA')
5122 if (this%iboundpak(n) /= 0)
then
5123 hlak = this%xnewpak(n)
5124 igwfnode = this%cellid(jj)
5125 hgwf = this%xnew(igwfnode)
5126 call this%lak_calculate_conn_warea(n, jj, hlak, hgwf, v)
5128 case (
'CONDUCTANCE')
5130 if (this%iboundpak(n) /= 0)
then
5131 hlak = this%xnewpak(n)
5132 igwfnode = this%cellid(jj)
5133 hgwf = this%xnew(igwfnode)
5134 call this%lak_calculate_conn_conductance(n, jj, hlak, hgwf, v)
5137 errmsg =
'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
5140 call this%obs%SaveOneSimval(obsrv, v)
5149 end subroutine lak_bd_obs
5156 subroutine lak_rp_obs(this)
5159 class(laktype),
intent(inout) :: this
5166 character(len=LENBOUNDNAME) :: bname
5167 logical(LGP) :: jfound
5170 10
format(
'Boundary "', a,
'" for observation "', a, &
5171 '" is invalid in package "', a,
'"')
5177 do i = 1, this%obs%npakobs
5178 obsrv => this%obs%pakobs(i)%obsrv
5181 nn1 = obsrv%NodeNumber
5183 bname = obsrv%FeatureName
5184 if (bname /=
'')
then
5189 if (obsrv%ObsTypeId ==
'LAK' .or. &
5190 obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5191 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5192 do j = 1, this%nlakes
5193 do jj = this%idxlakeconn(j), this%idxlakeconn(j + 1) - 1
5194 if (this%boundname(jj) == bname)
then
5196 call obsrv%AddObsIndex(jj)
5200 else if (obsrv%ObsTypeId ==
'EXT-OUTFLOW' .or. &
5201 obsrv%ObsTypeId ==
'TO-MVR' .or. &
5202 obsrv%ObsTypeId ==
'OUTLET')
then
5203 do j = 1, this%noutlets
5205 if (this%lakename(jj) == bname)
then
5207 call obsrv%AddObsIndex(j)
5211 do j = 1, this%nlakes
5212 if (this%lakename(j) == bname)
then
5214 call obsrv%AddObsIndex(j)
5218 if (.not. jfound)
then
5220 trim(bname), trim(obsrv%Name), trim(this%packName)
5225 if (obsrv%indxbnds_count == 0)
then
5226 if (obsrv%ObsTypeId ==
'LAK' .or. &
5227 obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5228 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5229 nn2 = obsrv%NodeNumber2
5230 j = this%idxlakeconn(nn1) + nn2 - 1
5231 call obsrv%AddObsIndex(j)
5233 call obsrv%AddObsIndex(nn1)
5236 errmsg =
'Programming error in lak_rp_obs'
5243 if (obsrv%ObsTypeId ==
'STAGE')
then
5244 if (obsrv%indxbnds_count > 1)
then
5245 write (
errmsg,
'(a,3(1x,a))') &
5246 trim(adjustl(obsrv%ObsTypeId)), &
5247 'for observation', trim(adjustl(obsrv%Name)), &
5248 ' must be assigned to a lake with a unique boundname.'
5254 if (obsrv%ObsTypeId ==
'TO-MVR' .or. &
5255 obsrv%ObsTypeId ==
'EXT-OUTFLOW' .or. &
5256 obsrv%ObsTypeId ==
'OUTLET')
then
5257 do j = 1, obsrv%indxbnds_count
5258 nn1 = obsrv%indxbnds(j)
5259 if (nn1 < 1 .or. nn1 > this%noutlets)
then
5260 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
5261 trim(adjustl(obsrv%ObsTypeId)), &
5262 ' outlet must be > 0 and <=', this%noutlets, &
5263 '(specified value is ', nn1,
')'
5267 else if (obsrv%ObsTypeId ==
'LAK' .or. &
5268 obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5269 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5270 do j = 1, obsrv%indxbnds_count
5271 nn1 = obsrv%indxbnds(j)
5272 if (nn1 < 1 .or. nn1 > this%maxbound)
then
5273 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
5274 trim(adjustl(obsrv%ObsTypeId)), &
5275 'lake connection number must be > 0 and <=', this%maxbound, &
5276 '(specified value is ', nn1,
')'
5281 do j = 1, obsrv%indxbnds_count
5282 nn1 = obsrv%indxbnds(j)
5283 if (nn1 < 1 .or. nn1 > this%nlakes)
then
5284 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
5285 trim(adjustl(obsrv%ObsTypeId)), &
5286 ' lake must be > 0 and <=', this%nlakes, &
5287 '(specified value is ', nn1,
')'
5299 end subroutine lak_rp_obs
5308 subroutine lak_process_obsid(obsrv, dis, inunitobs, iout)
5312 integer(I4B),
intent(in) :: inunitobs
5313 integer(I4B),
intent(in) :: iout
5315 integer(I4B) :: nn1, nn2
5316 integer(I4B) :: icol, istart, istop
5317 character(len=LINELENGTH) :: string
5318 character(len=LENBOUNDNAME) :: bndname
5320 string = obsrv%IDstring
5328 obsrv%FeatureName = bndname
5330 if (obsrv%ObsTypeId ==
'LAK' .or. obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5331 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5333 if (len_trim(bndname) < 1 .and. nn2 < 0)
then
5334 write (
errmsg,
'(a,1x,a,a,1x,a,1x,a)') &
5335 'For observation type', trim(adjustl(obsrv%ObsTypeId)), &
5336 ', ID given as an integer and not as boundname,', &
5337 'but ID2 (iconn) is missing. Either change ID to valid', &
5338 'boundname or supply valid entry for ID2.'
5342 obsrv%FeatureName = bndname
5346 obsrv%NodeNumber2 = nn2
5351 obsrv%NodeNumber = nn1
5352 end subroutine lak_process_obsid
5360 subroutine lak_accumulate_chterm(this, ilak, rrate, chratin, chratout)
5362 class(laktype) :: this
5363 integer(I4B),
intent(in) :: ilak
5364 real(DP),
intent(in) :: rrate
5365 real(DP),
intent(inout) :: chratin
5366 real(DP),
intent(inout) :: chratout
5371 if (this%iboundpak(ilak) < 0)
then
5373 this%chterm(ilak) = this%chterm(ilak) + q
5379 chratout = chratout - q
5383 chratin = chratin + q
5386 end subroutine lak_accumulate_chterm
5390 subroutine lak_bound_update(this)
5392 class(laktype),
intent(inout) :: this
5394 integer(I4B) :: j, n, node
5395 real(DP) :: hlak, head, clak
5398 if (this%nbound == 0)
return
5401 do n = 1, this%nlakes
5402 hlak = this%xnewpak(n)
5403 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5404 node = this%cellid(j)
5405 head = this%xnew(node)
5406 call this%lak_calculate_conn_conductance(n, j, hlak, head, clak)
5407 this%bound(1, j) = hlak
5408 this%bound(2, j) = clak
5411 end subroutine lak_bound_update
5421 subroutine lak_solve(this, update, only_legacy)
5425 class(laktype),
intent(inout) :: this
5426 logical(LGP),
intent(in),
optional :: update
5427 logical(LGP),
intent(in),
optional :: only_legacy
5429 logical(LGP) :: lupdate
5430 logical(LGP) :: legonly
5433 integer(I4B) :: iicnvg
5434 integer(I4B) :: iter
5435 integer(I4B) :: maxiter
5436 integer(I4B) :: ncnv
5448 if (
present(update))
then
5455 if (
present(only_legacy))
then
5456 legonly = only_legacy
5465 do n = 1, this%nlakes
5469 if (legonly .and. this%ilegacy(n) == 0)
then
5474 this%surfin(n) =
dzero
5475 this%surfout(n) =
dzero
5476 this%surfout1(n) =
dzero
5477 if (this%xnewpak(n) < this%lakebot(n))
then
5478 this%xnewpak(n) = this%lakebot(n)
5480 if (this%gwfiss /= 0)
then
5481 this%xoldpak(n) = this%xnewpak(n)
5486 this%en1(n) = this%lakebot(n)
5487 call this%lak_calculate_residual(n, this%en1(n), this%r1(n))
5488 this%en2(n) = this%laketop(n)
5489 call this%lak_calculate_residual(n, this%en2(n), this%r2(n))
5493 do n = 1, this%noutlets
5494 if (legonly .and. this%ilegacy(this%lakein(n)) == 0) cycle
5495 this%simoutrate(n) =
dzero
5499 do n = 1, this%nlakes
5500 call this%lak_calculate_outlet_inflow(n, this%surfin(n))
5505 do n = 1, this%nlakes
5506 hlak0 = this%xoldpak(n)
5507 hlak = this%xnewpak(n)
5508 call this%lak_calculate_runoff(n, ro)
5509 call this%lak_calculate_inflow(n, qinf)
5510 call this%lak_calculate_external(n, ex)
5511 call this%lak_calculate_vol(n, hlak0, v0)
5512 call this%lak_calculate_vol(n, hlak, v1)
5513 this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
5518 do n = 1, this%nlakes
5519 call this%lak_calculate_outlet_inflow(n, outinf)
5520 this%flwin(n) = this%flwin(n) + outinf
5524 maxiter = this%maxlakit
5527 converge:
do iter = 1, maxiter
5529 do n = 1, this%nlakes
5530 if (this%ncncvr(n) == 0) ncnv = 1
5532 if (iter == maxiter) ncnv = 0
5533 if (ncnv == 0) iicnvg = 1
5536 do n = 1, this%nlakes
5537 this%evap(n) =
dzero
5538 this%precip(n) =
dzero
5539 this%precip1(n) =
dzero
5540 this%seep(n) =
dzero
5541 this%seep1(n) =
dzero
5542 this%evap(n) =
dzero
5543 this%evap1(n) =
dzero
5544 this%evapo(n) =
dzero
5545 this%withr(n) =
dzero
5546 this%withr1(n) =
dzero
5547 this%flwiter(n) = this%flwin(n)
5548 this%flwiter1(n) = this%flwin(n)
5549 if (this%gwfiss /= 0)
then
5550 this%flwiter(n) =
dep20
5551 this%flwiter1(n) =
dep20
5553 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5554 this%hcof(j) =
dzero
5559 do n = 1, this%nlakes
5560 if (legonly .and. this%ilegacy(n) == 0) cycle
5561 call this%lak_estimate_seepage_single(n, ncnv)
5564 laklevel:
do n = 1, this%nlakes
5565 if (legonly .and. this%ilegacy(n) == 0) cycle laklevel
5566 call this%lak_solve_single(n, iter, maxiter, ncnv, lupdate)
5569 if (iicnvg == 1)
exit converge
5578 if (this%imover == 1 .and. .not. legonly)
then
5579 do n = 1, this%noutlets
5580 call this%pakmvrobj%accumulate_qformvr(n, -this%simoutrate(n))
5583 end subroutine lak_solve
5593 subroutine lak_solve_single(this, n, iter, maxiter, ncnv, lupdate)
5597 class(laktype),
intent(inout) :: this
5598 integer(I4B),
intent(in) :: n
5599 integer(I4B),
intent(in) :: iter
5600 integer(I4B),
intent(in) :: maxiter
5601 integer(I4B),
intent(in) :: ncnv
5602 logical(LGP),
intent(in) :: lupdate
5604 integer(I4B) :: ibflg
5605 integer(I4B) :: idhp
5625 real(DP) :: qtolfact
5631 if (this%iboundpak(n) == 0)
then
5636 hlak = this%xnewpak(n)
5637 if (iter < maxiter)
then
5638 this%stageiter(n) = this%xnewpak(n)
5640 call this%lak_calculate_rainfall(n, hlak, ra)
5642 this%flwiter(n) = this%flwiter(n) + ra
5643 call this%lak_calculate_rainfall(n, hlak + delh, ra)
5644 this%precip1(n) = ra
5645 this%flwiter1(n) = this%flwiter1(n) + ra
5648 call this%lak_calculate_withdrawal(n, this%flwiter(n), wr)
5650 call this%lak_calculate_withdrawal(n, this%flwiter1(n), wr)
5654 call this%lak_calculate_evaporation(n, hlak, this%flwiter(n), ev)
5656 call this%lak_calculate_evaporation(n, hlak + delh, this%flwiter1(n), ev)
5660 call this%lak_calculate_outlet_outflow(n, hlak + delh, &
5663 call this%lak_calculate_outlet_outflow(n, hlak, this%flwiter(n), &
5667 call this%lak_calculate_outlet_inflow(n, this%surfin(n))
5671 if (this%iboundpak(n) > 0 .and. lupdate .eqv. .true.)
then
5674 hlak0 = this%xoldpak(n)
5675 hlak = this%xnewpak(n)
5676 call this%lak_calculate_vol(n, hlak0, v0)
5677 call this%lak_calculate_vol(n, hlak, v1)
5678 call this%lak_calculate_runoff(n, ro)
5679 call this%lak_calculate_inflow(n, qinf)
5680 call this%lak_calculate_external(n, ex)
5681 this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
5685 resid = this%precip(n) + this%evap(n) + this%withr(n) + ro + &
5686 qinf + ex + this%surfin(n) + &
5687 this%surfout(n) + this%seep(n)
5688 resid1 = this%precip1(n) + this%evap1(n) + this%withr1(n) + ro + &
5689 qinf + ex + this%surfin(n) + &
5690 this%surfout1(n) + this%seep1(n)
5693 hlak = this%xnewpak(n)
5694 if (this%gwfiss /= 1)
then
5695 call this%lak_calculate_vol(n, hlak, v1)
5696 resid = resid + (v0 - v1) /
delt
5697 call this%lak_calculate_vol(n, hlak + delh, v1)
5698 resid1 = resid1 + (v0 - v1) /
delt
5702 if (abs(resid1 - resid) >
dzero)
then
5703 derv = (resid1 - resid) / delh
5705 if (abs(derv) >
dprec)
then
5709 if (resid <
dzero)
then
5712 call this%lak_vol2stage(n, resid, dh)
5719 if (iter == 1) this%dh0(n) = dh
5721 adh0 = abs(this%dh0(n))
5722 if ((ts >= this%en2(n)) .or. (ts < this%en1(n)))
then
5725 if ((adh > adh0) .or. (ts - this%lakebot(n)) <
dprec)
then
5727 call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
5733 this%seep0(n) = this%seep(n)
5737 if (this%seep(n) * this%seep0(n) <
dprec)
then
5738 this%iseepc(n) = this%iseepc(n) + 1
5744 if (dh * this%dh0(n) <
dprec) idhp = 1
5747 if (adh > adh0) idhp = 1
5750 this%idhc(n) = this%idhc(n) + 1
5755 if (ibflg == 1)
then
5756 if (this%iseepc(n) > 7 .or. this%idhc(n) > 12)
then
5757 call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
5766 if (hlak < this%lakebot(n))
then
5767 hlak = this%lakebot(n)
5771 call this%lak_calculate_sarea(n, hlak, area)
5774 if (area >
dzero)
then
5775 qtolfact =
delt / area
5781 call this%lak_calculate_residual(n, hlak, resid)
5785 if (abs(dh) < delh .and. abs(resid) * qtolfact < this%dmaxchg)
then
5788 this%xnewpak(n) = hlak
5791 this%seep0(n) = this%seep(n)
5794 end subroutine lak_solve_single
5807 subroutine lak_estimate_seepage_single(this, n, ncnv)
5809 class(laktype),
intent(inout) :: this
5810 integer(I4B),
intent(in) :: n
5811 integer(I4B),
intent(in) :: ncnv
5815 integer(I4B) :: igwfnode
5816 integer(I4B) :: idry
5817 integer(I4B) :: idry1
5829 if (this%iboundpak(n) == 0)
return
5831 estseep:
do i = 1, 2
5833 if (this%gwfiss /= 0)
then
5834 this%xoldpak(n) = this%xnewpak(n)
5836 hlak = this%xnewpak(n)
5837 calcconnseep:
do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5838 igwfnode = this%cellid(j)
5839 head = this%xnew(igwfnode)
5840 if (this%ncncvr(n) /= 2)
then
5841 if (this%ibound(igwfnode) > 0)
then
5842 call this%lak_estimate_conn_exchange(i, n, j, idry, hlak, &
5846 call this%lak_estimate_conn_exchange(i, n, j, idry1, &
5847 hlak + delh, head, qlakgw1, &
5851 if (ncnv == 0 .and. i == 2)
then
5852 if (j == this%maxbound)
then
5856 this%hcof(j) = gwfhcof
5857 this%rhs(j) = gwfrhs
5859 this%hcof(j) =
dzero
5860 this%rhs(j) = qlakgw
5864 this%seep(n) = this%seep(n) + qlakgw
5865 this%seep1(n) = this%seep1(n) + qlakgw1
5871 end subroutine lak_estimate_seepage_single
5877 subroutine lak_bisection(this, n, ibflg, hlak, temporary_stage, dh, residual)
5879 class(laktype),
intent(inout) :: this
5880 integer(I4B),
intent(in) :: n
5881 integer(I4B),
intent(inout) :: ibflg
5882 real(DP),
intent(in) :: hlak
5883 real(DP),
intent(inout) :: temporary_stage
5884 real(DP),
intent(inout) :: dh
5885 real(DP),
intent(inout) :: residual
5888 real(DP) :: temporary_stage0
5889 real(DP) :: residuala
5890 real(DP) :: endpoint1
5891 real(DP) :: endpoint2
5894 temporary_stage0 = hlak
5895 endpoint1 = this%en1(n)
5896 endpoint2 = this%en2(n)
5897 call this%lak_calculate_residual(n, temporary_stage, residuala)
5898 if (hlak > endpoint1 .and. hlak < endpoint2)
then
5901 do i = 1, this%maxlakit
5902 temporary_stage =
dhalf * (endpoint1 + endpoint2)
5903 call this%lak_calculate_residual(n, temporary_stage, residual)
5904 if (abs(residual) ==
dzero .or. &
5905 abs(temporary_stage0 - temporary_stage) < this%dmaxchg)
then
5908 call this%lak_calculate_residual(n, endpoint1, residuala)
5911 if (sign(
done, residuala) == sign(
done, residual))
then
5912 endpoint1 = temporary_stage
5915 endpoint2 = temporary_stage
5917 temporary_stage0 = temporary_stage
5919 dh = hlak - temporary_stage
5920 end subroutine lak_bisection
5925 subroutine lak_calculate_available(this, n, hlak, avail, &
5926 ra, ro, qinf, ex, headp)
5930 class(laktype),
intent(inout) :: this
5931 integer(I4B),
intent(in) :: n
5932 real(DP),
intent(in) :: hlak
5933 real(DP),
intent(inout) :: avail
5934 real(DP),
intent(inout) :: ra
5935 real(DP),
intent(inout) :: ro
5936 real(DP),
intent(inout) :: qinf
5937 real(DP),
intent(inout) :: ex
5938 real(DP),
intent(in),
optional :: headp
5941 integer(I4B) :: idry
5942 integer(I4B) :: igwfnode
5949 if (
present(headp))
then
5959 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5960 igwfnode = this%cellid(j)
5961 if (this%ibound(igwfnode) == 0) cycle
5962 head = this%xnew(igwfnode) + hp
5963 call this%lak_estimate_conn_exchange(1, n, j, idry, hlak, head, qlakgw, &
5968 call this%lak_calculate_rainfall(n, hlak, ra)
5972 call this%lak_calculate_runoff(n, ro)
5976 call this%lak_calculate_inflow(n, qinf)
5977 avail = avail + qinf
5980 call this%lak_calculate_external(n, ex)
5984 call this%lak_calculate_vol(n, this%xoldpak(n), v0)
5985 avail = avail + v0 /
delt
5986 end subroutine lak_calculate_available
5990 subroutine lak_calculate_residual(this, n, hlak, resid, headp)
5994 class(laktype),
intent(inout) :: this
5995 integer(I4B),
intent(in) :: n
5996 real(DP),
intent(in) :: hlak
5997 real(DP),
intent(inout) :: resid
5998 real(DP),
intent(in),
optional :: headp
6001 integer(I4B) :: idry
6002 integer(I4B) :: igwfnode
6021 if (
present(headp))
then
6033 call this%lak_calculate_available(n, hlak, avail, &
6034 ra, ro, qinf, ex, hp)
6037 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
6038 igwfnode = this%cellid(j)
6039 if (this%ibound(igwfnode) == 0) cycle
6040 head = this%xnew(igwfnode) + hp
6041 call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, head, qlakgw, &
6043 seep = seep + qlakgw
6047 call this%lak_calculate_withdrawal(n, avail, wr)
6050 call this%lak_calculate_evaporation(n, hlak, avail, ev)
6053 call this%lak_calculate_outlet_outflow(n, hlak, avail, sout)
6056 call this%lak_calculate_outlet_inflow(n, sin)
6059 resid = ra + ev + wr + ro + qinf + ex + sin + sout + seep
6062 if (this%gwfiss /= 1)
then
6063 hlak0 = this%xoldpak(n)
6064 call this%lak_calculate_vol(n, hlak0, v0)
6065 call this%lak_calculate_vol(n, hlak, v1)
6066 resid = resid + (v0 - v1) /
delt
6068 end subroutine lak_calculate_residual
6072 subroutine lak_setup_budobj(this)
6076 class(laktype) :: this
6078 integer(I4B) :: nbudterm
6079 integer(I4B) :: nlen
6080 integer(I4B) :: j, n, n1, n2
6081 integer(I4B) :: maxlist, naux
6084 character(len=LENBUDTXT) :: text
6085 character(len=LENBUDTXT),
dimension(1) :: auxtxt
6091 do n = 1, this%noutlets
6092 if (this%lakein(n) > 0 .and. this%lakeout(n) > 0)
then
6096 if (nlen > 0) nbudterm = nbudterm + 1
6097 if (this%imover == 1) nbudterm = nbudterm + 2
6098 if (this%naux > 0) nbudterm = nbudterm + 1
6102 call this%budobj%budgetobject_df(this%nlakes, nbudterm, 0, 0, &
6103 ibudcsv=this%ibudcsv)
6109 text =
' FLOW-JA-FACE'
6113 call this%budobj%budterm(idx)%initialize(text, &
6118 maxlist, .false., .false., &
6119 naux, ordered_id1=.false.)
6122 call this%budobj%budterm(idx)%reset(2 * nlen)
6124 do n = 1, this%noutlets
6126 n2 = this%lakeout(n)
6127 if (n1 > 0 .and. n2 > 0)
then
6128 call this%budobj%budterm(idx)%update_term(n1, n2, q)
6129 call this%budobj%budterm(idx)%update_term(n2, n1, -q)
6137 maxlist = this%maxbound
6139 auxtxt(1) =
' FLOW-AREA'
6140 call this%budobj%budterm(idx)%initialize(text, &
6145 maxlist, .false., .true., &
6147 call this%budobj%budterm(idx)%reset(this%maxbound)
6149 do n = 1, this%nlakes
6150 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
6152 call this%budobj%budterm(idx)%update_term(n, n2, q)
6159 maxlist = this%nlakes
6161 call this%budobj%budterm(idx)%initialize(text, &
6166 maxlist, .false., .false., &
6170 text =
' EVAPORATION'
6172 maxlist = this%nlakes
6174 call this%budobj%budterm(idx)%initialize(text, &
6179 maxlist, .false., .false., &
6185 maxlist = this%nlakes
6187 call this%budobj%budterm(idx)%initialize(text, &
6192 maxlist, .false., .false., &
6196 text =
' EXT-INFLOW'
6198 maxlist = this%nlakes
6200 call this%budobj%budterm(idx)%initialize(text, &
6205 maxlist, .false., .false., &
6209 text =
' WITHDRAWAL'
6211 maxlist = this%nlakes
6213 call this%budobj%budterm(idx)%initialize(text, &
6218 maxlist, .false., .false., &
6222 text =
' EXT-OUTFLOW'
6224 maxlist = this%nlakes
6226 call this%budobj%budterm(idx)%initialize(text, &
6231 maxlist, .false., .false., &
6237 maxlist = this%nlakes
6239 auxtxt(1) =
' VOLUME'
6240 call this%budobj%budterm(idx)%initialize(text, &
6245 maxlist, .false., .false., &
6251 maxlist = this%nlakes
6253 call this%budobj%budterm(idx)%initialize(text, &
6258 maxlist, .false., .false., &
6262 if (this%imover == 1)
then
6267 maxlist = this%nlakes
6269 call this%budobj%budterm(idx)%initialize(text, &
6274 maxlist, .false., .false., &
6280 maxlist = this%noutlets
6282 call this%budobj%budterm(idx)%initialize(text, &
6287 maxlist, .false., .false., &
6288 naux, ordered_id1=.false.)
6291 call this%budobj%budterm(idx)%reset(this%noutlets)
6293 do n = 1, this%noutlets
6295 call this%budobj%budterm(idx)%update_term(n1, n1, q)
6306 maxlist = this%nlakes
6307 call this%budobj%budterm(idx)%initialize(text, &
6312 maxlist, .false., .false., &
6317 if (this%iprflow /= 0)
then
6318 call this%budobj%flowtable_df(this%iout)
6320 end subroutine lak_setup_budobj
6324 subroutine lak_fill_budobj(this)
6326 class(laktype) :: this
6328 integer(I4B) :: naux
6329 real(DP),
dimension(:),
allocatable :: auxvartmp
6338 integer(I4B) :: nlen
6341 real(DP) :: lkstg, gwhead, wa
6348 do n = 1, this%noutlets
6349 if (this%lakein(n) > 0 .and. this%lakeout(n) > 0)
then
6355 call this%budobj%budterm(idx)%reset(2 * nlen)
6356 do n = 1, this%noutlets
6358 n2 = this%lakeout(n)
6359 if (n1 > 0 .and. n2 > 0)
then
6360 q = this%simoutrate(n)
6361 if (this%imover == 1)
then
6362 q = q + this%pakmvrobj%get_qtomvr(n)
6364 call this%budobj%budterm(idx)%update_term(n1, n2, q)
6365 call this%budobj%budterm(idx)%update_term(n2, n1, -q)
6372 call this%budobj%budterm(idx)%reset(this%maxbound)
6373 do n = 1, this%nlakes
6374 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
6377 lkstg = this%xnewpak(n)
6381 gwhead = this%xnew(n2)
6382 call this%lak_calculate_conn_warea(n, j, lkstg, gwhead, wa)
6386 if (this%belev(j) > lkstg) wa =
dzero
6387 this%qauxcbc(1) = wa
6388 call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
6394 call this%budobj%budterm(idx)%reset(this%nlakes)
6395 do n = 1, this%nlakes
6397 call this%budobj%budterm(idx)%update_term(n, n, q)
6402 call this%budobj%budterm(idx)%reset(this%nlakes)
6403 do n = 1, this%nlakes
6405 call this%budobj%budterm(idx)%update_term(n, n, q)
6410 call this%budobj%budterm(idx)%reset(this%nlakes)
6411 do n = 1, this%nlakes
6413 call this%budobj%budterm(idx)%update_term(n, n, q)
6418 call this%budobj%budterm(idx)%reset(this%nlakes)
6419 do n = 1, this%nlakes
6421 call this%budobj%budterm(idx)%update_term(n, n, q)
6426 call this%budobj%budterm(idx)%reset(this%nlakes)
6427 do n = 1, this%nlakes
6429 call this%budobj%budterm(idx)%update_term(n, n, q)
6434 call this%budobj%budterm(idx)%reset(this%nlakes)
6435 do n = 1, this%nlakes
6436 call this%lak_get_external_outlet(n, q)
6438 call this%lak_get_external_mover(n, v)
6440 call this%budobj%budterm(idx)%update_term(n, n, q)
6445 call this%budobj%budterm(idx)%reset(this%nlakes)
6446 do n = 1, this%nlakes
6447 call this%lak_calculate_vol(n, this%xnewpak(n), v1)
6449 this%qauxcbc(1) = v1
6450 call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc)
6455 call this%budobj%budterm(idx)%reset(this%nlakes)
6456 do n = 1, this%nlakes
6458 call this%budobj%budterm(idx)%update_term(n, n, q)
6462 if (this%imover == 1)
then
6466 call this%budobj%budterm(idx)%reset(this%nlakes)
6467 do n = 1, this%nlakes
6468 q = this%pakmvrobj%get_qfrommvr(n)
6469 call this%budobj%budterm(idx)%update_term(n, n, q)
6474 call this%budobj%budterm(idx)%reset(this%noutlets)
6475 do n = 1, this%noutlets
6477 q = this%pakmvrobj%get_qtomvr(n)
6481 call this%budobj%budterm(idx)%update_term(n1, n1, q)
6490 allocate (auxvartmp(naux))
6491 call this%budobj%budterm(idx)%reset(this%nlakes)
6492 do n = 1, this%nlakes
6496 auxvartmp(jj) = this%lauxvar(jj, ii)
6498 call this%budobj%budterm(idx)%update_term(n, n, q, auxvartmp)
6500 deallocate (auxvartmp)
6504 call this%budobj%accumulate_terms()
6505 end subroutine lak_fill_budobj
6512 subroutine lak_setup_tableobj(this)
6516 class(laktype) :: this
6518 integer(I4B) :: nterms
6519 character(len=LINELENGTH) :: title
6520 character(len=LINELENGTH) :: text
6523 if (this%iprhed > 0)
then
6530 if (this%inamedbound == 1)
then
6535 title = trim(adjustl(this%text))//
' PACKAGE ('// &
6536 trim(adjustl(this%packName))//
') STAGES FOR EACH CONTROL VOLUME'
6539 call table_cr(this%stagetab, this%packName, title)
6540 call this%stagetab%table_df(this%nlakes, nterms, this%iout, &
6544 if (this%inamedbound == 1)
then
6546 call this%stagetab%initialize_column(text, 20, alignment=tableft)
6551 call this%stagetab%initialize_column(text, 10, alignment=tabcenter)
6555 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6558 text =
'SURFACE AREA'
6559 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6562 text =
'WETTED AREA'
6563 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6567 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6569 end subroutine lak_setup_tableobj
6573 subroutine lak_activate_density(this)
6575 class(laktype),
intent(inout) :: this
6577 integer(I4B) :: i, j
6580 if (this%iimplicit /= 0)
then
6582 'The IMPLICIT option cannot be used with the BUY package until the &
6583 &implicit formulation includes density terms. Remove the IMPLICIT &
6584 &option from LAK package '//trim(this%packName)//
' to simulate &
6587 call this%parser%StoreErrorUnit()
6592 call mem_reallocate(this%denseterms, 3, this%MAXBOUND,
'DENSETERMS', &
6594 do i = 1, this%maxbound
6596 this%denseterms(j, i) =
dzero
6599 write (this%iout,
'(/1x,a)')
'DENSITY TERMS HAVE BEEN ACTIVATED FOR LAKE &
6600 &PACKAGE: '//trim(adjustl(this%packName))
6601 end subroutine lak_activate_density
6607 subroutine lak_activate_viscosity(this)
6611 class(laktype),
intent(inout) :: this
6618 call mem_reallocate(this%viscratios, 2, this%MAXBOUND,
'VISCRATIOS', &
6620 do i = 1, this%maxbound
6622 this%viscratios(j, i) =
done
6625 write (this%iout,
'(/1x,a)')
'VISCOSITY HAS BEEN ACTIVATED FOR LAK &
6626 &PACKAGE: '//trim(adjustl(this%packName))
6627 end subroutine lak_activate_viscosity
6647 subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, &
6648 botl, flow, gwfhcof, gwfrhs)
6650 class(laktype),
intent(inout) :: this
6651 integer(I4B),
intent(in) :: iconn
6652 real(DP),
intent(in) :: stage
6653 real(DP),
intent(in) :: head
6654 real(DP),
intent(in) :: cond
6655 real(DP),
intent(in) :: botl
6656 real(DP),
intent(inout) :: flow
6657 real(DP),
intent(inout) :: gwfhcof
6658 real(DP),
intent(inout) :: gwfrhs
6663 real(DP) :: rdenselak
6664 real(DP) :: rdensegwf
6665 real(DP) :: rdenseavg
6671 logical(LGP) :: stage_below_bot
6672 logical(LGP) :: head_below_bot
6675 if (stage >= botl)
then
6677 stage_below_bot = .false.
6678 rdenselak = this%denseterms(1, iconn)
6681 stage_below_bot = .true.
6682 rdenselak = this%denseterms(2, iconn)
6686 if (head >= botl)
then
6688 head_below_bot = .false.
6689 rdensegwf = this%denseterms(2, iconn)
6692 head_below_bot = .true.
6693 rdensegwf = this%denseterms(1, iconn)
6697 if (rdensegwf ==
dzero)
return
6700 if (stage_below_bot .and. head_below_bot)
then
6707 rdenseavg =
dhalf * (rdenselak + rdensegwf)
6711 d1 = cond * (rdenseavg -
done)
6712 gwfhcof = gwfhcof - d1
6713 gwfrhs = gwfrhs - d1 * ss
6718 if (.not. stage_below_bot .and. .not. head_below_bot)
then
6722 elevgwf = this%denseterms(3, iconn)
6723 if (this%ictype(iconn) == 0 .or. this%ictype(iconn) == 3)
then
6730 elevavg =
dhalf * (elevlak + elevgwf)
6731 havg =
dhalf * (hh + ss)
6732 d2 = cond * (havg - elevavg) * (rdensegwf - rdenselak)
6733 gwfrhs = gwfrhs + d2
6737 end subroutine lak_calculate_density_exchange
This module contains block parser methods.
This module contains the base boundary package.
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
@ tabucstring
upper case string table data
@ tabstring
string table data
@ tabinteger
integer table data
integer(i4b), parameter lenpackagename
maximum length of the package name
real(dp), parameter dp9
real constant 9/10
integer(i4b), parameter iwetlake
integer constant for a dry lake
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 dhnoflo
real no flow constant
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 dem1
real constant 1e-1
integer(i4b), parameter maxadpit
maximum advanced package Newton-Raphson iterations
integer(i4b), parameter lenvarname
maximum length of a variable name
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 dgravity
real constant gravitational acceleration (m/(s s))
real(dp), parameter dpi
real constant
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 dcd
real constant weir coefficient in SI units
real(dp), parameter dem5
real constant 1e-5
real(dp), parameter dten
real constant 10
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 dem9
real constant 1e-9
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 dthree
real constant 3
real(dp), parameter done
real constant 1
integer(i4b) function, public get_node(ilay, irow, icol, nlay, nrow, ncol)
Get node number, given layer, row, and column indices for a structured grid. If any argument is inval...
This module defines variable data types.
subroutine lak_activate_density(this)
Activate addition of density terms.
subroutine lak_cq(this, x, flowja, iadv)
Calculate flows.
subroutine lak_read_outlets(this)
Read the lake outlets for this package.
subroutine lak_nur(this, neqpak, x, xtemp, dx, inewtonur, dxmax, locmax)
Apply Newton under-relaxation to the lake stage.
subroutine lak_vol2stage(this, ilak, vol, stage)
Determine the stage from a provided volume.
subroutine lak_calculate_outlet_outflow(this, ilak, stage, avail, outoutf)
Calculate the outlet outflow from a lake.
subroutine lak_read_lake_connections(this)
Read the lake connections for this package.
subroutine lak_calculate_sarea(this, ilak, stage, sarea)
Calculate the surface area of a lake at a given stage.
subroutine lak_fn(this, rhs, ia, idxglo, matrix_sln)
Fill newton terms.
subroutine lak_read_tables(this)
Read the lake tables for this package.
subroutine lak_set_stressperiod(this, itemno)
Set a stress period attribute for lakweslls(itemno) using keywords.
subroutine lak_ac(this, moffset, sparse)
Add the lake rows and columns to the sparse matrix.
subroutine lak_ot_package_flows(this, icbcfl, ibudfl)
Output LAK package flow terms.
subroutine lak_accumulate_chterm(this, ilak, rrate, chratin, chratout)
Accumulate constant head terms for budget.
subroutine lak_get_external_mover(this, ilak, outoutf)
Get the mover outflow from a lake to an external boundary.
subroutine lak_cf(this)
Formulate the HCOF and RHS terms.
subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, botl, flow, gwfhcof, gwfrhs)
Calculate the groundwater-lake density exchange terms.
subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
Calculate the evaporation from a lake at a provided stage subject to an available volume.
subroutine lak_read_lakes(this)
Read the dimensions for this package.
subroutine lak_setup_budobj(this)
Set up the budget object that stores all the lake flows.
subroutine lak_calculate_external(this, ilak, ex)
Calculate the external flow terms to a lake.
subroutine lak_calculate_warea(this, ilak, stage, warea, hin)
Calculate the wetted area of a lake at a given stage.
subroutine lak_calculate_vol(this, ilak, stage, volume)
Calculate the volume of a lake at a given stage.
subroutine lak_da(this)
Deallocate objects.
subroutine lak_calculate_storagechange(this, ilak, stage, stage0, delt, dvr)
Calculate the storage change in a lake based on provided stages and a passed delt.
subroutine lak_get_internal_outlet(this, ilak, outoutf)
Get the outlet from a lake to another lake.
subroutine lak_calculate_conductance(this, ilak, stage, conductance)
Calculate the total conductance for a lake at a provided stage.
subroutine laktables_to_vectors(this, laketables)
Copy the laketables structure data into flattened vectors that are stored in the memory manager.
subroutine lak_calculate_residual(this, n, hlak, resid, headp)
Calculate the residual for a lake given a passed stage.
subroutine lak_get_outlet_tomover(this, ilak, outoutf)
Get the outlet to mover from a lake.
subroutine lak_allocate_arrays(this)
Allocate scalar members.
subroutine lak_calculate_available(this, n, hlak, avail, ra, ro, qinf, ex, headp)
Calculate the available volumetric rate for a lake given a passed stage.
logical function lak_obs_supported(this)
Procedures related to observations (type-bound)
subroutine lak_set_pointers(this, neq, ibound, xnew, xold, flowja)
Set pointers to model arrays and variables so that a package has access to these things.
subroutine lak_calculate_inflow(this, ilak, qin)
Calculate specified inflow to a lake.
subroutine lak_calculate_conn_exchange_deriv(this, ilak, iconn, stage, head, flow, dqds, dqdh)
Lakebed seepage and its derivatives for a connection (IMPLICIT)
character(len=lenpackagename) text
subroutine lak_mc(this, moffset, matrix_sln)
Find the matrix position of each lake row and connection.
subroutine lak_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
Write flows to binary file and/or print flows to budget.
subroutine lak_rp_obs(this)
Process each observation.
subroutine lak_calculate_conn_conductance(this, ilak, iconn, stage, head, cond)
Calculate the conductance for a lake connection at a provided stage and groundwater head.
subroutine lak_calculate_exchange(this, ilak, stage, totflow)
Calculate the total groundwater-lake flow at a provided stage.
subroutine lak_linear_interpolation(this, n, x, y, z, v)
Perform linear interpolation of two vectors.
subroutine lak_setup_tableobj(this)
Set up the table object that is used to write the lak stage data.
subroutine lak_ot_dv(this, idvsave, idvprint)
Save LAK-calculated values to binary file.
subroutine lak_options(this, option, found)
Set options specific to LakType.
subroutine lak_get_external_outlet(this, ilak, outoutf)
Get the outlet outflow from a lake to an external boundary.
subroutine lak_calculate_rainfall(this, ilak, stage, ra)
Calculate the rainfall for a lake.
subroutine lak_get_internal_inlet(this, ilak, outinf)
Get the outlet inflow to a lake from another lake.
subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
Final convergence check for package.
subroutine lak_rp(this)
Read and Prepare.
subroutine lak_read_dimensions(this)
Read the dimensions for this package.
subroutine lak_activate_viscosity(this)
Activate viscosity terms.
subroutine lak_read_initial_attr(this)
Read the initial parameters for this package.
integer(i4b) function lak_check_valid(this, itemno)
Determine if a valid lake or outlet number has been specified.
subroutine lak_bisection(this, n, ibflg, hlak, temporary_stage, dh, residual)
@ brief Lake package bisection method
subroutine lak_ad(this)
Add package connection to matrix.
character(len=lenftype) ftype
subroutine lak_fc(this, rhs, ia, idxglo, matrix_sln)
Copy rhs and hcof into solution rhs and amat.
subroutine lak_set_attribute_error(this, ilak, keyword, msg)
Issue a parameter error for lakweslls(ilak)
subroutine lak_calculate_runoff(this, ilak, ro)
Calculate runoff to a lake.
subroutine define_listlabel(this)
Define the list heading that is written to iout when PRINT_INPUT option is used.
subroutine lak_read_table(this, ilak, filename, laketable)
Read the lake table for this package.
subroutine lak_calculate_withdrawal(this, ilak, avail, wr)
Calculate the withdrawal from a lake subject to an available volume.
subroutine lak_calculate_conn_warea(this, ilak, iconn, stage, head, wa)
Calculate the wetted area of a lake connection at a given stage.
subroutine lak_solve(this, update, only_legacy)
Solve for lake stage.
subroutine lak_estimate_seepage_single(this, n, ncnv)
Estimate the lakebed seepage for a single lake.
subroutine lak_df_obs(this)
Store observation type supported by LAK package. Overrides BndTypebnd_df_obs.
subroutine lak_allocate_scalars(this)
Allocate scalar members.
subroutine lak_bound_update(this)
Store the lake head and connection conductance in the bound array.
subroutine lak_calculate_cond_head(this, iconn, stage, head, vv)
Calculate the controlling lake stage or groundwater head used to calculate the conductance for a lake...
subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl)
Write LAK budget to listing file.
subroutine, public lak_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
Create a new LAK Package and point bndobj to the new package.
subroutine lak_bd_obs(this)
Calculate observations this time step and call ObsTypeSaveOneSimval for each LakType observation.
subroutine lak_estimate_conn_exchange(this, iflag, ilak, iconn, idry, stage, head, flow, source, gwfhcof, gwfrhs)
Calculate the groundwater-lake flow at a provided stage and groundwater head.
subroutine lak_fill_budobj(this)
Copy flow terms into thisbudobj.
subroutine lak_get_internal_mover(this, ilak, outoutf)
Get the mover outflow from a lake to another lake.
subroutine lak_calculate_outlet_inflow(this, ilak, outinf)
Calculate the outlet inflow to a lake.
subroutine lak_outlet_outflow_rate(this, ilak, stage, qout)
Total uncapped outlet outflow rate from a lake at a provided stage.
subroutine lak_calculate_conn_exchange(this, ilak, iconn, stage, head, flow, gwfhcof, gwfrhs)
Calculate the groundwater-lake flow at a provided stage and groundwater head.
subroutine lak_ar(this)
Allocate and Read.
subroutine lak_solve_single(this, n, iter, maxiter, ncnv, lupdate)
Advance one lake stage by a single substitution iteration.
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
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.
This module contains the derived type ObsType.
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
integer(i4b) ifailedstepretry
current retry for this time step
character(len=maxcharlen) warnmsg
warning message string
real(dp) function squadraticsaturation(top, bot, x, eps)
@ brief sQuadraticSaturation
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
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).