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 :: iforcefb => null()
83 integer(I4B),
pointer :: bditems => null()
86 integer(I4B),
dimension(:),
pointer,
contiguous :: nlakeconn => null()
87 integer(I4B),
dimension(:),
pointer,
contiguous :: idxlakeconn => null()
88 integer(I4B),
dimension(:),
pointer,
contiguous :: ntabrow => null()
89 real(dp),
dimension(:),
pointer,
contiguous :: strt => null()
90 real(dp),
dimension(:),
pointer,
contiguous :: laketop => null()
91 real(dp),
dimension(:),
pointer,
contiguous :: lakebot => null()
92 real(dp),
dimension(:),
pointer,
contiguous :: sareamax => null()
93 character(len=LENBOUNDNAME),
dimension(:),
pointer, &
94 contiguous :: lakename => null()
95 character(len=8),
dimension(:),
pointer,
contiguous :: status => null()
96 real(dp),
dimension(:),
pointer,
contiguous :: avail => null()
97 real(dp),
dimension(:),
pointer,
contiguous :: lkgwsink => null()
98 real(dp),
dimension(:),
pointer,
contiguous :: stage => null()
99 real(dp),
dimension(:),
pointer,
contiguous :: rainfall => null()
100 real(dp),
dimension(:),
pointer,
contiguous :: evaporation => null()
101 real(dp),
dimension(:),
pointer,
contiguous :: runoff => null()
102 real(dp),
dimension(:),
pointer,
contiguous :: inflow => null()
103 real(dp),
dimension(:),
pointer,
contiguous :: withdrawal => null()
104 real(dp),
dimension(:, :),
pointer,
contiguous :: lauxvar => null()
107 integer(I4B),
dimension(:),
pointer,
contiguous :: ialaktab => null()
108 real(dp),
dimension(:),
pointer,
contiguous :: tabstage => null()
109 real(dp),
dimension(:),
pointer,
contiguous :: tabvolume => null()
110 real(dp),
dimension(:),
pointer,
contiguous :: tabsarea => null()
111 real(dp),
dimension(:),
pointer,
contiguous :: tabwarea => null()
114 integer(I4B),
dimension(:),
pointer,
contiguous :: ncncvr => null()
118 integer(I4B),
dimension(:),
pointer,
contiguous :: ifallback => null()
119 integer(I4B),
dimension(:),
pointer,
contiguous :: nstuck => null()
120 real(dp),
dimension(:),
pointer,
contiguous :: surfin => null()
121 real(dp),
dimension(:),
pointer,
contiguous :: surfout => null()
122 real(dp),
dimension(:),
pointer,
contiguous :: surfout1 => null()
123 real(dp),
dimension(:),
pointer,
contiguous :: precip => null()
124 real(dp),
dimension(:),
pointer,
contiguous :: precip1 => null()
125 real(dp),
dimension(:),
pointer,
contiguous :: evap => null()
126 real(dp),
dimension(:),
pointer,
contiguous :: evap1 => null()
127 real(dp),
dimension(:),
pointer,
contiguous :: evapo => null()
128 real(dp),
dimension(:),
pointer,
contiguous :: withr => null()
129 real(dp),
dimension(:),
pointer,
contiguous :: withr1 => null()
130 real(dp),
dimension(:),
pointer,
contiguous :: flwin => null()
131 real(dp),
dimension(:),
pointer,
contiguous :: flwiter => null()
132 real(dp),
dimension(:),
pointer,
contiguous :: flwiter1 => null()
133 real(dp),
dimension(:),
pointer,
contiguous :: seep => null()
134 real(dp),
dimension(:),
pointer,
contiguous :: seep1 => null()
135 real(dp),
dimension(:),
pointer,
contiguous :: seep0 => null()
136 real(dp),
dimension(:),
pointer,
contiguous :: stageiter => null()
137 real(dp),
dimension(:),
pointer,
contiguous :: chterm => null()
140 integer(I4B),
dimension(:),
pointer,
contiguous :: iseepc => null()
141 integer(I4B),
dimension(:),
pointer,
contiguous :: idhc => null()
142 real(dp),
dimension(:),
pointer,
contiguous :: en1 => null()
143 real(dp),
dimension(:),
pointer,
contiguous :: en2 => null()
144 real(dp),
dimension(:),
pointer,
contiguous :: r1 => null()
145 real(dp),
dimension(:),
pointer,
contiguous :: r2 => null()
146 real(dp),
dimension(:),
pointer,
contiguous :: dh0 => null()
147 real(dp),
dimension(:),
pointer,
contiguous :: s0 => null()
148 real(dp),
dimension(:),
pointer,
contiguous :: qgwf0 => null()
150 integer(I4B),
dimension(:),
pointer,
contiguous :: idxlocnode => null()
151 integer(I4B),
dimension(:),
pointer,
contiguous :: idxdiag => null()
152 integer(I4B),
dimension(:),
pointer,
contiguous :: idxoffdglo => null()
153 integer(I4B),
dimension(:),
pointer,
contiguous :: idxsymdglo => null()
154 integer(I4B),
dimension(:),
pointer,
contiguous :: idxsymoffdglo => null()
157 integer(I4B),
dimension(:),
pointer,
contiguous :: imap => null()
158 integer(I4B),
dimension(:),
pointer,
contiguous :: cellid => null()
159 integer(I4B),
dimension(:),
pointer,
contiguous :: nodesontop => null()
160 integer(I4B),
dimension(:),
pointer,
contiguous :: ictype => null()
161 real(dp),
dimension(:),
pointer,
contiguous :: bedleak => null()
162 real(dp),
dimension(:),
pointer,
contiguous :: belev => null()
163 real(dp),
dimension(:),
pointer,
contiguous :: telev => null()
164 real(dp),
dimension(:),
pointer,
contiguous :: connlength => null()
165 real(dp),
dimension(:),
pointer,
contiguous :: connwidth => null()
166 real(dp),
dimension(:),
pointer,
contiguous :: sarea => null()
167 real(dp),
dimension(:),
pointer,
contiguous :: warea => null()
168 real(dp),
dimension(:),
pointer,
contiguous :: satcond => null()
169 real(dp),
dimension(:),
pointer,
contiguous :: simcond => null()
170 real(dp),
dimension(:),
pointer,
contiguous :: simlakgw => null()
173 integer(I4B),
dimension(:),
pointer,
contiguous :: lakein => null()
174 integer(I4B),
dimension(:),
pointer,
contiguous :: lakeout => null()
175 integer(I4B),
dimension(:),
pointer,
contiguous :: iouttype => null()
176 real(dp),
dimension(:),
pointer,
contiguous :: outrate => null()
177 real(dp),
dimension(:),
pointer,
contiguous :: outinvert => null()
178 real(dp),
dimension(:),
pointer,
contiguous :: outwidth => null()
179 real(dp),
dimension(:),
pointer,
contiguous :: outrough => null()
180 real(dp),
dimension(:),
pointer,
contiguous :: outslope => null()
181 real(dp),
dimension(:),
pointer,
contiguous :: simoutrate => null()
184 real(dp),
dimension(:),
pointer,
contiguous :: qauxcbc => null()
185 real(dp),
dimension(:),
pointer,
contiguous :: dbuff => null()
186 real(dp),
dimension(:),
pointer,
contiguous :: qleak => null()
190 real(dp),
dimension(:),
pointer,
contiguous :: holdconn => null()
191 real(dp),
dimension(:),
pointer,
contiguous :: qsto => null()
194 integer(I4B),
pointer :: gwfiss => null()
195 real(dp),
dimension(:),
pointer,
contiguous :: gwfk11 => null()
196 real(dp),
dimension(:),
pointer,
contiguous :: gwfk33 => null()
197 real(dp),
dimension(:),
pointer,
contiguous :: gwfsat => null()
198 integer(I4B),
pointer :: gwfik33 => null()
201 integer(I4B),
dimension(:),
pointer,
contiguous :: iboundpak => null()
202 real(dp),
dimension(:),
pointer,
contiguous :: xnewpak => null()
203 real(dp),
dimension(:),
pointer,
contiguous :: xoldpak => null()
213 integer(I4B),
pointer :: idense
214 real(dp),
dimension(:, :),
pointer,
contiguous :: denseterms => null()
217 real(dp),
dimension(:, :),
pointer,
contiguous :: viscratios => null()
239 procedure,
private :: lak_set_fallback
240 procedure,
private :: lak_check_disconnected
294 procedure,
private :: lak_fc_implicit
295 procedure,
private :: lak_budget_nogwf
316 module subroutine lak_budget_nogwf(this, n, stage, b)
317 class(laktype),
intent(inout) :: this
318 integer(I4B),
intent(in) :: n
319 real(dp),
intent(in) :: stage
320 real(dp),
intent(inout) :: b
325 module subroutine lak_fc_implicit(this, rhs, matrix_sln)
326 class(laktype) :: this
327 real(dp),
dimension(:),
intent(inout) :: rhs
333 module subroutine lak_set_fallback(this, kiter, icnvgmod)
334 class(laktype),
intent(inout) :: this
335 integer(I4B),
intent(in) :: kiter
336 integer(I4B),
intent(in) :: icnvgmod
341 module subroutine lak_check_disconnected(this)
342 class(laktype),
intent(inout) :: this
350 subroutine lak_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
352 class(
bndtype),
pointer :: packobj
353 integer(I4B),
intent(in) :: id
354 integer(I4B),
intent(in) :: ibcnum
355 integer(I4B),
intent(in) :: inunit
356 integer(I4B),
intent(in) :: iout
357 character(len=*),
intent(in) :: namemodel
358 character(len=*),
intent(in) :: pakname
360 type(laktype),
pointer :: lakobj
367 call packobj%set_names(ibcnum, namemodel, pakname, ftype)
371 call lakobj%lak_allocate_scalars()
374 call packobj%pack_initialize()
376 packobj%inunit = inunit
379 packobj%ibcnum = ibcnum
384 end subroutine lak_create
388 subroutine lak_allocate_scalars(this)
390 class(laktype),
intent(inout) :: this
393 call this%BndType%allocate_scalars()
396 call mem_allocate(this%iprhed,
'IPRHED', this%memoryPath)
397 call mem_allocate(this%istageout,
'ISTAGEOUT', this%memoryPath)
398 call mem_allocate(this%ibudgetout,
'IBUDGETOUT', this%memoryPath)
399 call mem_allocate(this%ibudcsv,
'IBUDCSV', this%memoryPath)
400 call mem_allocate(this%ipakcsv,
'IPAKCSV', this%memoryPath)
401 call mem_allocate(this%nlakes,
'NLAKES', this%memoryPath)
402 call mem_allocate(this%noutlets,
'NOUTLETS', this%memoryPath)
403 call mem_allocate(this%ntables,
'NTABLES', this%memoryPath)
404 call mem_allocate(this%convlength,
'CONVLENGTH', this%memoryPath)
405 call mem_allocate(this%convtime,
'CONVTIME', this%memoryPath)
406 call mem_allocate(this%outdmax,
'OUTDMAX', this%memoryPath)
407 call mem_allocate(this%igwhcopt,
'IGWHCOPT', this%memoryPath)
408 call mem_allocate(this%iconvchk,
'ICONVCHK', this%memoryPath)
409 call mem_allocate(this%maxlakit,
'MAXLAKIT', this%memoryPath)
410 call mem_allocate(this%surfdep,
'SURFDEP', this%memoryPath)
411 call mem_allocate(this%dmaxchg,
'DMAXCHG', this%memoryPath)
413 call mem_allocate(this%check_attr,
'CHECK_ATTR', this%memoryPath)
414 call mem_allocate(this%iimplicit,
'IIMPLICIT', this%memoryPath)
415 call mem_allocate(this%iforcefb,
'IFORCEFB', this%memoryPath)
416 call mem_allocate(this%bditems,
'BDITEMS', this%memoryPath)
417 call mem_allocate(this%cbcauxitems,
'CBCAUXITEMS', this%memoryPath)
418 call mem_allocate(this%idense,
'IDENSE', this%memoryPath)
429 this%convlength =
done
437 this%delh =
dp999 * this%dmaxchg
444 end subroutine lak_allocate_scalars
448 subroutine lak_allocate_arrays(this)
451 class(laktype),
intent(inout) :: this
456 call this%BndType%allocate_arrays()
459 allocate (this%clakbudget(this%bditems))
462 this%clakbudget(1) =
' GWF'
463 this%clakbudget(2) =
' RAINFALL'
464 this%clakbudget(3) =
' EVAPORATION'
465 this%clakbudget(4) =
' RUNOFF'
466 this%clakbudget(5) =
' EXT-INFLOW'
467 this%clakbudget(6) =
' WITHDRAWAL'
468 this%clakbudget(7) =
' EXT-OUTFLOW'
469 this%clakbudget(8) =
' STORAGE'
470 this%clakbudget(9) =
' CONSTANT'
471 this%clakbudget(10) =
' FROM-MVR'
472 this%clakbudget(11) =
' TO-MVR'
475 if (this%istageout > 0)
then
476 call mem_allocate(this%dbuff, this%nlakes,
'DBUFF', this%memoryPath)
477 do i = 1, this%nlakes
478 this%dbuff(i) =
dzero
481 call mem_allocate(this%dbuff, 0,
'DBUFF', this%memoryPath)
485 allocate (this%cauxcbc(this%cbcauxitems))
488 call mem_allocate(this%qauxcbc, this%cbcauxitems,
'QAUXCBC', this%memoryPath)
489 do i = 1, this%cbcauxitems
490 this%qauxcbc(i) =
dzero
494 call mem_allocate(this%qleak, this%maxbound,
'QLEAK', this%memoryPath)
495 do i = 1, this%maxbound
496 this%qleak(i) =
dzero
500 if (this%iimplicit /= 0)
then
501 call mem_allocate(this%holdconn, this%maxbound,
'HOLDCONN', this%memoryPath)
502 do i = 1, this%maxbound
503 this%holdconn(i) =
dzero
506 call mem_allocate(this%holdconn, 0,
'HOLDCONN', this%memoryPath)
508 call mem_allocate(this%qsto, this%nlakes,
'QSTO', this%memoryPath)
509 do i = 1, this%nlakes
514 call mem_allocate(this%denseterms, 3, 0,
'DENSETERMS', this%memoryPath)
517 call mem_allocate(this%viscratios, 2, 0,
'VISCRATIOS', this%memoryPath)
518 end subroutine lak_allocate_arrays
522 subroutine lak_read_lakes(this)
528 class(laktype),
intent(inout) :: this
530 character(len=LINELENGTH) :: text
531 character(len=LENBOUNDNAME) :: bndName, bndNameTemp
532 character(len=9) :: cno
533 character(len=50),
dimension(:),
allocatable :: caux
534 integer(I4B) :: ierr, ival
535 logical(LGP) :: isfound, endOfBlock
537 integer(I4B) :: ii, jj
541 integer(I4B) :: nconn
542 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
543 real(DP),
pointer :: bndElem => null()
549 call mem_allocate(this%nlakeconn, this%nlakes,
'NLAKECONN', this%memoryPath)
550 call mem_allocate(this%idxlakeconn, this%nlakes + 1,
'IDXLAKECONN', &
552 call mem_allocate(this%ntabrow, this%nlakes,
'NTABROW', this%memoryPath)
553 call mem_allocate(this%strt, this%nlakes,
'STRT', this%memoryPath)
554 call mem_allocate(this%laketop, this%nlakes,
'LAKETOP', this%memoryPath)
555 call mem_allocate(this%lakebot, this%nlakes,
'LAKEBOT', this%memoryPath)
556 call mem_allocate(this%sareamax, this%nlakes,
'SAREAMAX', this%memoryPath)
557 call mem_allocate(this%stage, this%nlakes,
'STAGE', this%memoryPath)
558 call mem_allocate(this%rainfall, this%nlakes,
'RAINFALL', this%memoryPath)
559 call mem_allocate(this%evaporation, this%nlakes,
'EVAPORATION', &
561 call mem_allocate(this%runoff, this%nlakes,
'RUNOFF', this%memoryPath)
562 call mem_allocate(this%inflow, this%nlakes,
'INFLOW', this%memoryPath)
563 call mem_allocate(this%withdrawal, this%nlakes,
'WITHDRAWAL', this%memoryPath)
564 call mem_allocate(this%lauxvar, this%naux, this%nlakes,
'LAUXVAR', &
566 call mem_allocate(this%avail, this%nlakes,
'AVAIL', this%memoryPath)
567 call mem_allocate(this%lkgwsink, this%nlakes,
'LKGWSINK', this%memoryPath)
568 call mem_allocate(this%ncncvr, this%nlakes,
'NCNCVR', this%memoryPath)
569 call mem_allocate(this%ifallback, this%nlakes,
'IFALLBACK', this%memoryPath)
570 call mem_allocate(this%nstuck, this%nlakes,
'NSTUCK', this%memoryPath)
571 call mem_allocate(this%surfin, this%nlakes,
'SURFIN', this%memoryPath)
572 call mem_allocate(this%surfout, this%nlakes,
'SURFOUT', this%memoryPath)
573 call mem_allocate(this%surfout1, this%nlakes,
'SURFOUT1', this%memoryPath)
574 call mem_allocate(this%precip, this%nlakes,
'PRECIP', this%memoryPath)
575 call mem_allocate(this%precip1, this%nlakes,
'PRECIP1', this%memoryPath)
576 call mem_allocate(this%evap, this%nlakes,
'EVAP', this%memoryPath)
577 call mem_allocate(this%evap1, this%nlakes,
'EVAP1', this%memoryPath)
578 call mem_allocate(this%evapo, this%nlakes,
'EVAPO', this%memoryPath)
579 call mem_allocate(this%withr, this%nlakes,
'WITHR', this%memoryPath)
580 call mem_allocate(this%withr1, this%nlakes,
'WITHR1', this%memoryPath)
581 call mem_allocate(this%flwin, this%nlakes,
'FLWIN', this%memoryPath)
582 call mem_allocate(this%flwiter, this%nlakes,
'FLWITER', this%memoryPath)
583 call mem_allocate(this%flwiter1, this%nlakes,
'FLWITER1', this%memoryPath)
584 call mem_allocate(this%seep, this%nlakes,
'SEEP', this%memoryPath)
585 call mem_allocate(this%seep1, this%nlakes,
'SEEP1', this%memoryPath)
586 call mem_allocate(this%seep0, this%nlakes,
'SEEP0', this%memoryPath)
587 call mem_allocate(this%stageiter, this%nlakes,
'STAGEITER', this%memoryPath)
588 call mem_allocate(this%chterm, this%nlakes,
'CHTERM', this%memoryPath)
594 if (this%iimplicit == 0)
then
595 call mem_allocate(this%iboundpak, this%nlakes,
'IBOUND', this%memoryPath)
596 call mem_allocate(this%xnewpak, this%nlakes,
'XNEWPAK', this%memoryPath)
598 call mem_allocate(this%xoldpak, this%nlakes,
'XOLDPAK', this%memoryPath)
601 call mem_allocate(this%iseepc, this%nlakes,
'ISEEPC', this%memoryPath)
602 call mem_allocate(this%idhc, this%nlakes,
'IDHC', this%memoryPath)
603 call mem_allocate(this%en1, this%nlakes,
'EN1', this%memoryPath)
604 call mem_allocate(this%en2, this%nlakes,
'EN2', this%memoryPath)
605 call mem_allocate(this%r1, this%nlakes,
'R1', this%memoryPath)
606 call mem_allocate(this%r2, this%nlakes,
'R2', this%memoryPath)
607 call mem_allocate(this%dh0, this%nlakes,
'DH0', this%memoryPath)
608 call mem_allocate(this%s0, this%nlakes,
'S0', this%memoryPath)
609 call mem_allocate(this%qgwf0, this%nlakes,
'QGWF0', this%memoryPath)
612 allocate (this%lakename(this%nlakes))
613 allocate (this%status(this%nlakes))
615 do n = 1, this%nlakes
617 this%status(n) =
'ACTIVE'
618 this%laketop(n) = -dep20
619 this%lakebot(n) = dep20
620 this%sareamax(n) = dzero
625 if (this%iimplicit == 0)
then
626 this%iboundpak(n) = 1
627 this%xnewpak(n) = dep20
629 this%xoldpak(n) = dep20
632 this%rainfall(n) = dzero
633 this%evaporation(n) = dzero
634 this%runoff(n) = dzero
635 this%inflow(n) = dzero
636 this%withdrawal(n) = dzero
637 this%ifallback(n) = 0
642 if (this%naux > 0)
then
643 allocate (caux(this%naux))
647 allocate (nboundchk(this%nlakes))
648 do n = 1, this%nlakes
654 call this%parser%GetBlock(
'PACKAGEDATA', isfound, ierr, &
655 supportopenclose=.true.)
659 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
664 call this%parser%GetNextLine(endofblock)
666 n = this%parser%GetInteger()
668 if (n < 1 .or. n > this%nlakes)
then
669 write (
errmsg,
'(a,1x,i0)')
'lakeno MUST BE > 0 and <= ', this%nlakes
675 nboundchk(n) = nboundchk(n) + 1
678 this%strt(n) = this%parser%GetDouble()
681 ival = this%parser%GetInteger()
684 write (
errmsg,
'(a,1x,i0)')
'nlakeconn MUST BE >= 0 for lake ', n
691 if (this%iimplicit /= 0 .and. ival == 0)
then
692 write (
errmsg,
'(a,1x,i0,1x,a)') &
693 'lake', n,
'has no connections; the IMPLICIT option requires &
694 &each lake to have at least one GWF connection.'
699 this%nlakeconn(n) = ival
702 do iaux = 1, this%naux
703 call this%parser%GetString(caux(iaux))
707 write (cno,
'(i9.9)') n
708 bndname =
'Lake'//cno
711 if (this%inamedbound /= 0)
then
712 call this%parser%GetStringCaps(bndnametemp)
713 if (bndnametemp /=
'')
then
714 bndname = bndnametemp
717 this%lakename(n) = bndname
724 bndelem => this%lauxvar(jj, ii)
726 this%packName,
'AUX', &
727 this%tsManager, this%iprpak, &
735 do n = 1, this%nlakes
736 if (nboundchk(n) == 0)
then
737 write (
errmsg,
'(a,1x,i0)')
'NO DATA SPECIFIED FOR LAKE', n
739 else if (nboundchk(n) > 1)
then
740 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
741 'DATA FOR LAKE', n,
'SPECIFIED', nboundchk(n),
'TIMES'
746 write (this%iout,
'(1x,a)')
'END OF '//trim(adjustl(this%text))// &
749 call store_error(
'REQUIRED PACKAGEDATA BLOCK NOT FOUND.')
754 call this%parser%StoreErrorUnit()
758 this%MAXBOUND = nconn
759 write (this%iout,
'(//4x,a,i7)')
'MAXBOUND = ', this%maxbound
762 this%idxlakeconn(1) = 1
763 do n = 1, this%nlakes
764 this%idxlakeconn(n + 1) = this%idxlakeconn(n) + this%nlakeconn(n)
768 if (this%naux > 0)
then
773 deallocate (nboundchk)
774 end subroutine lak_read_lakes
778 subroutine lak_read_lake_connections(this)
782 class(laktype),
intent(inout) :: this
784 character(len=LINELENGTH) :: keyword, cellid
785 integer(I4B) :: ierr, ival
786 logical(LGP) :: isfound, endOfBlock
787 logical(LGP) :: is_lake_bed
791 integer(I4B) :: ipos, ipos0
792 integer(I4B) :: icellid, icellid0
795 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
796 character(len=LENVARNAME) :: ctypenm
799 allocate (nboundchk(this%MAXBOUND))
800 do n = 1, this%MAXBOUND
805 call this%parser%GetBlock(
'CONNECTIONDATA', isfound, ierr, &
806 supportopenclose=.true.)
811 call mem_allocate(this%imap, this%MAXBOUND,
'IMAP', this%memoryPath)
812 call mem_allocate(this%cellid, this%MAXBOUND,
'CELLID', this%memoryPath)
813 call mem_allocate(this%nodesontop, this%MAXBOUND,
'NODESONTOP', &
815 call mem_allocate(this%ictype, this%MAXBOUND,
'ICTYPE', this%memoryPath)
816 call mem_allocate(this%bedleak, this%MAXBOUND,
'BEDLEAK', this%memoryPath)
817 call mem_allocate(this%belev, this%MAXBOUND,
'BELEV', this%memoryPath)
818 call mem_allocate(this%telev, this%MAXBOUND,
'TELEV', this%memoryPath)
819 call mem_allocate(this%connlength, this%MAXBOUND,
'CONNLENGTH', &
821 call mem_allocate(this%connwidth, this%MAXBOUND,
'CONNWIDTH', &
823 call mem_allocate(this%sarea, this%MAXBOUND,
'SAREA', this%memoryPath)
824 call mem_allocate(this%warea, this%MAXBOUND,
'WAREA', this%memoryPath)
825 call mem_allocate(this%satcond, this%MAXBOUND,
'SATCOND', this%memoryPath)
826 call mem_allocate(this%simcond, this%MAXBOUND,
'SIMCOND', this%memoryPath)
827 call mem_allocate(this%simlakgw, this%MAXBOUND,
'SIMLAKGW', this%memoryPath)
830 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
833 call this%parser%GetNextLine(endofblock)
835 n = this%parser%GetInteger()
837 if (n < 1 .or. n > this%nlakes)
then
838 write (
errmsg,
'(a,1x,i0)')
'lakeno MUST BE > 0 and <= ', this%nlakes
844 ival = this%parser%GetInteger()
845 if (ival < 1 .or. ival > this%nlakeconn(n))
then
846 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
847 'iconn FOR LAKE ', n,
'MUST BE > 1 and <= ', this%nlakeconn(n)
853 ipos = this%idxlakeconn(n) + ival - 1
860 nboundchk(ipos) = nboundchk(ipos) + 1
863 call this%parser%GetCellid(this%dis%ndim, cellid)
864 nn = this%dis%noder_from_cellid(cellid, &
865 this%parser%iuactive, this%iout)
869 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
870 'INVALID cellid FOR LAKE ', n,
'connection', j
875 this%cellid(ipos) = nn
876 this%nodesontop(ipos) = nn
879 call this%parser%GetStringCaps(keyword)
880 select case (keyword)
882 this%ictype(ipos) = 0
884 this%ictype(ipos) = 1
886 this%ictype(ipos) = 2
888 this%ictype(ipos) = 3
890 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,a,a)') &
891 'UNKNOWN ctype FOR LAKE ', n,
'connection', j, &
892 '(', trim(keyword),
')'
895 write (ctypenm,
'(a16)') keyword
899 call this%parser%GetStringCaps(keyword)
900 select case (keyword)
902 is_lake_bed = .false.
903 this%bedleak(ipos) = dnodata
906 write (
warnmsg,
'(2(a,1x,i0,1x),a,1pe8.1,a)') &
907 'BEDLEAK for connection', j,
'in lake', n,
'is specified to '// &
908 'be NONE. Lake connections where the lake-GWF connection '// &
909 'conductance is solely a function of aquifer properties '// &
910 'in the connected GWF cell should be specified with a '// &
911 'DNODATA (', dnodata,
') value.'
914 call deprecation_warning(
'CONNECTIONDATA',
'bedleak=NONE',
'6.4.3', &
915 warnmsg, this%parser%GetUnit())
917 read (keyword, *) rval
919 is_lake_bed = .false.
923 this%bedleak(ipos) = rval
926 if (is_lake_bed .and. this%bedleak(ipos) < dzero)
then
927 write (
errmsg,
'(a,1x,i0,1x,a)')
'bedleak FOR LAKE ', n,
'MUST BE >= 0'
932 this%belev(ipos) = this%parser%GetDouble()
935 this%telev(ipos) = this%parser%GetDouble()
938 rval = this%parser%GetDouble()
939 if (rval <= dzero)
then
940 if (this%ictype(ipos) == 1 .or. this%ictype(ipos) == 2 .or. &
941 this%ictype(ipos) == 3)
then
942 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,a,1x,a)') &
943 'connection length (connlen) FOR LAKE ', n, &
944 ', CONNECTION NO.', j,
', MUST BE > 0 FOR SPECIFIED ', &
945 'connection type (ctype)', ctypenm
951 this%connlength(ipos) = rval
954 rval = this%parser%GetDouble()
955 if (rval < dzero)
then
956 if (this%ictype(ipos) == 1)
then
957 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
958 'cell width (connwidth) FOR LAKE ', n, &
959 ' HORIZONTAL CONNECTION ', j,
'MUST BE >= 0'
965 this%connwidth(ipos) = rval
967 write (this%iout,
'(1x,a)') &
968 'END OF '//trim(adjustl(this%text))//
' CONNECTIONDATA'
970 call store_error(
'REQUIRED CONNECTIONDATA BLOCK NOT FOUND.')
975 call this%parser%StoreErrorUnit()
979 do n = 1, this%nlakes
981 do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
982 if (this%ictype(ipos) /= 2 .and. this%ictype(ipos) /= 3) cycle
985 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
986 'nlakeconn FOR LAKE', n,
'EMBEDDED CONNECTION', j,
' EXCEEDS 1.'
993 do n = 1, this%nlakes
994 ipos0 = this%idxlakeconn(n)
995 icellid0 = this%cellid(ipos0)
996 if (this%ictype(ipos0) /= 2 .and. this%ictype(ipos0) /= 3) cycle
997 do nn = 1, this%nlakes
1000 do ipos = this%idxlakeconn(nn), this%idxlakeconn(nn + 1) - 1
1002 icellid = this%cellid(ipos)
1003 if (icellid == icellid0)
then
1004 if (this%ictype(ipos) == 0)
then
1005 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
1006 'EMBEDDED LAKE', n, &
1007 'CANNOT COINCIDE WITH VERTICAL CONNECTION', j, &
1017 do n = 1, this%nlakes
1019 do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1021 nn = this%cellid(ipos)
1022 top = this%dis%top(nn)
1023 bot = this%dis%bot(nn)
1025 if (this%ictype(ipos) == 0)
then
1026 this%telev(ipos) = top + this%surfdep
1027 this%belev(ipos) = top
1028 this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
1030 else if (this%ictype(ipos) == 1)
then
1031 if (this%belev(ipos) == this%telev(ipos))
then
1032 this%telev(ipos) = top
1033 this%belev(ipos) = bot
1035 if (this%belev(ipos) >= this%telev(ipos))
then
1036 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1037 'telev FOR LAKE ', n,
' HORIZONTAL CONNECTION ', j, &
1040 else if (this%belev(ipos) < bot)
then
1041 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
1042 'belev FOR LAKE ', n,
' HORIZONTAL CONNECTION ', j, &
1043 'MUST BE >= cell bottom (', bot,
')'
1045 else if (this%telev(ipos) > top)
then
1046 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
1047 'telev FOR LAKE ', n,
' HORIZONTAL CONNECTION ', j, &
1048 'MUST BE <= cell top (', top,
')'
1052 this%laketop(n) = max(this%telev(ipos), this%laketop(n))
1053 this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
1055 else if (this%ictype(ipos) == 2 .or. this%ictype(ipos) == 3)
then
1056 this%telev(ipos) = top
1057 this%belev(ipos) = bot
1058 this%lakebot(n) = bot
1062 if (nboundchk(ipos) == 0)
then
1063 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1064 'NO DATA SPECIFIED FOR LAKE', n,
'CONNECTION', j
1066 else if (nboundchk(ipos) > 1)
then
1067 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
1068 'DATA FOR LAKE', n,
'CONNECTION', j, &
1069 'SPECIFIED', nboundchk(ipos),
'TIMES'
1075 if (this%laketop(n) == -dep20)
then
1076 this%laketop(n) = this%lakebot(n) + 100.
1081 deallocate (nboundchk)
1085 call this%parser%StoreErrorUnit()
1087 end subroutine lak_read_lake_connections
1091 subroutine lak_read_tables(this)
1095 class(laktype),
intent(inout) :: this
1097 type(laktabtype),
dimension(:),
allocatable :: laketables
1098 character(len=LINELENGTH) :: line
1099 character(len=LINELENGTH) :: keyword
1100 integer(I4B) :: ierr
1101 logical(LGP) :: isfound, endOfBlock
1103 integer(I4B) :: iconn
1104 integer(I4B) :: ntabs
1105 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
1108 if (this%ntables < 1)
return
1111 allocate (nboundchk(this%nlakes))
1112 do n = 1, this%nlakes
1117 allocate (laketables(this%nlakes))
1120 call this%parser%GetBlock(
'TABLES', isfound, ierr, &
1121 supportopenclose=.true.)
1127 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
1130 call this%parser%GetNextLine(endofblock)
1131 if (endofblock)
exit
1132 n = this%parser%GetInteger()
1134 if (n < 1 .or. n > this%nlakes)
then
1135 write (
errmsg,
'(a,1x,i0)')
'lakeno MUST BE > 0 and <= ', this%nlakes
1142 nboundchk(n) = nboundchk(n) + 1
1145 call this%parser%GetStringCaps(keyword)
1146 select case (keyword)
1148 call this%parser%GetStringCaps(keyword)
1149 if (trim(adjustl(keyword)) /=
'FILEIN')
then
1150 errmsg =
'TAB6 keyword must be followed by "FILEIN" '// &
1155 call this%parser%GetString(line)
1156 call this%lak_read_table(n, line, laketables(n))
1158 write (
errmsg,
'(a,1x,i0,1x,a)') &
1159 'LAKE TABLE ENTRY for LAKE ', n,
'MUST INCLUDE TAB6 KEYWORD'
1165 write (this%iout,
'(1x,a)') &
1166 'END OF '//trim(adjustl(this%text))//
' LAKE_TABLES'
1169 if (ntabs < this%ntables)
then
1170 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1171 'TABLE DATA ARE SPECIFIED', ntabs, &
1172 'TIMES BUT NTABLES IS SET TO', this%ntables
1175 do n = 1, this%nlakes
1176 if (this%ntabrow(n) > 0 .and. nboundchk(n) > 1)
then
1177 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1178 'TABLE DATA FOR LAKE', n,
'SPECIFIED', nboundchk(n),
'TIMES'
1183 call store_error(
'REQUIRED TABLES BLOCK NOT FOUND.')
1187 deallocate (nboundchk)
1191 call this%parser%StoreErrorUnit()
1195 call this%laktables_to_vectors(laketables)
1198 do n = 1, this%nlakes
1199 if (this%ntabrow(n) > 0)
then
1200 deallocate (laketables(n)%tabstage)
1201 deallocate (laketables(n)%tabvolume)
1202 deallocate (laketables(n)%tabsarea)
1203 iconn = this%idxlakeconn(n)
1204 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1205 deallocate (laketables(n)%tabwarea)
1209 deallocate (laketables)
1210 end subroutine lak_read_tables
1215 subroutine laktables_to_vectors(this, laketables)
1216 class(laktype),
intent(inout) :: this
1217 type(laktabtype),
intent(in),
dimension(:),
contiguous :: laketables
1219 integer(I4B) :: ntabrows
1221 integer(I4B) :: ipos
1222 integer(I4B) :: iconn
1225 call mem_allocate(this%ialaktab, this%nlakes + 1,
'IALAKTAB', this%memoryPath)
1228 this%ialaktab(1) = 1
1229 do n = 1, this%nlakes
1231 this%ialaktab(n + 1) = this%ialaktab(n) + this%ntabrow(n)
1235 ntabrows = this%ialaktab(this%nlakes + 1) - 1
1236 call mem_allocate(this%tabstage, ntabrows,
'TABSTAGE', this%memoryPath)
1237 call mem_allocate(this%tabvolume, ntabrows,
'TABVOLUME', this%memoryPath)
1238 call mem_allocate(this%tabsarea, ntabrows,
'TABSAREA', this%memoryPath)
1239 call mem_allocate(this%tabwarea, ntabrows,
'TABWAREA', this%memoryPath)
1242 do n = 1, this%nlakes
1244 do ipos = this%ialaktab(n), this%ialaktab(n + 1) - 1
1245 this%tabstage(ipos) = laketables(n)%tabstage(j)
1246 this%tabvolume(ipos) = laketables(n)%tabvolume(j)
1247 this%tabsarea(ipos) = laketables(n)%tabsarea(j)
1248 iconn = this%idxlakeconn(n)
1249 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1252 this%tabwarea(ipos) = laketables(n)%tabwarea(j)
1254 this%tabwarea(ipos) =
dzero
1259 end subroutine laktables_to_vectors
1263 subroutine lak_read_table(this, ilak, filename, laketable)
1268 class(laktype),
intent(inout) :: this
1269 integer(I4B),
intent(in) :: ilak
1270 character(len=*),
intent(in) :: filename
1271 type(laktabtype),
intent(inout) :: laketable
1273 character(len=LINELENGTH) :: keyword
1274 integer(I4B) :: ierr
1275 logical(LGP) :: isfound, endOfBlock
1278 integer(I4B) :: ipos
1280 integer(I4B) :: jmin
1281 integer(I4B) :: iconn
1289 character(len=*),
parameter :: fmttaberr = &
1290 &
'(a,1x,i0,1x,a,1x,g15.6,1x,a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.6,1x,a)'
1298 call openfile(iu, this%iout, filename,
'LAKE TABLE')
1299 call parser%Initialize(iu, this%iout)
1302 call parser%GetBlock(
'DIMENSIONS', isfound, ierr, supportopenclose=.true.)
1307 if (this%iprpak /= 0)
then
1308 write (this%iout,
'(/1x,a)') &
1309 'PROCESSING '//trim(adjustl(this%text))//
' DIMENSIONS'
1312 call parser%GetNextLine(endofblock)
1313 if (endofblock)
exit
1314 call parser%GetStringCaps(keyword)
1315 select case (keyword)
1317 n = parser%GetInteger()
1320 write (
errmsg,
'(a)')
'LAKE TABLE NROW MUST BE > 0'
1324 j = parser%GetInteger()
1326 if (this%ictype(ilak) == 2 .or. this%ictype(ilak) == 3)
then
1332 write (
errmsg,
'(a,1x,i0)')
'LAKE TABLE NCOL MUST BE >= ', jmin
1337 write (
errmsg,
'(a,a)') &
1338 'UNKNOWN '//trim(this%text)//
' DIMENSIONS KEYWORD: ', trim(keyword)
1342 if (this%iprpak /= 0)
then
1343 write (this%iout,
'(1x,a)') &
1344 'END OF '//trim(adjustl(this%text))//
' DIMENSIONS'
1347 call store_error(
'REQUIRED DIMENSIONS BLOCK NOT FOUND.')
1353 'NROW NOT SPECIFIED IN THE LAKE TABLE DIMENSIONS BLOCK'
1358 'NCOL NOT SPECIFIED IN THE LAKE TABLE DIMENSIONS BLOCK'
1367 this%ntabrow(ilak) = n
1368 allocate (laketable%tabstage(n))
1369 allocate (laketable%tabvolume(n))
1370 allocate (laketable%tabsarea(n))
1371 ipos = this%idxlakeconn(ilak)
1372 if (this%ictype(ipos) == 2 .or. this%ictype(ipos) == 3)
then
1373 allocate (laketable%tabwarea(n))
1377 call parser%GetBlock(
'TABLE', isfound, ierr, supportopenclose=.true.)
1383 if (this%iprpak /= 0)
then
1384 write (this%iout,
'(/1x,a)') &
1385 'PROCESSING '//trim(adjustl(this%text))//
' TABLE'
1387 iconn = this%idxlakeconn(ilak)
1390 call parser%GetNextLine(endofblock)
1391 if (endofblock)
exit
1393 if (ipos > this%ntabrow(ilak))
then
1396 laketable%tabstage(ipos) = parser%GetDouble()
1397 laketable%tabvolume(ipos) = parser%GetDouble()
1398 laketable%tabsarea(ipos) = parser%GetDouble()
1399 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1400 laketable%tabwarea(ipos) = parser%GetDouble()
1402 end do readtabledata
1404 if (this%iprpak /= 0)
then
1405 write (this%iout,
'(1x,a)') &
1406 'END OF '//trim(adjustl(this%text))//
' TABLE'
1409 call store_error(
'REQUIRED TABLE BLOCK NOT FOUND.')
1413 if (ipos /= this%ntabrow(ilak))
then
1414 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1415 'NROW SET TO', this%ntabrow(ilak),
'BUT', ipos,
'ROWS WERE READ'
1420 iconn = this%idxlakeconn(ilak)
1421 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1422 do n = 1, this%ntabrow(ilak)
1423 vol = laketable%tabvolume(n)
1424 sa = laketable%tabsarea(n)
1425 wa = laketable%tabwarea(n)
1428 if (vol > dzero)
exit
1430 this%lakebot(ilak) = laketable%tabstage(n)
1431 this%belev(ilak) = laketable%tabstage(n)
1434 n = this%ntabrow(ilak)
1435 this%sareamax(ilak) = laketable%tabsarea(n)
1439 do n = 2, this%ntabrow(ilak)
1440 v = laketable%tabstage(n)
1441 v0 = laketable%tabstage(n - 1)
1443 write (
errmsg, fmttaberr) &
1444 'TABLE STAGE ENTRY', n,
'(', laketable%tabstage(n),
') FOR LAKE ', &
1445 ilak,
'MUST BE GREATER THAN THE PREVIOUS STAGE ENTRY', &
1446 n - 1,
'(', laketable%tabstage(n - 1),
')'
1449 v = laketable%tabvolume(n)
1450 v0 = laketable%tabvolume(n - 1)
1452 write (
errmsg, fmttaberr) &
1453 'TABLE VOLUME ENTRY', n,
'(', laketable%tabvolume(n), &
1455 ilak,
'MUST BE GREATER THAN THE PREVIOUS VOLUME ENTRY', &
1456 n - 1,
'(', laketable%tabvolume(n - 1),
')'
1459 v = laketable%tabsarea(n)
1460 v0 = laketable%tabsarea(n - 1)
1462 write (
errmsg, fmttaberr) &
1463 'TABLE SURFACE AREA ENTRY', n,
'(', &
1464 laketable%tabsarea(n),
') FOR LAKE ', ilak, &
1465 'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS SURFACE AREA ENTRY', &
1466 n - 1,
'(', laketable%tabsarea(n - 1),
')'
1469 iconn = this%idxlakeconn(ilak)
1470 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
1471 v = laketable%tabwarea(n)
1472 v0 = laketable%tabwarea(n - 1)
1474 write (
errmsg, fmttaberr) &
1475 'TABLE EXCHANGE AREA ENTRY', n,
'(', &
1476 laketable%tabwarea(n),
') FOR LAKE ', ilak, &
1477 'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS EXCHANGE AREA '// &
1478 'ENTRY', n - 1,
'(', laketable%tabwarea(n - 1),
')'
1487 call parser%StoreErrorUnit()
1492 end subroutine lak_read_table
1496 subroutine lak_read_outlets(this)
1501 class(laktype),
intent(inout) :: this
1503 character(len=LINELENGTH) :: text, keyword
1504 character(len=LENBOUNDNAME) :: bndName
1505 character(len=9) :: citem
1506 integer(I4B) :: ierr, ival
1507 logical(LGP) :: isfound, endOfBlock
1510 integer(I4B),
dimension(:),
pointer,
contiguous :: nboundchk
1511 real(DP),
pointer :: bndElem => null()
1514 call this%parser%GetBlock(
'OUTLETS', isfound, ierr, &
1515 supportopenclose=.true., blockrequired=.false.)
1519 if (this%noutlets > 0)
then
1522 allocate (nboundchk(this%noutlets))
1523 do n = 1, this%noutlets
1528 call mem_allocate(this%lakein, this%NOUTLETS,
'LAKEIN', this%memoryPath)
1529 call mem_allocate(this%lakeout, this%NOUTLETS,
'LAKEOUT', this%memoryPath)
1530 call mem_allocate(this%iouttype, this%NOUTLETS,
'IOUTTYPE', &
1532 call mem_allocate(this%outrate, this%NOUTLETS,
'OUTRATE', this%memoryPath)
1533 call mem_allocate(this%outinvert, this%NOUTLETS,
'OUTINVERT', &
1535 call mem_allocate(this%outwidth, this%NOUTLETS,
'OUTWIDTH', &
1537 call mem_allocate(this%outrough, this%NOUTLETS,
'OUTROUGH', &
1539 call mem_allocate(this%outslope, this%NOUTLETS,
'OUTSLOPE', &
1541 call mem_allocate(this%simoutrate, this%NOUTLETS,
'SIMOUTRATE', &
1545 do n = 1, this%noutlets
1546 this%outrate(n) = dzero
1550 write (this%iout,
'(/1x,a)') &
1551 'PROCESSING '//trim(adjustl(this%text))//
' OUTLETS'
1553 call this%parser%GetNextLine(endofblock)
1554 if (endofblock)
exit
1555 n = this%parser%GetInteger()
1557 if (n < 1 .or. n > this%noutlets)
then
1558 write (
errmsg,
'(a,1x,i0)') &
1559 'outletno MUST BE > 0 and <= ', this%noutlets
1565 nboundchk(n) = nboundchk(n) + 1
1568 ival = this%parser%GetInteger()
1569 if (ival < 1 .or. ival > this%nlakes)
then
1570 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1571 'lakein FOR OUTLET ', n,
'MUST BE > 0 and <= ', this%nlakes
1575 this%lakein(n) = ival
1578 ival = this%parser%GetInteger()
1579 if (ival < 0 .or. ival > this%nlakes)
then
1580 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0)') &
1581 'lakeout FOR OUTLET ', n,
'MUST BE >= 0 and <= ', this%nlakes
1585 this%lakeout(n) = ival
1588 call this%parser%GetStringCaps(keyword)
1589 select case (keyword)
1591 this%iouttype(n) = 0
1593 this%iouttype(n) = 1
1595 this%iouttype(n) = 2
1597 write (
errmsg,
'(a,1x,i0,1x,a,a,a)') &
1598 'UNKNOWN couttype FOR OUTLET ', n,
'(', trim(keyword),
')'
1604 write (citem,
'(i9.9)') n
1605 bndname =
'OUTLET'//citem
1611 call this%parser%GetString(text)
1612 bndelem => this%outinvert(n)
1614 this%packName,
'BND', &
1615 this%tsManager, this%iprpak, &
1619 call this%parser%GetString(text)
1620 bndelem => this%outwidth(n)
1622 this%packName,
'BND', &
1623 this%tsManager, this%iprpak,
'WIDTH')
1626 call this%parser%GetString(text)
1627 bndelem => this%outrough(n)
1629 this%packName,
'BND', &
1630 this%tsManager, this%iprpak,
'ROUGH')
1633 call this%parser%GetString(text)
1634 bndelem => this%outslope(n)
1636 this%packName,
'BND', &
1637 this%tsManager, this%iprpak,
'SLOPE')
1639 write (this%iout,
'(1x,a)')
'END OF '//trim(adjustl(this%text))// &
1643 do n = 1, this%noutlets
1644 if (nboundchk(n) == 0)
then
1645 write (
errmsg,
'(a,1x,i0)')
'NO DATA SPECIFIED FOR OUTLET', n
1647 else if (nboundchk(n) > 1)
then
1648 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
1649 'DATA FOR OUTLET', n,
'SPECIFIED', nboundchk(n),
'TIMES'
1655 deallocate (nboundchk)
1657 write (
errmsg,
'(a,1x,a)') &
1658 'AN OUTLETS BLOCK SHOULD NOT BE SPECIFIED IF NOUTLETS IS NOT', &
1659 'SPECIFIED OR IS SPECIFIED TO BE 0.'
1664 if (this%noutlets > 0)
then
1665 call store_error(
'REQUIRED OUTLETS BLOCK NOT FOUND.')
1672 call this%parser%StoreErrorUnit()
1674 end subroutine lak_read_outlets
1678 subroutine lak_read_dimensions(this)
1682 class(laktype),
intent(inout) :: this
1684 character(len=LINELENGTH) :: keyword
1685 integer(I4B) :: ierr
1686 logical(LGP) :: isfound, endOfBlock
1693 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
1694 supportopenclose=.true.)
1698 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text))// &
1701 call this%parser%GetNextLine(endofblock)
1702 if (endofblock)
exit
1703 call this%parser%GetStringCaps(keyword)
1704 select case (keyword)
1706 this%nlakes = this%parser%GetInteger()
1707 write (this%iout,
'(4x,a,i7)')
'NLAKES = ', this%nlakes
1709 this%noutlets = this%parser%GetInteger()
1710 write (this%iout,
'(4x,a,i7)')
'NOUTLETS = ', this%noutlets
1712 this%ntables = this%parser%GetInteger()
1713 write (this%iout,
'(4x,a,i7)')
'NTABLES = ', this%ntables
1715 write (
errmsg,
'(a,a)') &
1716 'UNKNOWN '//trim(this%text)//
' DIMENSION: ', trim(keyword)
1720 write (this%iout,
'(1x,a)') &
1721 'END OF '//trim(adjustl(this%text))//
' DIMENSIONS'
1723 call store_error(
'REQUIRED DIMENSIONS BLOCK NOT FOUND.')
1726 if (this%nlakes < 0)
then
1728 'NLAKES WAS NOT SPECIFIED OR WAS SPECIFIED INCORRECTLY.'
1734 call this%parser%StoreErrorUnit()
1739 if (this%iimplicit /= 0)
then
1740 this%npakeq = this%nlakes
1749 call this%lak_read_lakes()
1752 call this%lak_read_lake_connections()
1755 call this%lak_read_tables()
1758 call this%lak_read_outlets()
1762 call this%define_listlabel()
1765 call this%lak_setup_budobj()
1768 call this%lak_setup_tableobj()
1769 end subroutine lak_read_dimensions
1773 subroutine lak_read_initial_attr(this)
1779 class(laktype),
intent(inout) :: this
1781 character(len=LINELENGTH) :: text
1782 integer(I4B) :: j, jj, n
1799 real(DP),
allocatable,
dimension(:) :: clb, caq
1800 character(len=14) :: cbedleak
1801 character(len=14) :: cbedcond
1802 character(len=10),
dimension(0:3) :: ctype
1803 character(len=15) :: nodestr
1804 real(DP),
pointer :: bndElem => null()
1806 data ctype(0)/
'VERTICAL '/
1807 data ctype(1)/
'HORIZONTAL'/
1808 data ctype(2)/
'EMBEDDEDH '/
1809 data ctype(3)/
'EMBEDDEDV '/
1812 do n = 1, this%nlakes
1813 this%xnewpak(n) = this%strt(n)
1814 write (text,
'(g15.7)') this%strt(n)
1816 bndelem => this%stage(n)
1818 'BND', this%tsManager, this%iprpak, &
1823 do n = 1, this%nlakes
1824 if (this%status(n) ==
'CONSTANT')
then
1825 this%iboundpak(n) = -1
1826 else if (this%status(n) ==
'INACTIVE')
then
1827 this%iboundpak(n) = 0
1828 else if (this%status(n) ==
'ACTIVE ')
then
1829 this%iboundpak(n) = 1
1834 if (this%inamedbound /= 0)
then
1835 do n = 1, this%nlakes
1836 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1837 this%boundname(j) = this%lakename(n)
1843 call this%copy_boundname()
1853 allocate (clb(this%MAXBOUND))
1854 allocate (caq(this%MAXBOUND))
1857 do n = 1, this%nlakes
1858 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1860 top = this%dis%top(nn)
1861 bot = this%dis%bot(nn)
1863 if (this%ictype(j) == 0)
then
1864 area = this%dis%area(nn)
1865 this%sarea(j) = area
1866 this%warea(j) = area
1867 this%sareamax(n) = this%sareamax(n) + area
1868 if (this%gwfik33 == 0)
then
1873 length = dhalf * (top - bot)
1875 else if (this%ictype(j) == 1)
then
1876 area = (this%telev(j) - this%belev(j)) * this%connwidth(j)
1879 if (top == this%telev(j) .and. bot == this%belev(j))
then
1880 if (this%icelltype(nn) == 0)
then
1881 area = this%gwfsat(nn) * (top - bot) * this%connwidth(j)
1884 this%sarea(j) = dzero
1885 this%warea(j) = area
1886 this%sareamax(n) = this%sareamax(n) + dzero
1888 length = this%connlength(j)
1890 else if (this%ictype(j) == 2)
then
1892 this%sarea(j) = dzero
1893 this%warea(j) = area
1894 this%sareamax(n) = this%sareamax(n) + dzero
1896 length = this%connlength(j)
1898 else if (this%ictype(j) == 3)
then
1900 this%sarea(j) = dzero
1901 this%warea(j) = area
1902 this%sareamax(n) = this%sareamax(n) + dzero
1903 if (this%gwfik33 == 0)
then
1908 length = this%connlength(j)
1910 if (
is_close(this%bedleak(j), dnodata))
then
1912 else if (this%bedleak(j) > dzero)
then
1913 clb(j) = done / this%bedleak(j)
1922 if (
is_close(this%bedleak(j), dnodata))
then
1923 this%satcond(j) = area / caq(j)
1924 else if (clb(j) * caq(j) > dzero)
then
1925 this%satcond(j) = area / (clb(j) + caq(j))
1927 this%satcond(j) = dzero
1933 if (this%iprpak > 0)
then
1934 write (this%iout,
'(//,29x,a,/)') &
1935 'INTERFACE CONDUCTANCE BETWEEN LAKE AND AQUIFER CELLS'
1936 write (this%iout,
'(1x,a)') &
1937 &
' LAKE CONNECTION CONNECTION LAKEBED'// &
1938 &
' C O N D U C T A N C E S '
1939 write (this%iout,
'(1x,a)') &
1940 &
' NUMBER NUMBER CELLID DIRECTION LEAKANCE'// &
1941 &
' LAKEBED AQUIFER COMBINED'
1942 write (this%iout,
"(1x,108('-'))")
1943 do n = 1, this%nlakes
1945 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1948 if (this%ictype(j) == 1)
then
1949 fact = this%telev(j) - this%belev(j)
1950 if (abs(fact) > dzero)
then
1955 area = this%warea(j)
1957 if (
is_close(clb(j), dnodata))
then
1960 else if (clb(j) > dzero)
then
1961 c1 = area * fact / clb(j)
1962 write (cbedleak,
'(g14.5)') this%bedleak(j)
1963 write (cbedcond,
'(g14.5)') c1
1965 write (cbedleak,
'(g14.5)') c1
1966 write (cbedcond,
'(g14.5)') c1
1969 if (caq(j) > dzero)
then
1970 c2 = area * fact / caq(j)
1972 call this%dis%noder_to_string(nn, nodestr)
1974 '(1x,i10,1x,i10,1x,a15,1x,a10,2(1x,a14),2(1x,g14.5))') &
1975 n, idx, nodestr, ctype(this%ictype(j)), cbedleak, &
1976 cbedcond, c2, this%satcond(j) * fact
1979 write (this%iout,
"(1x,108('-'))")
1980 write (this%iout,
'(1x,a)') &
1981 'IF VERTICAL CONNECTION, CONDUCTANCE (L^2/T) IS &
1982 &BETWEEN AQUIFER CELL AND OVERLYING LAKE CELL.'
1983 write (this%iout,
'(1x,a)') &
1984 'IF HORIZONTAL CONNECTION, CONDUCTANCES ARE PER &
1985 &UNIT SATURATED THICKNESS (L/T).'
1986 write (this%iout,
'(1x,a)') &
1987 'IF EMBEDDED CONNECTION, CONDUCTANCES ARE PER &
1988 &UNIT EXCHANGE AREA (1/T).'
1993 do n = 1, this%nlakes
1994 write (this%iout,
'(//1x,a,1x,i10)')
'STAGE/VOLUME RELATION FOR LAKE ', n
1995 write (this%iout,
'(/1x,5(a14))')
' STAGE',
' SURFACE AREA', &
1996 &
' WETTED AREA',
' CONDUCTANCE', &
1998 write (this%iout,
"(1x,70('-'))")
1999 dx = (this%laketop(n) - this%lakebot(n)) / 150.
2002 call this%lak_calculate_conductance(n, s, c)
2003 call this%lak_calculate_sarea(n, s, sa)
2004 call this%lak_calculate_warea(n, s, wa, s)
2005 call this%lak_calculate_vol(n, s, v)
2006 write (this%iout,
'(1x,5(E14.5))') s, sa, wa, c, v
2009 write (this%iout,
"(1x,70('-'))")
2011 write (this%iout,
'(//1x,a,1x,i10)')
'STAGE/VOLUME RELATION FOR LAKE ', n
2012 write (this%iout,
'(/1x,4(a14))')
' ',
' ', &
2013 &
' CALCULATED',
' STAGE'
2014 write (this%iout,
'(1x,4(a14))')
' STAGE',
' VOLUME', &
2015 &
' STAGE',
' DIFFERENCE'
2016 write (this%iout,
"(1x,56('-'))")
2017 s = this%lakebot(n) - dx
2019 call this%lak_calculate_vol(n, s, v)
2020 call this%lak_vol2stage(n, v, c)
2021 write (this%iout,
'(1x,4(E14.5))') s, v, c, s - c
2024 write (this%iout,
"(1x,56('-'))")
2029 this%gwfk11 => null()
2030 this%gwfk33 => null()
2031 this%gwfsat => null()
2032 this%gwfik33 => null()
2037 end subroutine lak_read_initial_attr
2043 subroutine lak_linear_interpolation(this, n, x, y, z, v)
2045 class(laktype),
intent(inout) :: this
2046 integer(I4B),
intent(in) :: n
2047 real(DP),
dimension(n),
intent(in) :: x
2048 real(DP),
dimension(n),
intent(in) :: y
2049 real(DP),
intent(in) :: z
2050 real(DP),
intent(inout) :: v
2053 real(DP) :: dx, dydx
2061 else if (z > x(n))
then
2062 dx = x(n) - x(n - 1)
2064 if (abs(dx) >
dzero)
then
2065 dydx = (y(n) - y(n - 1)) / dx
2068 v = y(n) + dydx * dx
2072 dx = x(i) - x(i - 1)
2074 if (z >= x(i - 1) .and. z <= x(i))
then
2075 if (abs(dx) >
dzero)
then
2076 dydx = (y(i) - y(i - 1)) / dx
2079 v = y(i - 1) + dydx * dx
2084 end subroutine lak_linear_interpolation
2088 subroutine lak_calculate_sarea(this, ilak, stage, sarea)
2090 class(laktype),
intent(inout) :: this
2091 integer(I4B),
intent(in) :: ilak
2092 real(DP),
intent(in) :: stage
2093 real(DP),
intent(inout) :: sarea
2096 integer(I4B) :: ifirst
2097 integer(I4B) :: ilast
2104 i = this%ntabrow(ilak)
2106 ifirst = this%ialaktab(ilak)
2107 ilast = this%ialaktab(ilak + 1) - 1
2108 if (stage <= this%tabstage(ifirst))
then
2109 sarea = this%tabsarea(ifirst)
2110 else if (stage >= this%tabstage(ilast))
then
2111 sarea = this%tabsarea(ilast)
2113 call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2114 this%tabsarea(ifirst:ilast), &
2118 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2119 topl = this%telev(i)
2120 botl = this%belev(i)
2122 sa = sat * this%sarea(i)
2126 end subroutine lak_calculate_sarea
2130 subroutine lak_calculate_warea(this, ilak, stage, warea, hin)
2132 class(laktype),
intent(inout) :: this
2133 integer(I4B),
intent(in) :: ilak
2134 real(DP),
intent(in) :: stage
2135 real(DP),
intent(inout) :: warea
2136 real(DP),
optional,
intent(inout) :: hin
2139 integer(I4B) :: igwfnode
2144 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2145 if (
present(hin))
then
2148 igwfnode = this%cellid(i)
2149 head = this%xnew(igwfnode)
2151 call this%lak_calculate_conn_warea(ilak, i, stage, head, wa)
2154 end subroutine lak_calculate_warea
2158 subroutine lak_calculate_conn_warea(this, ilak, iconn, stage, head, wa)
2160 class(laktype),
intent(inout) :: this
2161 integer(I4B),
intent(in) :: ilak
2162 integer(I4B),
intent(in) :: iconn
2163 real(DP),
intent(in) :: stage
2164 real(DP),
intent(in) :: head
2165 real(DP),
intent(inout) :: wa
2168 integer(I4B) :: ifirst
2169 integer(I4B) :: ilast
2170 integer(I4B) :: node
2177 topl = this%telev(iconn)
2178 botl = this%belev(iconn)
2179 call this%lak_calculate_cond_head(iconn, stage, head, vv)
2180 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
2181 if (vv > topl) vv = topl
2182 i = this%ntabrow(ilak)
2183 ifirst = this%ialaktab(ilak)
2184 ilast = this%ialaktab(ilak + 1) - 1
2185 if (vv <= this%tabstage(ifirst))
then
2186 wa = this%tabwarea(ifirst)
2187 else if (vv >= this%tabstage(ilast))
then
2188 wa = this%tabwarea(ilast)
2190 call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2191 this%tabwarea(ifirst:ilast), &
2195 node = this%cellid(iconn)
2197 if (this%icelltype(node) == 0)
then
2203 wa = sat * this%warea(iconn)
2205 end subroutine lak_calculate_conn_warea
2209 subroutine lak_calculate_vol(this, ilak, stage, volume)
2211 class(laktype),
intent(inout) :: this
2212 integer(I4B),
intent(in) :: ilak
2213 real(DP),
intent(in) :: stage
2214 real(DP),
intent(inout) :: volume
2217 integer(I4B) :: ifirst
2218 integer(I4B) :: ilast
2227 i = this%ntabrow(ilak)
2229 ifirst = this%ialaktab(ilak)
2230 ilast = this%ialaktab(ilak + 1) - 1
2231 if (stage <= this%tabstage(ifirst))
then
2232 volume = this%tabvolume(ifirst)
2233 else if (stage >= this%tabstage(ilast))
then
2234 ds = stage - this%tabstage(ilast)
2235 sa = this%tabsarea(ilast)
2236 volume = this%tabvolume(ilast) + ds * sa
2238 call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2239 this%tabvolume(ifirst:ilast), &
2243 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2244 topl = this%telev(i)
2245 botl = this%belev(i)
2247 sa = sat * this%sarea(i)
2248 if (stage < botl)
then
2250 else if (stage > botl .and. stage < topl)
then
2251 v = sa * (stage - botl)
2253 v = sa * (topl - botl) + sa * (stage - topl)
2258 end subroutine lak_calculate_vol
2262 subroutine lak_calculate_conductance(this, ilak, stage, conductance)
2264 class(laktype),
intent(inout) :: this
2265 integer(I4B),
intent(in) :: ilak
2266 real(DP),
intent(in) :: stage
2267 real(DP),
intent(inout) :: conductance
2273 do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2274 call this%lak_calculate_conn_conductance(ilak, i, stage, stage, c)
2275 conductance = conductance + c
2277 end subroutine lak_calculate_conductance
2283 subroutine lak_calculate_cond_head(this, iconn, stage, head, vv)
2285 class(laktype),
intent(inout) :: this
2286 integer(I4B),
intent(in) :: iconn
2287 real(DP),
intent(in) :: stage
2288 real(DP),
intent(in) :: head
2289 real(DP),
intent(inout) :: vv
2296 topl = this%telev(iconn)
2297 botl = this%belev(iconn)
2298 ss = min(stage, topl)
2299 hh = min(head, topl)
2300 if (this%igwhcopt > 0)
then
2302 else if (this%inewton > 0)
then
2305 vv =
dhalf * (ss + hh)
2307 end subroutine lak_calculate_cond_head
2312 subroutine lak_calculate_conn_conductance(this, ilak, iconn, stage, head, cond)
2314 class(laktype),
intent(inout) :: this
2315 integer(I4B),
intent(in) :: ilak
2316 integer(I4B),
intent(in) :: iconn
2317 real(DP),
intent(in) :: stage
2318 real(DP),
intent(in) :: head
2319 real(DP),
intent(inout) :: cond
2321 integer(I4B) :: node
2329 real(DP) :: vscratio
2333 topl = this%telev(iconn)
2334 botl = this%belev(iconn)
2335 call this%lak_calculate_cond_head(iconn, stage, head, vv)
2340 if (this%ictype(iconn) == 0)
then
2341 if (abs(topl - botl) <
dprec)
then
2346 else if (this%ictype(iconn) == 1)
then
2347 node = this%cellid(iconn)
2348 if (this%icelltype(node) == 0)
then
2352 else if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
2353 node = this%cellid(iconn)
2354 if (this%icelltype(node) == 0)
then
2355 vv = this%telev(iconn)
2356 call this%lak_calculate_conn_warea(ilak, iconn, vv, vv, wa)
2358 call this%lak_calculate_conn_warea(ilak, iconn, stage, head, wa)
2364 if (this%ivsc == 1)
then
2366 if (stage > head)
then
2367 vscratio = this%viscratios(1, iconn)
2370 vscratio = this%viscratios(2, iconn)
2373 cond = sat * this%satcond(iconn) * vscratio
2374 end subroutine lak_calculate_conn_conductance
2378 subroutine lak_calculate_exchange(this, ilak, stage, totflow)
2380 class(laktype),
intent(inout) :: this
2381 integer(I4B),
intent(in) :: ilak
2382 real(DP),
intent(in) :: stage
2383 real(DP),
intent(inout) :: totflow
2386 integer(I4B) :: igwfnode
2391 do j = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2392 igwfnode = this%cellid(j)
2393 hgwf = this%xnew(igwfnode)
2394 call this%lak_calculate_conn_exchange(ilak, j, stage, hgwf, flow)
2395 totflow = totflow + flow
2397 end subroutine lak_calculate_exchange
2402 subroutine lak_calculate_conn_exchange(this, ilak, iconn, stage, head, flow, &
2405 class(laktype),
intent(inout) :: this
2406 integer(I4B),
intent(in) :: ilak
2407 integer(I4B),
intent(in) :: iconn
2408 real(DP),
intent(in) :: stage
2409 real(DP),
intent(in) :: head
2410 real(DP),
intent(inout) :: flow
2411 real(DP),
intent(inout),
optional :: gwfhcof
2412 real(DP),
intent(inout),
optional :: gwfrhs
2418 real(DP) :: gwfhcof0
2422 call this%lak_calculate_conn_conductance(ilak, iconn, stage, head, cond)
2423 botl = this%belev(iconn)
2426 if (stage >= botl)
then
2433 if (head >= botl)
then
2440 flow = cond * (hh - ss)
2443 if (head >= botl)
then
2445 gwfrhs0 = -cond * ss
2452 if (this%idense /= 0)
then
2453 call this%lak_calculate_density_exchange(iconn, stage, head, cond, botl, &
2454 flow, gwfhcof0, gwfrhs0)
2458 if (
present(gwfhcof)) gwfhcof = gwfhcof0
2459 if (
present(gwfrhs)) gwfrhs = gwfrhs0
2460 end subroutine lak_calculate_conn_exchange
2471 subroutine lak_calculate_conn_exchange_deriv(this, ilak, iconn, stage, &
2472 head, flow, dqds, dqdh)
2474 class(laktype),
intent(inout) :: this
2475 integer(I4B),
intent(in) :: ilak
2476 integer(I4B),
intent(in) :: iconn
2477 real(DP),
intent(in) :: stage
2478 real(DP),
intent(in) :: head
2479 real(DP),
intent(inout) :: flow
2480 real(DP),
intent(inout),
optional :: dqds
2481 real(DP),
intent(inout),
optional :: dqdh
2483 real(DP) :: cond, botl, dps, dph, ss, hh
2485 call this%lak_calculate_conn_conductance(ilak, iconn, stage, head, cond)
2486 botl = this%belev(iconn)
2487 if (stage >= botl)
then
2494 if (head >= botl)
then
2501 flow = cond * (hh - ss)
2502 if (
present(dqds)) dqds = -cond * dps
2503 if (
present(dqdh)) dqdh = cond * dph
2504 end subroutine lak_calculate_conn_exchange_deriv
2509 subroutine lak_estimate_conn_exchange(this, iflag, ilak, iconn, idry, stage, &
2510 head, flow, source, gwfhcof, gwfrhs)
2512 class(laktype),
intent(inout) :: this
2513 integer(I4B),
intent(in) :: iflag
2514 integer(I4B),
intent(in) :: ilak
2515 integer(I4B),
intent(in) :: iconn
2516 integer(I4B),
intent(inout) :: idry
2517 real(DP),
intent(in) :: stage
2518 real(DP),
intent(in) :: head
2519 real(DP),
intent(inout) :: flow
2520 real(DP),
intent(inout) :: source
2521 real(DP),
intent(inout),
optional :: gwfhcof
2522 real(DP),
intent(inout),
optional :: gwfrhs
2524 real(DP) :: gwfhcof0, gwfrhs0
2528 call this%lak_calculate_conn_exchange(ilak, iconn, stage, head, flow, &
2530 if (iflag == 1)
then
2531 if (flow >
dzero)
then
2532 source = source + flow
2534 else if (iflag == 2)
then
2535 if (-flow > source)
then
2539 else if (flow <
dzero)
then
2540 source = source + flow
2545 if (
present(gwfhcof)) gwfhcof = gwfhcof0
2546 if (
present(gwfrhs)) gwfrhs = gwfrhs0
2547 end subroutine lak_estimate_conn_exchange
2552 subroutine lak_calculate_storagechange(this, ilak, stage, stage0, delt, dvr)
2554 class(laktype),
intent(inout) :: this
2555 integer(I4B),
intent(in) :: ilak
2556 real(DP),
intent(in) :: stage
2557 real(DP),
intent(in) :: stage0
2558 real(DP),
intent(in) :: delt
2559 real(DP),
intent(inout) :: dvr
2565 if (this%gwfiss /= 1)
then
2566 call this%lak_calculate_vol(ilak, stage, v)
2567 call this%lak_calculate_vol(ilak, stage0, v0)
2568 dvr = (v0 - v) / delt
2570 end subroutine lak_calculate_storagechange
2574 subroutine lak_calculate_rainfall(this, ilak, stage, ra)
2576 class(laktype),
intent(inout) :: this
2577 integer(I4B),
intent(in) :: ilak
2578 real(DP),
intent(in) :: stage
2579 real(DP),
intent(inout) :: ra
2581 integer(I4B) :: iconn
2585 iconn = this%idxlakeconn(ilak)
2586 if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3)
then
2587 sa = this%sareamax(ilak)
2589 call this%lak_calculate_sarea(ilak, stage, sa)
2591 ra = this%rainfall(ilak) * sa
2592 end subroutine lak_calculate_rainfall
2596 subroutine lak_calculate_runoff(this, ilak, ro)
2598 class(laktype),
intent(inout) :: this
2599 integer(I4B),
intent(in) :: ilak
2600 real(DP),
intent(inout) :: ro
2603 ro = this%runoff(ilak)
2604 end subroutine lak_calculate_runoff
2608 subroutine lak_calculate_inflow(this, ilak, qin)
2610 class(laktype),
intent(inout) :: this
2611 integer(I4B),
intent(in) :: ilak
2612 real(DP),
intent(inout) :: qin
2615 qin = this%inflow(ilak)
2616 end subroutine lak_calculate_inflow
2620 subroutine lak_calculate_external(this, ilak, ex)
2622 class(laktype),
intent(inout) :: this
2623 integer(I4B),
intent(in) :: ilak
2624 real(DP),
intent(inout) :: ex
2629 if (this%imover == 1)
then
2630 ex = this%pakmvrobj%get_qfrommvr(ilak)
2632 end subroutine lak_calculate_external
2636 subroutine lak_calculate_withdrawal(this, ilak, avail, wr)
2638 class(laktype),
intent(inout) :: this
2639 integer(I4B),
intent(in) :: ilak
2640 real(DP),
intent(inout) :: avail
2641 real(DP),
intent(inout) :: wr
2644 wr = this%withdrawal(ilak)
2645 if (wr > avail)
then
2648 if (wr >
dzero)
then
2653 end subroutine lak_calculate_withdrawal
2658 subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
2660 class(laktype),
intent(inout) :: this
2661 integer(I4B),
intent(in) :: ilak
2662 real(DP),
intent(in) :: stage
2663 real(DP),
intent(inout) :: avail
2664 real(DP),
intent(inout) :: ev
2669 call this%lak_calculate_sarea(ilak, stage, sa)
2670 ev = sa * this%evaporation(ilak)
2671 if (ev > avail)
then
2681 end subroutine lak_calculate_evaporation
2685 subroutine lak_calculate_outlet_inflow(this, ilak, outinf)
2687 class(laktype),
intent(inout) :: this
2688 integer(I4B),
intent(in) :: ilak
2689 real(DP),
intent(inout) :: outinf
2694 do n = 1, this%noutlets
2695 if (this%lakeout(n) == ilak)
then
2696 outinf = outinf - this%simoutrate(n)
2697 if (this%imover == 1)
then
2698 outinf = outinf - this%pakmvrobj%get_qtomvr(n)
2702 end subroutine lak_calculate_outlet_inflow
2706 subroutine lak_calculate_outlet_outflow(this, ilak, stage, avail, outoutf)
2708 class(laktype),
intent(inout) :: this
2709 integer(I4B),
intent(in) :: ilak
2710 real(DP),
intent(in) :: stage
2711 real(DP),
intent(inout) :: avail
2712 real(DP),
intent(inout) :: outoutf
2722 do n = 1, this%noutlets
2723 if (this%lakein(n) == ilak)
then
2725 d = stage - this%outinvert(n)
2726 if (this%outdmax >
dzero)
then
2727 if (d > this%outdmax) d = this%outdmax
2729 g =
dgravity * this%convlength * this%convtime * this%convtime
2730 select case (this%iouttype(n))
2733 rate = this%outrate(n)
2734 if (-rate > avail)
then
2740 c = (this%convlength**
donethird) * this%convtime
2742 if (this%outrough(n) >
dzero)
then
2743 gsm =
done / this%outrough(n)
2745 rate = -c * gsm * this%outwidth(n) * (d**
dfivethirds) * &
2746 sqrt(this%outslope(n))
2755 this%simoutrate(n) = rate
2756 avail = avail + rate
2757 outoutf = outoutf + rate
2760 end subroutine lak_calculate_outlet_outflow
2770 subroutine lak_outlet_outflow_rate(this, ilak, stage, qout)
2772 class(laktype),
intent(inout) :: this
2773 integer(I4B),
intent(in) :: ilak
2774 real(DP),
intent(in) :: stage
2775 real(DP),
intent(inout) :: qout
2778 real(DP) :: g, d, c, gsm, rate
2781 do n = 1, this%noutlets
2782 if (this%lakein(n) /= ilak) cycle
2784 d = stage - this%outinvert(n)
2785 if (this%outdmax >
dzero .and. d > this%outdmax) d = this%outdmax
2786 g =
dgravity * this%convlength * this%convtime * this%convtime
2787 select case (this%iouttype(n))
2789 rate = this%outrate(n)
2792 c = (this%convlength**
donethird) * this%convtime
2794 if (this%outrough(n) >
dzero) gsm =
done / this%outrough(n)
2795 rate = -c * gsm * this%outwidth(n) * (d**
dfivethirds) * &
2796 sqrt(this%outslope(n))
2805 end subroutine lak_outlet_outflow_rate
2809 subroutine lak_get_internal_inlet(this, ilak, outinf)
2811 class(laktype),
intent(inout) :: this
2812 integer(I4B),
intent(in) :: ilak
2813 real(DP),
intent(inout) :: outinf
2818 do n = 1, this%noutlets
2819 if (this%lakeout(n) == ilak)
then
2820 outinf = outinf - this%simoutrate(n)
2821 if (this%imover == 1)
then
2822 outinf = outinf - this%pakmvrobj%get_qtomvr(n)
2826 end subroutine lak_get_internal_inlet
2830 subroutine lak_get_internal_outlet(this, ilak, outoutf)
2832 class(laktype),
intent(inout) :: this
2833 integer(I4B),
intent(in) :: ilak
2834 real(DP),
intent(inout) :: outoutf
2839 do n = 1, this%noutlets
2840 if (this%lakein(n) == ilak)
then
2841 if (this%lakeout(n) < 1) cycle
2842 outoutf = outoutf + this%simoutrate(n)
2845 end subroutine lak_get_internal_outlet
2849 subroutine lak_get_external_outlet(this, ilak, outoutf)
2851 class(laktype),
intent(inout) :: this
2852 integer(I4B),
intent(in) :: ilak
2853 real(DP),
intent(inout) :: outoutf
2858 do n = 1, this%noutlets
2859 if (this%lakein(n) == ilak)
then
2860 if (this%lakeout(n) > 0) cycle
2861 outoutf = outoutf + this%simoutrate(n)
2864 end subroutine lak_get_external_outlet
2868 subroutine lak_get_external_mover(this, ilak, outoutf)
2870 class(laktype),
intent(inout) :: this
2871 integer(I4B),
intent(in) :: ilak
2872 real(DP),
intent(inout) :: outoutf
2877 if (this%imover == 1)
then
2878 do n = 1, this%noutlets
2879 if (this%lakein(n) == ilak)
then
2880 if (this%lakeout(n) > 0) cycle
2881 outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2885 end subroutine lak_get_external_mover
2889 subroutine lak_get_internal_mover(this, ilak, outoutf)
2891 class(laktype),
intent(inout) :: this
2892 integer(I4B),
intent(in) :: ilak
2893 real(DP),
intent(inout) :: outoutf
2898 if (this%imover == 1)
then
2899 do n = 1, this%noutlets
2900 if (this%lakein(n) == ilak)
then
2901 if (this%lakeout(n) < 1) cycle
2902 outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2906 end subroutine lak_get_internal_mover
2910 subroutine lak_get_outlet_tomover(this, ilak, outoutf)
2912 class(laktype),
intent(inout) :: this
2913 integer(I4B),
intent(in) :: ilak
2914 real(DP),
intent(inout) :: outoutf
2919 if (this%imover == 1)
then
2920 do n = 1, this%noutlets
2921 if (this%lakein(n) == ilak)
then
2922 outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2926 end subroutine lak_get_outlet_tomover
2930 subroutine lak_vol2stage(this, ilak, vol, stage)
2932 class(laktype),
intent(inout) :: this
2933 integer(I4B),
intent(in) :: ilak
2934 real(DP),
intent(in) :: vol
2935 real(DP),
intent(inout) :: stage
2939 real(DP) :: s0, s1, sm
2940 real(DP) :: v0, v1, vm
2941 real(DP) :: f0, f1, fm
2943 real(DP) :: en0, en1
2947 s0 = this%lakebot(ilak)
2948 call this%lak_calculate_vol(ilak, s0, v0)
2949 s1 = this%laketop(ilak)
2950 call this%lak_calculate_vol(ilak, s1, v1)
2955 else if (vol >= v1)
then
2956 call this%lak_calculate_sarea(ilak, s1, sa)
2957 stage = s1 + (vol - v1) / sa
2968 secantbisection:
do i = 1, 150
2970 if (denom /=
dzero)
then
2971 ds = f1 * (s1 - s0) / denom
2979 if (sm < en0 .or. sm > en1) ibs = 13
2983 if (ds * ds0 <
dprec .or. abs(ds) > abs(ds0)) ibs = ibs + 1
2985 ds =
dhalf * (s1 - s0)
2989 if (abs(ds) <
dem6)
then
2990 exit secantbisection
2992 call this%lak_calculate_vol(ilak, sm, vm)
2999 end do secantbisection
3001 if (abs(ds) >=
dem6)
then
3002 write (this%iout,
'(1x,a,1x,i0,4(1x,a,1x,g15.6))') &
3003 &
'LAK_VOL2STAGE failed for lake', ilak,
'volume error =', fm, &
3004 &
'finding stage (', stage,
') for volume =', vol, &
3005 &
'final change in stage =', ds
3008 end subroutine lak_vol2stage
3011 function lak_check_valid(this, itemno)
result(ierr)
3015 integer(I4B) :: ierr
3017 class(laktype),
intent(inout) :: this
3018 integer(I4B),
intent(in) :: itemno
3020 integer(I4B) :: ival
3024 if (itemno > 0)
then
3025 if (ival < 1 .or. ival > this%nlakes)
then
3026 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,a)') &
3027 'LAKENO', itemno,
'must be greater than 0 and less than or equal to', &
3033 if (ival < 1 .or. ival > this%noutlets)
then
3034 write (
errmsg,
'(a,1x,i0,1x,a,1x,i0,a)') &
3035 'IOUTLET', itemno,
'must be greater than 0 and less than or equal to', &
3041 end function lak_check_valid
3045 subroutine lak_set_stressperiod(this, itemno)
3050 class(laktype),
intent(inout) :: this
3051 integer(I4B),
intent(in) :: itemno
3053 character(len=LINELENGTH) :: text
3054 character(len=LINELENGTH) :: caux
3055 character(len=LINELENGTH) :: keyword
3056 integer(I4B) :: ierr
3059 real(DP),
pointer :: bndElem => null()
3062 call this%parser%GetStringCaps(keyword)
3063 select case (keyword)
3065 ierr = this%lak_check_valid(itemno)
3069 call this%parser%GetStringCaps(text)
3070 this%status(itemno) = text(1:8)
3071 if (text ==
'CONSTANT')
then
3072 this%iboundpak(itemno) = -1
3073 else if (text ==
'INACTIVE')
then
3074 this%iboundpak(itemno) = 0
3075 else if (text ==
'ACTIVE')
then
3076 this%iboundpak(itemno) = 1
3078 write (
errmsg,
'(a,a)') &
3079 'Unknown '//trim(this%text)//
' lak status keyword: ', text//
'.'
3083 ierr = this%lak_check_valid(itemno)
3087 call this%parser%GetString(text)
3089 bndelem => this%stage(itemno)
3091 this%packName,
'BND', this%tsManager, &
3092 this%iprpak,
'STAGE')
3094 ierr = this%lak_check_valid(itemno)
3098 call this%parser%GetString(text)
3100 bndelem => this%rainfall(itemno)
3102 this%packName,
'BND', this%tsManager, &
3103 this%iprpak,
'RAINFALL')
3104 if (this%rainfall(itemno) <
dzero)
then
3105 write (
errmsg,
'(a,i0,a,G0,a)') &
3106 'Lake ', itemno,
' was assigned a rainfall value of ', &
3107 this%rainfall(itemno),
'. Rainfall must be positive.'
3110 case (
'EVAPORATION')
3111 ierr = this%lak_check_valid(itemno)
3115 call this%parser%GetString(text)
3117 bndelem => this%evaporation(itemno)
3119 this%packName,
'BND', this%tsManager, &
3120 this%iprpak,
'EVAPORATION')
3121 if (this%evaporation(itemno) <
dzero)
then
3122 write (
errmsg,
'(a,i0,a,G0,a)') &
3123 'Lake ', itemno,
' was assigned an evaporation value of ', &
3124 this%evaporation(itemno),
'. Evaporation must be positive.'
3128 ierr = this%lak_check_valid(itemno)
3132 call this%parser%GetString(text)
3134 bndelem => this%runoff(itemno)
3136 this%packName,
'BND', this%tsManager, &
3137 this%iprpak,
'RUNOFF')
3138 if (this%runoff(itemno) <
dzero)
then
3139 write (
errmsg,
'(a,i0,a,G0,a)') &
3140 'Lake ', itemno,
' was assigned a runoff value of ', &
3141 this%runoff(itemno),
'. Runoff must be positive.'
3145 ierr = this%lak_check_valid(itemno)
3149 call this%parser%GetString(text)
3151 bndelem => this%inflow(itemno)
3153 this%packName,
'BND', this%tsManager, &
3154 this%iprpak,
'INFLOW')
3155 if (this%inflow(itemno) <
dzero)
then
3156 write (
errmsg,
'(a,i0,a,G0,a)') &
3157 'Lake ', itemno,
' was assigned an inflow value of ', &
3158 this%inflow(itemno),
'. Inflow must be positive.'
3162 ierr = this%lak_check_valid(itemno)
3166 call this%parser%GetString(text)
3168 bndelem => this%withdrawal(itemno)
3170 this%packName,
'BND', this%tsManager, &
3171 this%iprpak,
'WITHDRAWAL')
3172 if (this%withdrawal(itemno) <
dzero)
then
3173 write (
errmsg,
'(a,i0,a,G0,a)') &
3174 'Lake ', itemno,
' was assigned a withdrawal value of ', &
3175 this%withdrawal(itemno),
'. Withdrawal must be positive.'
3179 ierr = this%lak_check_valid(-itemno)
3183 call this%parser%GetString(text)
3185 bndelem => this%outrate(itemno)
3187 this%packName,
'BND', this%tsManager, &
3188 this%iprpak,
'RATE')
3190 ierr = this%lak_check_valid(-itemno)
3194 call this%parser%GetString(text)
3196 bndelem => this%outinvert(itemno)
3198 this%packName,
'BND', this%tsManager, &
3199 this%iprpak,
'INVERT')
3201 ierr = this%lak_check_valid(-itemno)
3205 call this%parser%GetString(text)
3207 bndelem => this%outwidth(itemno)
3209 this%packName,
'BND', this%tsManager, &
3210 this%iprpak,
'WIDTH')
3212 ierr = this%lak_check_valid(-itemno)
3216 call this%parser%GetString(text)
3218 bndelem => this%outrough(itemno)
3220 this%packName,
'BND', this%tsManager, &
3221 this%iprpak,
'ROUGH')
3223 ierr = this%lak_check_valid(-itemno)
3227 call this%parser%GetString(text)
3229 bndelem => this%outslope(itemno)
3231 this%packName,
'BND', this%tsManager, &
3232 this%iprpak,
'SLOPE')
3234 ierr = this%lak_check_valid(itemno)
3238 call this%parser%GetStringCaps(caux)
3239 do jj = 1, this%naux
3240 if (trim(adjustl(caux)) /= trim(adjustl(this%auxname(jj)))) cycle
3241 call this%parser%GetString(text)
3243 bndelem => this%lauxvar(jj, ii)
3245 this%packName,
'AUX', &
3246 this%tsManager, this%iprpak, &
3252 'Unknown '//trim(this%text)//
' lak data keyword: ', &
3258 end subroutine lak_set_stressperiod
3264 subroutine lak_set_attribute_error(this, ilak, keyword, msg)
3268 class(laktype),
intent(inout) :: this
3269 integer(I4B),
intent(in) :: ilak
3270 character(len=*),
intent(in) :: keyword
3271 character(len=*),
intent(in) :: msg
3273 if (len(msg) == 0)
then
3274 write (
errmsg,
'(a,1x,a,1x,i0,1x,a)') &
3275 keyword,
' for LAKE', ilak,
'has already been set.'
3277 write (
errmsg,
'(a,1x,a,1x,i0,1x,a)') keyword,
' for LAKE', ilak, msg
3280 end subroutine lak_set_attribute_error
3286 subroutine lak_options(this, option, found)
3293 class(laktype),
intent(inout) :: this
3294 character(len=*),
intent(inout) :: option
3295 logical(LGP),
intent(inout) :: found
3297 character(len=MAXCHARLEN) :: fname, keyword
3300 character(len=*),
parameter :: fmtlengthconv = &
3301 &
"(4x, 'LENGTH CONVERSION VALUE (',g15.7,') SPECIFIED.')"
3302 character(len=*),
parameter :: fmttimeconv = &
3303 &
"(4x, 'TIME CONVERSION VALUE (',g15.7,') SPECIFIED.')"
3304 character(len=*),
parameter :: fmtoutdmax = &
3305 &
"(4x, 'MAXIMUM OUTLET WATER DEPTH (',g15.7,') SPECIFIED.')"
3306 character(len=*),
parameter :: fmtlakeopt = &
3307 &
"(4x, 'LAKE ', a, ' VALUE (',g15.7,') SPECIFIED.')"
3308 character(len=*),
parameter :: fmtlakbin = &
3309 "(4x, 'LAK ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', &
3310 &a, /4x, 'OPENED ON UNIT: ', I0)"
3311 character(len=*),
parameter :: fmtiter = &
3312 &
"(4x, 'MAXIMUM LAK ITERATION VALUE (',i0,') SPECIFIED.')"
3313 character(len=*),
parameter :: fmtdmaxchg = &
3314 &
"(4x, 'MAXIMUM STAGE CHANGE VALUE (',g0,') SPECIFIED.')"
3317 select case (option)
3318 case (
'PRINT_STAGE')
3320 write (this%iout,
'(4x,a)') trim(adjustl(this%text))// &
3321 ' STAGES WILL BE PRINTED TO LISTING FILE.'
3323 call this%parser%GetStringCaps(keyword)
3324 if (keyword ==
'FILEOUT')
then
3325 call this%parser%GetString(fname)
3327 call openfile(this%istageout, this%iout, fname,
'DATA(BINARY)', &
3329 write (this%iout, fmtlakbin)
'STAGE', trim(adjustl(fname)), &
3332 call store_error(
'OPTIONAL STAGE KEYWORD MUST BE FOLLOWED BY FILEOUT')
3335 call this%parser%GetStringCaps(keyword)
3336 if (keyword ==
'FILEOUT')
then
3337 call this%parser%GetString(fname)
3338 call assign_iounit(this%ibudgetout, this%inunit,
"BUDGET fileout")
3339 call openfile(this%ibudgetout, this%iout, fname,
'DATA(BINARY)', &
3341 write (this%iout, fmtlakbin)
'BUDGET', trim(adjustl(fname)), &
3344 call store_error(
'OPTIONAL BUDGET 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%ibudcsv, this%inunit,
"BUDGETCSV fileout")
3351 call openfile(this%ibudcsv, this%iout, fname,
'CSV', &
3352 filstat_opt=
'REPLACE')
3353 write (this%iout, fmtlakbin)
'BUDGET CSV', trim(adjustl(fname)), &
3356 call store_error(
'OPTIONAL BUDGETCSV KEYWORD MUST BE FOLLOWED BY &
3359 case (
'PACKAGE_CONVERGENCE')
3360 call this%parser%GetStringCaps(keyword)
3361 if (keyword ==
'FILEOUT')
then
3362 call this%parser%GetString(fname)
3366 this%pakcsvfile = trim(adjustl(fname))
3368 call store_error(
'OPTIONAL PACKAGE_CONVERGENCE KEYWORD MUST BE '// &
3369 'FOLLOWED BY FILEOUT')
3373 write (this%iout,
'(4x,A)')
'MOVER OPTION ENABLED'
3374 case (
'LENGTH_CONVERSION')
3375 this%convlength = this%parser%GetDouble()
3376 write (this%iout, fmtlengthconv) this%convlength
3377 case (
'TIME_CONVERSION')
3378 this%convtime = this%parser%GetDouble()
3379 write (this%iout, fmttimeconv) this%convtime
3381 r = this%parser%GetDouble()
3386 write (this%iout, fmtlakeopt)
'SURFDEP', this%surfdep
3387 case (
'MAXIMUM_ITERATIONS')
3388 this%maxlakit = this%parser%GetInteger()
3389 write (this%iout, fmtiter) this%maxlakit
3390 case (
'MAXIMUM_STAGE_CHANGE')
3391 r = this%parser%GetDouble()
3393 this%delh = dp999 * r
3394 write (this%iout, fmtdmaxchg) this%dmaxchg
3400 case (
'DEV_GROUNDWATER_HEAD_CONDUCTANCE')
3401 call this%parser%DevOpt()
3403 write (this%iout,
'(4x,a)') &
3404 'CONDUCTANCE FOR HORIZONTAL CONNECTIONS WILL BE CALCULATED &
3405 &USING THE GROUNDWATER HEAD'
3406 case (
'DEV_MAXIMUM_OUTLET_DEPTH')
3407 call this%parser%DevOpt()
3408 this%outdmax = this%parser%GetDouble()
3409 write (this%iout, fmtoutdmax) this%outdmax
3412 write (this%iout,
'(4x,a)') &
3413 'LAKE STAGE WILL BE SOLVED AS AN UNKNOWN IN THE GROUNDWATER FLOW '// &
3414 'MATRIX (IMPLICIT FORMULATION)'
3415 case (
'DEV_FORCE_FALLBACK')
3416 call this%parser%DevOpt()
3418 write (this%iout,
'(4x,a)') &
3419 'EVERY ACTIVE LAKE WILL BE SOLVED WITH THE SUBSTITUTION FALLBACK '// &
3420 'UNDER THE IMPLICIT FORMULATION'
3421 case (
'DEV_NO_FINAL_CHECK')
3422 call this%parser%DevOpt()
3424 write (this%iout,
'(4x,a)') &
3425 'A FINAL CONVERGENCE CHECK OF THE CHANGE IN LAKE STAGES &
3432 end subroutine lak_options
3438 subroutine lak_ar(this)
3443 class(laktype),
intent(inout) :: this
3445 character(len=*),
parameter :: fmtlakbin = &
3446 "(4x, 'LAK ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', &
3447 &a, /4x, 'OPENED ON UNIT: ', I0)"
3454 if (
allocated(this%pakcsvfile))
then
3455 if (this%iimplicit /= 0)
then
3457 'PACKAGE_CONVERGENCE output file "'//trim(this%pakcsvfile)// &
3458 '" is not written when the IMPLICIT option is active; the lake '// &
3459 'stage is part of the solver (IMS) convergence check.'
3463 call openfile(this%ipakcsv, this%iout, this%pakcsvfile,
'CSV', &
3464 filstat_opt=
'REPLACE', mode_opt=
mnormal)
3465 write (this%iout, fmtlakbin)
'PACKAGE_CONVERGENCE', &
3466 trim(this%pakcsvfile), this%ipakcsv
3468 deallocate (this%pakcsvfile)
3471 call this%obs%obs_ar()
3474 call this%lak_allocate_arrays()
3477 call this%read_initial_attr()
3480 if (this%imover /= 0)
then
3481 allocate (this%pakmvrobj)
3482 call this%pakmvrobj%ar(this%noutlets, this%nlakes, this%memoryPath)
3484 end subroutine lak_ar
3490 subroutine lak_rp(this)
3496 class(laktype),
intent(inout) :: this
3498 character(len=LINELENGTH) :: title
3499 character(len=LINELENGTH) :: line
3500 character(len=LINELENGTH) :: text
3501 logical(LGP) :: isfound
3502 logical(LGP) :: endOfBlock
3503 integer(I4B) :: ierr
3504 integer(I4B) :: node
3506 integer(I4B) :: itemno
3509 character(len=*),
parameter :: fmtblkerr = &
3510 &
"('Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
3511 character(len=*),
parameter :: fmtlsp = &
3512 &
"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
3515 this%nbound = this%maxbound
3519 if (this%inunit == 0)
return
3522 if (this%ionper <
kper)
then
3525 call this%parser%GetBlock(
'PERIOD', isfound, ierr, &
3526 supportopenclose=.true., &
3527 blockrequired=.false.)
3531 call this%read_check_ionper()
3537 this%ionper =
nper + 1
3540 call this%parser%GetCurrentLine(line)
3541 write (
errmsg, fmtblkerr) adjustl(trim(line))
3543 call this%parser%StoreErrorUnit()
3549 if (this%ionper ==
kper)
then
3552 if (this%iprpak /= 0)
then
3555 title = trim(adjustl(this%text))//
' PACKAGE ('// &
3556 trim(adjustl(this%packName))//
') DATA FOR PERIOD'
3557 write (title,
'(a,1x,i6)') trim(adjustl(title)),
kper
3558 call table_cr(this%inputtab, this%packName, title)
3559 call this%inputtab%table_df(1, 4, this%iout, finalize=.false.)
3561 call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
3563 call this%inputtab%initialize_column(text, 20, alignment=tableft)
3565 write (text,
'(a,1x,i6)')
'VALUE', n
3566 call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
3573 call this%parser%GetNextLine(endofblock)
3574 if (endofblock)
exit
3577 itemno = this%parser%GetInteger()
3580 call this%lak_set_stressperiod(itemno)
3583 if (this%iprpak /= 0)
then
3584 call this%parser%GetCurrentLine(line)
3585 call this%inputtab%line_to_columns(line)
3589 if (this%iprpak /= 0)
then
3590 call this%inputtab%finalize_table()
3595 write (this%iout, fmtlsp) trim(this%filtyp)
3600 call this%parser%StoreErrorUnit()
3604 do n = 1, this%nlakes
3605 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3606 node = this%cellid(j)
3607 this%nodelist(j) = node
3608 this%bound(1, j) = this%xnewpak(n)
3609 this%bound(2, j) = this%satcond(j)
3610 this%bound(3, j) = this%belev(j)
3615 if (this%imover == 1)
then
3616 do n = 1, this%noutlets
3617 this%pakmvrobj%iprmap(n) = this%lakein(n)
3620 end subroutine lak_rp
3624 subroutine lak_ad(this)
3628 class(laktype) :: this
3632 integer(I4B) :: iaux
3635 call this%TsManager%ad()
3640 if (this%naux > 0)
then
3641 do n = 1, this%nlakes
3642 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3643 do iaux = 1, this%naux
3644 if (this%noupdateauxvar(iaux) /= 0) cycle
3645 this%auxvar(iaux, j) = this%lauxvar(iaux, n)
3656 do n = 1, this%nlakes
3657 this%xoldpak(n) = this%xnewpak(n)
3658 this%stageiter(n) = this%xnewpak(n)
3659 if (this%iboundpak(n) < 0)
then
3660 this%xnewpak(n) = this%stage(n)
3662 this%seep0(n) =
dzero
3668 do n = 1, this%nlakes
3669 this%xnewpak(n) = this%xoldpak(n)
3670 this%stageiter(n) = this%xnewpak(n)
3671 if (this%iboundpak(n) < 0)
then
3672 this%xnewpak(n) = this%stage(n)
3674 this%seep0(n) =
dzero
3679 if (this%imover == 1)
then
3680 call this%pakmvrobj%ad()
3686 call this%obs%obs_ad()
3687 end subroutine lak_ad
3693 subroutine lak_cf(this)
3695 class(laktype) :: this
3697 integer(I4B) :: j, n
3698 integer(I4B) :: igwfnode
3699 real(DP) :: hlak, bottom_lake
3702 do n = 1, this%nlakes
3703 this%seep0(n) = this%seep(n)
3707 do n = 1, this%nlakes
3708 this%s0(n) = this%xnewpak(n)
3709 call this%lak_calculate_exchange(n, this%s0(n), this%qgwf0(n))
3713 do n = 1, this%nlakes
3714 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3716 if (this%ictype(j) /= 0)
then
3719 igwfnode = this%nodesontop(j)
3720 if (this%ibound(igwfnode) == 0)
then
3721 call this%dis%highest_active(igwfnode, this%ibound)
3723 this%nodelist(j) = igwfnode
3724 this%cellid(j) = igwfnode
3731 do n = 1, this%nlakes
3733 hlak = this%xnewpak(n)
3736 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3739 igwfnode = this%cellid(j)
3742 if (this%ibound(igwfnode) < 1)
then
3747 if (this%ictype(j) /= 0)
then
3752 if (this%ictype(j) == 2 .or. this%ictype(j) == 3)
then
3757 bottom_lake = this%belev(j)
3758 if (hlak > bottom_lake .or. this%iboundpak(n) == 0)
then
3761 this%ibound(igwfnode) = 1
3769 call this%lak_bound_update()
3770 end subroutine lak_cf
3774 subroutine lak_fc(this, rhs, ia, idxglo, matrix_sln)
3776 class(laktype) :: this
3777 real(DP),
dimension(:),
intent(inout) :: rhs
3778 integer(I4B),
dimension(:),
intent(in) :: ia
3779 integer(I4B),
dimension(:),
intent(in) :: idxglo
3782 integer(I4B) :: j, n
3783 integer(I4B) :: igwfnode
3784 integer(I4B) :: ipossymd
3787 if (this%imover == 1)
then
3788 call this%pakmvrobj%fc()
3793 if (this%iimplicit /= 0)
then
3794 call this%lak_fc_implicit(rhs, matrix_sln)
3800 call this%lak_solve()
3801 do n = 1, this%nlakes
3802 if (this%iboundpak(n) == 0) cycle
3803 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3804 igwfnode = this%cellid(j)
3805 if (this%ibound(igwfnode) < 1) cycle
3806 ipossymd = idxglo(ia(igwfnode))
3807 call matrix_sln%add_value_pos(ipossymd, this%hcof(j))
3808 rhs(igwfnode) = rhs(igwfnode) + this%rhs(j)
3811 end subroutine lak_fc
3815 subroutine lak_fn(this, rhs, ia, idxglo, matrix_sln)
3817 class(laktype) :: this
3818 real(DP),
dimension(:),
intent(inout) :: rhs
3819 integer(I4B),
dimension(:),
intent(in) :: ia
3820 integer(I4B),
dimension(:),
intent(in) :: idxglo
3823 integer(I4B) :: j, n
3824 integer(I4B) :: ipos
3825 integer(I4B) :: igwfnode
3826 integer(I4B) :: idry
3843 if (this%iimplicit /= 0)
then
3847 do n = 1, this%nlakes
3848 if (this%iboundpak(n) == 0) cycle
3849 hlak = this%xnewpak(n)
3850 call this%lak_calculate_available(n, hlak, avail, &
3851 ra, ro, qinf, ex, this%delh)
3852 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3853 igwfnode = this%cellid(j)
3855 head = this%xnew(igwfnode)
3856 if (-this%hcof(j) >
dzero)
then
3857 if (this%ibound(igwfnode) > 0)
then
3861 call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, &
3862 head + this%delh, q1, avail)
3865 q = this%hcof(j) * head - this%rhs(j)
3867 rterm = this%hcof(j) * head
3869 drterm = (q1 - q) / this%delh
3872 call matrix_sln%add_value_pos(idxglo(ipos), drterm - this%hcof(j))
3873 rhs(igwfnode) = rhs(igwfnode) - rterm + drterm * head
3878 end subroutine lak_fn
3888 subroutine lak_nur(this, neqpak, x, xtemp, dx, inewtonur, dxmax, locmax)
3890 class(laktype),
intent(inout) :: this
3891 integer(I4B),
intent(in) :: neqpak
3892 real(DP),
dimension(neqpak),
intent(inout) :: x
3893 real(DP),
dimension(neqpak),
intent(in) :: xtemp
3894 real(DP),
dimension(neqpak),
intent(inout) :: dx
3895 integer(I4B),
intent(inout) :: inewtonur
3896 real(DP),
intent(inout) :: dxmax
3897 integer(I4B),
intent(inout) :: locmax
3905 if (this%iimplicit == 0)
return
3908 do n = 1, this%nlakes
3909 if (this%iboundpak(n) < 1) cycle
3910 botl = this%lakebot(n)
3914 if (x(n) < botl)
then
3918 if (abs(dxx) > abs(dxmax))
then
3926 end subroutine lak_nur
3930 subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
3934 class(laktype),
intent(inout) :: this
3935 integer(I4B),
intent(in) :: innertot
3936 integer(I4B),
intent(in) :: kiter
3937 integer(I4B),
intent(in) :: iend
3938 integer(I4B),
intent(in) :: icnvgmod
3939 character(len=LENPAKLOC),
intent(inout) :: cpak
3940 integer(I4B),
intent(inout) :: ipak
3941 real(DP),
intent(inout) :: dpak
3943 character(len=LENPAKLOC) :: cloc
3944 character(len=LINELENGTH) :: tag
3945 integer(I4B) :: icheck
3946 integer(I4B) :: ipakfail
3947 integer(I4B) :: locdhmax
3948 integer(I4B) :: locresidmax
3949 integer(I4B) :: locdgwfmax
3950 integer(I4B) :: locdqoutmax
3951 integer(I4B) :: locdqfrommvrmax
3952 integer(I4B) :: ntabrows
3953 integer(I4B) :: ntabcols
3957 real(DP) :: qtolfact
3975 real(DP) :: residmax
3977 real(DP) :: dqoutmax
3978 real(DP) :: dqfrommvr
3979 real(DP) :: dqfrommvrmax
3981 call this%lak_set_fallback(kiter, icnvgmod)
3986 if (iend /= 0 .and. icnvgmod == 0)
then
3987 call this%lak_check_disconnected()
3993 if (this%iimplicit /= 0)
then
3998 icheck = this%iconvchk
4009 dqfrommvrmax =
dzero
4013 if (this%ipakcsv == 0)
then
4014 if (icnvgmod == 0)
then
4022 if (.not.
associated(this%pakcsvtab))
then
4027 if (this%noutlets > 0)
then
4028 ntabcols = ntabcols + 2
4030 if (this%imover == 1)
then
4031 ntabcols = ntabcols + 2
4035 call table_cr(this%pakcsvtab, this%packName,
'')
4036 call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
4037 lineseparator=.false., separator=
',', &
4041 tag =
'total_inner_iterations'
4042 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4044 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4046 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4048 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4050 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
4052 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4054 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4056 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4057 tag =
'residmax_loc'
4058 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4060 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4062 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4063 if (this%noutlets > 0)
then
4065 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4066 tag =
'dqoutmax_loc'
4067 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4069 if (this%imover == 1)
then
4070 tag =
'dqfrommvrmax'
4071 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
4072 tag =
'dqfrommvrmax_loc'
4073 call this%pakcsvtab%initialize_column(tag, 16, alignment=
tableft)
4079 if (icheck /= 0)
then
4080 final_check:
do n = 1, this%nlakes
4081 if (this%iboundpak(n) < 1) cycle
4085 hlak = this%xnewpak(n)
4091 call this%lak_calculate_sarea(n, hlak, area)
4094 if (area >
dzero)
then
4095 qtolfact =
delt / area
4101 call this%lak_calculate_residual(n, hlak, resid)
4102 resid = resid * qtolfact
4106 if (area >
dzero)
then
4107 gwf0 = this%qgwf0(n)
4108 call this%lak_calculate_exchange(n, hlak, gwf)
4109 dgwf = (gwf0 - gwf) * qtolfact
4114 if (this%noutlets > 0)
then
4115 if (area >
dzero)
then
4116 call this%lak_calculate_available(n, hlak0, inf, ra, ro, qinf, ex)
4117 call this%lak_calculate_outlet_outflow(n, hlak0, inf, qout0)
4118 call this%lak_calculate_available(n, hlak, inf, ra, ro, qinf, ex)
4119 call this%lak_calculate_outlet_outflow(n, hlak, inf, qout)
4120 dqout = (qout0 - qout) * qtolfact
4126 if (this%imover == 1)
then
4127 q = this%pakmvrobj%get_qfrommvr(n)
4128 q0 = this%pakmvrobj%get_qfrommvr0(n)
4129 dqfrommvr = qtolfact * (q0 - q)
4142 dqfrommvrmax = dqfrommvr
4145 if (abs(dh) > abs(dhmax))
then
4149 if (abs(resid) > abs(residmax))
then
4153 if (abs(dgwf) > abs(dgwfmax))
then
4157 if (abs(dqout) > abs(dqoutmax))
then
4161 if (abs(dqfrommvr) > abs(dqfrommvrmax))
then
4162 dqfrommvrmax = dqfrommvr
4169 if (abs(dhmax) > abs(dpak))
then
4172 write (cloc,
"(a,'-',a)") &
4173 trim(this%packName),
'stage'
4176 if (abs(residmax) > abs(dpak))
then
4179 write (cloc,
"(a,'-',a)") &
4180 trim(this%packName),
'residual'
4183 if (abs(dgwfmax) > abs(dpak))
then
4186 write (cloc,
"(a,'-',a)") &
4187 trim(this%packName),
'gwf'
4190 if (this%noutlets > 0)
then
4191 if (abs(dqoutmax) > abs(dpak))
then
4194 write (cloc,
"(a,'-',a)") &
4195 trim(this%packName),
'outlet'
4199 if (this%imover == 1)
then
4200 if (abs(dqfrommvrmax) > abs(dpak))
then
4201 ipak = locdqfrommvrmax
4203 write (cloc,
"(a,'-',a)") trim(this%packName),
'qfrommvr'
4209 if (this%ipakcsv /= 0)
then
4212 call this%pakcsvtab%add_term(innertot)
4213 call this%pakcsvtab%add_term(
totim)
4214 call this%pakcsvtab%add_term(
kper)
4215 call this%pakcsvtab%add_term(
kstp)
4216 call this%pakcsvtab%add_term(kiter)
4217 call this%pakcsvtab%add_term(dhmax)
4218 call this%pakcsvtab%add_term(locdhmax)
4219 call this%pakcsvtab%add_term(residmax)
4220 call this%pakcsvtab%add_term(locresidmax)
4221 call this%pakcsvtab%add_term(dgwfmax)
4222 call this%pakcsvtab%add_term(locdgwfmax)
4223 if (this%noutlets > 0)
then
4224 call this%pakcsvtab%add_term(dqoutmax)
4225 call this%pakcsvtab%add_term(locdqoutmax)
4227 if (this%imover == 1)
then
4228 call this%pakcsvtab%add_term(dqfrommvrmax)
4229 call this%pakcsvtab%add_term(locdqfrommvrmax)
4234 call this%pakcsvtab%finalize_table()
4238 end subroutine lak_cc
4242 subroutine lak_cq(this, x, flowja, iadv)
4246 class(laktype),
intent(inout) :: this
4247 real(DP),
dimension(:),
intent(in) :: x
4248 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
4249 integer(I4B),
optional,
intent(in) :: iadv
4252 real(DP) :: chratin, chratout
4254 integer(I4B) :: j, n, igwfnode
4255 real(DP) :: hlak, head, flow, dqdh
4256 real(DP) :: v0, v1, sa, sf
4258 call this%lak_solve(update=.false.)
4267 if (this%iimplicit /= 0)
then
4268 do n = 1, this%nlakes
4269 if (this%iboundpak(n) < 1 .or. this%ifallback(n) /= 0) cycle
4270 hlak = this%xnewpak(n)
4273 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4274 igwfnode = this%cellid(j)
4275 if (this%ibound(igwfnode) < 1) cycle
4276 head = this%xnew(igwfnode)
4277 call this%lak_calculate_conn_exchange_deriv(n, j, hlak, head, &
4279 this%hcof(j) = -dqdh
4280 this%rhs(j) = -dqdh * head + flow
4290 call this%lak_calculate_sarea(n, hlak, sa)
4292 if (this%surfdep >
dzero)
then
4294 this%lakebot(n), hlak)
4296 this%evap(n) = -this%evaporation(n) * sa * sf
4297 this%withr(n) = -this%withdrawal(n) * sf
4303 call this%BndType%bnd_cq(x, flowja, iadv=1)
4308 do n = 1, this%nlakes
4309 this%chterm(n) =
dzero
4310 if (this%iboundpak(n) == 0) cycle
4311 hlak = this%xnewpak(n)
4312 call this%lak_calculate_vol(n, hlak, v1)
4315 if (this%iboundpak(n) /= 0)
then
4318 rrate = this%precip(n)
4319 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4322 rrate = this%evap(n)
4323 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4326 rrate = this%runoff(n)
4327 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4330 rrate = this%inflow(n)
4331 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4334 rrate = this%withr(n)
4335 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4339 if (this%iboundpak(n) > 0)
then
4340 if (this%gwfiss /= 1)
then
4341 call this%lak_calculate_vol(n, this%xoldpak(n), v0)
4342 rrate = -(v1 - v0) /
delt
4343 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4346 this%qsto(n) = rrate
4349 call this%lak_get_external_outlet(n, rrate)
4350 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4353 if (this%imover == 1)
then
4354 if (this%iboundpak(n) /= 0)
then
4355 rrate = this%pakmvrobj%get_qfrommvr(n)
4359 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4365 do n = 1, this%nlakes
4367 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4371 rrate = -this%simvals(j)
4372 this%qleak(j) = rrate
4373 if (this%iboundpak(n) /= 0)
then
4374 call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4380 call this%lak_fill_budobj()
4381 end subroutine lak_cq
4385 subroutine lak_ot_package_flows(this, icbcfl, ibudfl)
4387 class(laktype) :: this
4388 integer(I4B),
intent(in) :: icbcfl
4389 integer(I4B),
intent(in) :: ibudfl
4390 integer(I4B) :: ibinun
4394 if (this%ibudgetout /= 0)
then
4395 ibinun = this%ibudgetout
4397 if (icbcfl == 0) ibinun = 0
4398 if (ibinun > 0)
then
4399 call this%budobj%save_flows(this%dis, ibinun,
kstp,
kper,
delt, &
4404 if (ibudfl /= 0 .and. this%iprflow /= 0)
then
4405 call this%budobj%write_flowtable(this%dis,
kstp,
kper)
4407 end subroutine lak_ot_package_flows
4411 subroutine lak_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
4412 class(laktype) :: this
4413 integer(I4B),
intent(in) :: icbcfl
4414 integer(I4B),
intent(in) :: ibudfl
4415 integer(I4B),
intent(in) :: icbcun
4416 integer(I4B),
dimension(:),
optional,
intent(in) :: imap
4419 call this%BndType%bnd_ot_model_flows(icbcfl, ibudfl, icbcun, this%imap)
4420 end subroutine lak_ot_model_flows
4424 subroutine lak_ot_dv(this, idvsave, idvprint)
4428 class(laktype) :: this
4429 integer(I4B),
intent(in) :: idvsave
4430 integer(I4B),
intent(in) :: idvprint
4431 integer(I4B) :: ibinun
4441 if (this%istageout /= 0)
then
4442 ibinun = this%istageout
4444 if (idvsave == 0) ibinun = 0
4447 if (ibinun > 0)
then
4448 do n = 1, this%nlakes
4450 d = v - this%lakebot(n)
4451 if (this%iboundpak(n) == 0)
then
4453 else if (d <= dzero)
then
4459 this%nlakes, 1, 1, ibinun)
4463 if (idvprint /= 0 .and. this%iprhed /= 0)
then
4466 call this%stagetab%set_kstpkper(
kstp,
kper)
4469 do n = 1, this%nlakes
4470 if (this%iboundpak(n) == 0)
then
4476 stage = this%xnewpak(n)
4477 call this%lak_calculate_sarea(n, stage, sa)
4478 call this%lak_calculate_warea(n, stage, wa)
4479 call this%lak_calculate_vol(n, stage, v)
4481 if (this%inamedbound == 1)
then
4482 call this%stagetab%add_term(this%lakename(n))
4484 call this%stagetab%add_term(n)
4485 call this%stagetab%add_term(stage)
4486 call this%stagetab%add_term(sa)
4487 call this%stagetab%add_term(wa)
4488 call this%stagetab%add_term(v)
4491 end subroutine lak_ot_dv
4495 subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl)
4499 class(laktype) :: this
4500 integer(I4B),
intent(in) :: kstp
4501 integer(I4B),
intent(in) :: kper
4502 integer(I4B),
intent(in) :: iout
4503 integer(I4B),
intent(in) :: ibudfl
4505 call this%budobj%write_budtable(kstp, kper, iout, ibudfl,
totim,
delt)
4506 end subroutine lak_ot_bdsummary
4510 subroutine lak_da(this)
4514 class(laktype) :: this
4517 deallocate (this%lakename)
4518 deallocate (this%status)
4519 deallocate (this%clakbudget)
4521 deallocate (this%cauxcbc)
4530 if (this%ntables > 0)
then
4539 call this%budobj%budgetobject_da()
4540 deallocate (this%budobj)
4541 nullify (this%budobj)
4544 if (this%noutlets > 0)
then
4557 if (this%iprhed > 0)
then
4558 call this%stagetab%table_da()
4559 deallocate (this%stagetab)
4560 nullify (this%stagetab)
4564 if (this%ipakcsv > 0)
then
4565 if (
associated(this%pakcsvtab))
then
4566 call this%pakcsvtab%table_da()
4567 deallocate (this%pakcsvtab)
4568 nullify (this%pakcsvtab)
4578 if (
allocated(this%pakcsvfile))
deallocate (this%pakcsvfile)
4637 if (this%iimplicit == 0)
then
4682 nullify (this%gwfiss)
4685 call this%BndType%bnd_da()
4686 end subroutine lak_da
4691 subroutine define_listlabel(this)
4693 class(laktype),
intent(inout) :: this
4696 this%listlabel = trim(this%filtyp)//
' NO.'
4697 if (this%dis%ndim == 3)
then
4698 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
4699 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'ROW'
4700 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'COL'
4701 elseif (this%dis%ndim == 2)
then
4702 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
4703 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'CELL2D'
4705 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'NODE'
4707 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'STRESS RATE'
4708 if (this%inamedbound == 1)
then
4709 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'BOUNDARY NAME'
4711 end subroutine define_listlabel
4716 subroutine lak_set_pointers(this, neq, ibound, xnew, xold, flowja)
4720 class(laktype) :: this
4721 integer(I4B),
pointer :: neq
4722 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
4723 real(DP),
dimension(:),
pointer,
contiguous :: xnew
4724 real(DP),
dimension(:),
pointer,
contiguous :: xold
4725 real(DP),
dimension(:),
pointer,
contiguous :: flowja
4728 integer(I4B) :: istart, iend
4731 call this%BndType%set_pointers(neq, ibound, xnew, xold, flowja)
4737 if (this%iimplicit /= 0)
then
4738 istart = this%dis%nodes + this%ioffset + 1
4739 iend = istart + this%nlakes - 1
4740 this%iboundpak => this%ibound(istart:iend)
4741 this%xnewpak => this%xnew(istart:iend)
4742 call mem_checkin(this%xnewpak,
'XNEWPAK', this%memoryPath,
'X', &
4743 this%memoryPathModel)
4746 do n = 1, this%nlakes
4747 this%xnewpak(n) =
dep20
4750 end subroutine lak_set_pointers
4758 subroutine lak_ac(this, moffset, sparse)
4761 class(laktype),
intent(inout) :: this
4762 integer(I4B),
intent(in) :: moffset
4765 integer(I4B) :: j, n
4767 integer(I4B) :: jglo
4768 integer(I4B) :: nglo
4771 if (this%iimplicit == 0)
return
4774 do n = 1, this%nlakes
4775 nglo = moffset + this%dis%nodes + this%ioffset + n
4776 call sparse%addconnection(nglo, nglo, 1)
4777 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4780 call sparse%addconnection(nglo, jglo, 1)
4781 call sparse%addconnection(jglo, nglo, 1)
4784 end subroutine lak_ac
4792 subroutine lak_mc(this, moffset, matrix_sln)
4795 class(laktype),
intent(inout) :: this
4796 integer(I4B),
intent(in) :: moffset
4801 integer(I4B) :: iglo
4802 integer(I4B) :: jglo
4803 integer(I4B) :: ipos
4808 if (this%iimplicit == 0)
then
4809 call mem_allocate(this%idxlocnode, 0,
'IDXLOCNODE', this%memoryPath)
4810 call mem_allocate(this%idxdiag, 0,
'IDXDIAG', this%memoryPath)
4811 call mem_allocate(this%idxoffdglo, 0,
'IDXOFFDGLO', this%memoryPath)
4812 call mem_allocate(this%idxsymdglo, 0,
'IDXSYMDGLO', this%memoryPath)
4813 call mem_allocate(this%idxsymoffdglo, 0,
'IDXSYMOFFDGLO', this%memoryPath)
4816 call mem_allocate(this%idxlocnode, this%nlakes,
'IDXLOCNODE', &
4818 call mem_allocate(this%idxdiag, this%nlakes,
'IDXDIAG', this%memoryPath)
4819 call mem_allocate(this%idxoffdglo, this%maxbound,
'IDXOFFDGLO', &
4821 call mem_allocate(this%idxsymdglo, this%maxbound,
'IDXSYMDGLO', &
4823 call mem_allocate(this%idxsymoffdglo, this%maxbound,
'IDXSYMOFFDGLO', &
4830 do n = 1, this%nlakes
4831 iglo = moffset + this%dis%nodes + this%ioffset + n
4832 this%idxlocnode(n) = this%dis%nodes + this%ioffset + n
4833 this%idxdiag(n) = matrix_sln%get_position_diag(iglo)
4834 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4835 jglo = this%cellid(j) + moffset
4836 this%idxoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
4843 do n = 1, this%nlakes
4844 jglo = moffset + this%dis%nodes + this%ioffset + n
4845 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4846 iglo = this%cellid(j) + moffset
4847 this%idxsymdglo(ipos) = matrix_sln%get_position_diag(iglo)
4848 this%idxsymoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
4852 end subroutine lak_mc
4859 logical function lak_obs_supported(this)
4861 class(laktype) :: this
4863 lak_obs_supported = .true.
4864 end function lak_obs_supported
4869 subroutine lak_df_obs(this)
4871 class(laktype) :: this
4873 integer(I4B) :: indx
4877 call this%obs%StoreObsType(
'stage', .false., indx)
4878 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4882 call this%obs%StoreObsType(
'ext-inflow', .true., indx)
4883 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4887 call this%obs%StoreObsType(
'outlet-inflow', .true., indx)
4888 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4892 call this%obs%StoreObsType(
'inflow', .true., indx)
4893 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4897 call this%obs%StoreObsType(
'from-mvr', .true., indx)
4898 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4902 call this%obs%StoreObsType(
'rainfall', .true., indx)
4903 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4907 call this%obs%StoreObsType(
'runoff', .true., indx)
4908 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4912 call this%obs%StoreObsType(
'lak', .true., indx)
4913 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4917 call this%obs%StoreObsType(
'evaporation', .true., indx)
4918 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4922 call this%obs%StoreObsType(
'withdrawal', .true., indx)
4923 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4927 call this%obs%StoreObsType(
'ext-outflow', .true., indx)
4928 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4932 call this%obs%StoreObsType(
'to-mvr', .true., indx)
4933 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4937 call this%obs%StoreObsType(
'storage', .true., indx)
4938 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4942 call this%obs%StoreObsType(
'constant', .true., indx)
4943 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4947 call this%obs%StoreObsType(
'outlet', .true., indx)
4948 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4952 call this%obs%StoreObsType(
'volume', .true., indx)
4953 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4957 call this%obs%StoreObsType(
'surface-area', .true., indx)
4958 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4962 call this%obs%StoreObsType(
'wetted-area', .true., indx)
4963 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4967 call this%obs%StoreObsType(
'conductance', .true., indx)
4968 this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4969 end subroutine lak_df_obs
4974 subroutine lak_bd_obs(this)
4976 class(laktype) :: this
4979 integer(I4B) :: igwfnode
4990 if (this%obs%npakobs > 0)
then
4991 call this%obs%obs_bd_clear()
4992 do i = 1, this%obs%npakobs
4993 obsrv => this%obs%pakobs(i)%obsrv
4994 do j = 1, obsrv%indxbnds_count
4996 jj = obsrv%indxbnds(j)
4997 select case (obsrv%ObsTypeId)
4999 if (this%iboundpak(jj) /= 0)
then
5000 v = this%xnewpak(jj)
5003 if (this%iboundpak(jj) /= 0)
then
5004 call this%lak_calculate_inflow(jj, v)
5006 case (
'OUTLET-INFLOW')
5007 if (this%iboundpak(jj) /= 0)
then
5008 call this%lak_calculate_outlet_inflow(jj, v)
5011 if (this%iboundpak(jj) /= 0)
then
5012 call this%lak_calculate_inflow(jj, v)
5013 call this%lak_calculate_outlet_inflow(jj, v2)
5017 if (this%iboundpak(jj) /= 0)
then
5018 if (this%imover == 1)
then
5019 v = this%pakmvrobj%get_qfrommvr(jj)
5023 if (this%iboundpak(jj) /= 0)
then
5027 if (this%iboundpak(jj) /= 0)
then
5032 if (this%iboundpak(n) /= 0)
then
5033 igwfnode = this%cellid(jj)
5034 hgwf = this%xnew(igwfnode)
5035 if (this%hcof(jj) /=
dzero)
then
5036 v = -(this%hcof(jj) * (this%xnewpak(n) - hgwf))
5041 case (
'EVAPORATION')
5042 if (this%iboundpak(jj) /= 0)
then
5046 if (this%iboundpak(jj) /= 0)
then
5049 case (
'EXT-OUTFLOW')
5051 if (this%iboundpak(n) /= 0)
then
5052 if (this%lakeout(jj) == 0)
then
5053 v = this%simoutrate(jj)
5055 if (this%imover == 1)
then
5056 v = v + this%pakmvrobj%get_qtomvr(jj)
5063 if (this%iboundpak(n) /= 0)
then
5064 if (this%imover == 1)
then
5065 v = this%pakmvrobj%get_qtomvr(jj)
5072 if (this%iboundpak(jj) /= 0)
then
5076 if (this%iboundpak(jj) /= 0)
then
5081 if (this%iboundpak(n) /= 0)
then
5082 v = this%simoutrate(jj)
5085 if (this%iboundpak(jj) /= 0)
then
5086 call this%lak_calculate_vol(jj, this%xnewpak(jj), v)
5088 case (
'SURFACE-AREA')
5089 if (this%iboundpak(jj) /= 0)
then
5090 hlak = this%xnewpak(jj)
5091 call this%lak_calculate_sarea(jj, hlak, v)
5093 case (
'WETTED-AREA')
5095 if (this%iboundpak(n) /= 0)
then
5096 hlak = this%xnewpak(n)
5097 igwfnode = this%cellid(jj)
5098 hgwf = this%xnew(igwfnode)
5099 call this%lak_calculate_conn_warea(n, jj, hlak, hgwf, v)
5101 case (
'CONDUCTANCE')
5103 if (this%iboundpak(n) /= 0)
then
5104 hlak = this%xnewpak(n)
5105 igwfnode = this%cellid(jj)
5106 hgwf = this%xnew(igwfnode)
5107 call this%lak_calculate_conn_conductance(n, jj, hlak, hgwf, v)
5110 errmsg =
'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
5113 call this%obs%SaveOneSimval(obsrv, v)
5122 end subroutine lak_bd_obs
5129 subroutine lak_rp_obs(this)
5132 class(laktype),
intent(inout) :: this
5139 character(len=LENBOUNDNAME) :: bname
5140 logical(LGP) :: jfound
5143 10
format(
'Boundary "', a,
'" for observation "', a, &
5144 '" is invalid in package "', a,
'"')
5150 do i = 1, this%obs%npakobs
5151 obsrv => this%obs%pakobs(i)%obsrv
5154 nn1 = obsrv%NodeNumber
5156 bname = obsrv%FeatureName
5157 if (bname /=
'')
then
5162 if (obsrv%ObsTypeId ==
'LAK' .or. &
5163 obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5164 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5165 do j = 1, this%nlakes
5166 do jj = this%idxlakeconn(j), this%idxlakeconn(j + 1) - 1
5167 if (this%boundname(jj) == bname)
then
5169 call obsrv%AddObsIndex(jj)
5173 else if (obsrv%ObsTypeId ==
'EXT-OUTFLOW' .or. &
5174 obsrv%ObsTypeId ==
'TO-MVR' .or. &
5175 obsrv%ObsTypeId ==
'OUTLET')
then
5176 do j = 1, this%noutlets
5178 if (this%lakename(jj) == bname)
then
5180 call obsrv%AddObsIndex(j)
5184 do j = 1, this%nlakes
5185 if (this%lakename(j) == bname)
then
5187 call obsrv%AddObsIndex(j)
5191 if (.not. jfound)
then
5193 trim(bname), trim(obsrv%Name), trim(this%packName)
5198 if (obsrv%indxbnds_count == 0)
then
5199 if (obsrv%ObsTypeId ==
'LAK' .or. &
5200 obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5201 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5202 nn2 = obsrv%NodeNumber2
5203 j = this%idxlakeconn(nn1) + nn2 - 1
5204 call obsrv%AddObsIndex(j)
5206 call obsrv%AddObsIndex(nn1)
5209 errmsg =
'Programming error in lak_rp_obs'
5216 if (obsrv%ObsTypeId ==
'STAGE')
then
5217 if (obsrv%indxbnds_count > 1)
then
5218 write (
errmsg,
'(a,3(1x,a))') &
5219 trim(adjustl(obsrv%ObsTypeId)), &
5220 'for observation', trim(adjustl(obsrv%Name)), &
5221 ' must be assigned to a lake with a unique boundname.'
5227 if (obsrv%ObsTypeId ==
'TO-MVR' .or. &
5228 obsrv%ObsTypeId ==
'EXT-OUTFLOW' .or. &
5229 obsrv%ObsTypeId ==
'OUTLET')
then
5230 do j = 1, obsrv%indxbnds_count
5231 nn1 = obsrv%indxbnds(j)
5232 if (nn1 < 1 .or. nn1 > this%noutlets)
then
5233 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
5234 trim(adjustl(obsrv%ObsTypeId)), &
5235 ' outlet must be > 0 and <=', this%noutlets, &
5236 '(specified value is ', nn1,
')'
5240 else if (obsrv%ObsTypeId ==
'LAK' .or. &
5241 obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5242 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5243 do j = 1, obsrv%indxbnds_count
5244 nn1 = obsrv%indxbnds(j)
5245 if (nn1 < 1 .or. nn1 > this%maxbound)
then
5246 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
5247 trim(adjustl(obsrv%ObsTypeId)), &
5248 'lake connection number must be > 0 and <=', this%maxbound, &
5249 '(specified value is ', nn1,
')'
5254 do j = 1, obsrv%indxbnds_count
5255 nn1 = obsrv%indxbnds(j)
5256 if (nn1 < 1 .or. nn1 > this%nlakes)
then
5257 write (
errmsg,
'(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
5258 trim(adjustl(obsrv%ObsTypeId)), &
5259 ' lake must be > 0 and <=', this%nlakes, &
5260 '(specified value is ', nn1,
')'
5272 end subroutine lak_rp_obs
5281 subroutine lak_process_obsid(obsrv, dis, inunitobs, iout)
5285 integer(I4B),
intent(in) :: inunitobs
5286 integer(I4B),
intent(in) :: iout
5288 integer(I4B) :: nn1, nn2
5289 integer(I4B) :: icol, istart, istop
5290 character(len=LINELENGTH) :: string
5291 character(len=LENBOUNDNAME) :: bndname
5293 string = obsrv%IDstring
5301 obsrv%FeatureName = bndname
5303 if (obsrv%ObsTypeId ==
'LAK' .or. obsrv%ObsTypeId ==
'CONDUCTANCE' .or. &
5304 obsrv%ObsTypeId ==
'WETTED-AREA')
then
5306 if (len_trim(bndname) < 1 .and. nn2 < 0)
then
5307 write (
errmsg,
'(a,1x,a,a,1x,a,1x,a)') &
5308 'For observation type', trim(adjustl(obsrv%ObsTypeId)), &
5309 ', ID given as an integer and not as boundname,', &
5310 'but ID2 (iconn) is missing. Either change ID to valid', &
5311 'boundname or supply valid entry for ID2.'
5315 obsrv%FeatureName = bndname
5319 obsrv%NodeNumber2 = nn2
5324 obsrv%NodeNumber = nn1
5325 end subroutine lak_process_obsid
5333 subroutine lak_accumulate_chterm(this, ilak, rrate, chratin, chratout)
5335 class(laktype) :: this
5336 integer(I4B),
intent(in) :: ilak
5337 real(DP),
intent(in) :: rrate
5338 real(DP),
intent(inout) :: chratin
5339 real(DP),
intent(inout) :: chratout
5344 if (this%iboundpak(ilak) < 0)
then
5346 this%chterm(ilak) = this%chterm(ilak) + q
5352 chratout = chratout - q
5356 chratin = chratin + q
5359 end subroutine lak_accumulate_chterm
5363 subroutine lak_bound_update(this)
5365 class(laktype),
intent(inout) :: this
5367 integer(I4B) :: j, n, node
5368 real(DP) :: hlak, head, clak
5371 if (this%nbound == 0)
return
5374 do n = 1, this%nlakes
5375 hlak = this%xnewpak(n)
5376 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5377 node = this%cellid(j)
5378 head = this%xnew(node)
5379 call this%lak_calculate_conn_conductance(n, j, hlak, head, clak)
5380 this%bound(1, j) = hlak
5381 this%bound(2, j) = clak
5384 end subroutine lak_bound_update
5394 subroutine lak_solve(this, update, only_fallback)
5398 class(laktype),
intent(inout) :: this
5399 logical(LGP),
intent(in),
optional :: update
5400 logical(LGP),
intent(in),
optional :: only_fallback
5402 logical(LGP) :: lupdate
5403 logical(LGP) :: fbonly
5406 integer(I4B) :: iicnvg
5407 integer(I4B) :: iter
5408 integer(I4B) :: maxiter
5409 integer(I4B) :: ncnv
5421 if (
present(update))
then
5428 if (
present(only_fallback))
then
5429 fbonly = only_fallback
5438 do n = 1, this%nlakes
5442 if (fbonly .and. this%ifallback(n) == 0)
then
5447 this%surfin(n) =
dzero
5448 this%surfout(n) =
dzero
5449 this%surfout1(n) =
dzero
5450 if (this%xnewpak(n) < this%lakebot(n))
then
5451 this%xnewpak(n) = this%lakebot(n)
5453 if (this%gwfiss /= 0)
then
5454 this%xoldpak(n) = this%xnewpak(n)
5459 this%en1(n) = this%lakebot(n)
5460 call this%lak_calculate_residual(n, this%en1(n), this%r1(n))
5461 this%en2(n) = this%laketop(n)
5462 call this%lak_calculate_residual(n, this%en2(n), this%r2(n))
5464 do n = 1, this%noutlets
5465 this%simoutrate(n) =
dzero
5469 do n = 1, this%nlakes
5470 call this%lak_calculate_outlet_inflow(n, this%surfin(n))
5475 do n = 1, this%nlakes
5476 hlak0 = this%xoldpak(n)
5477 hlak = this%xnewpak(n)
5478 call this%lak_calculate_runoff(n, ro)
5479 call this%lak_calculate_inflow(n, qinf)
5480 call this%lak_calculate_external(n, ex)
5481 call this%lak_calculate_vol(n, hlak0, v0)
5482 call this%lak_calculate_vol(n, hlak, v1)
5483 this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
5488 do n = 1, this%nlakes
5489 call this%lak_calculate_outlet_inflow(n, outinf)
5490 this%flwin(n) = this%flwin(n) + outinf
5494 maxiter = this%maxlakit
5497 converge:
do iter = 1, maxiter
5499 do n = 1, this%nlakes
5500 if (this%ncncvr(n) == 0) ncnv = 1
5502 if (iter == maxiter) ncnv = 0
5503 if (ncnv == 0) iicnvg = 1
5506 do n = 1, this%nlakes
5507 this%evap(n) =
dzero
5508 this%precip(n) =
dzero
5509 this%precip1(n) =
dzero
5510 this%seep(n) =
dzero
5511 this%seep1(n) =
dzero
5512 this%evap(n) =
dzero
5513 this%evap1(n) =
dzero
5514 this%evapo(n) =
dzero
5515 this%withr(n) =
dzero
5516 this%withr1(n) =
dzero
5517 this%flwiter(n) = this%flwin(n)
5518 this%flwiter1(n) = this%flwin(n)
5519 if (this%gwfiss /= 0)
then
5520 this%flwiter(n) =
dep20
5521 this%flwiter1(n) =
dep20
5523 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5524 this%hcof(j) =
dzero
5529 do n = 1, this%nlakes
5530 if (fbonly .and. this%ifallback(n) == 0) cycle
5531 call this%lak_estimate_seepage_single(n, ncnv)
5534 laklevel:
do n = 1, this%nlakes
5535 if (fbonly .and. this%ifallback(n) == 0) cycle laklevel
5536 call this%lak_solve_single(n, iter, maxiter, ncnv, lupdate)
5539 if (iicnvg == 1)
exit converge
5546 if (this%imover == 1)
then
5547 do n = 1, this%noutlets
5548 call this%pakmvrobj%accumulate_qformvr(n, -this%simoutrate(n))
5551 end subroutine lak_solve
5561 subroutine lak_solve_single(this, n, iter, maxiter, ncnv, lupdate)
5565 class(laktype),
intent(inout) :: this
5566 integer(I4B),
intent(in) :: n
5567 integer(I4B),
intent(in) :: iter
5568 integer(I4B),
intent(in) :: maxiter
5569 integer(I4B),
intent(in) :: ncnv
5570 logical(LGP),
intent(in) :: lupdate
5572 integer(I4B) :: ibflg
5573 integer(I4B) :: idhp
5593 real(DP) :: qtolfact
5599 if (this%iboundpak(n) == 0)
then
5604 hlak = this%xnewpak(n)
5605 if (iter < maxiter)
then
5606 this%stageiter(n) = this%xnewpak(n)
5608 call this%lak_calculate_rainfall(n, hlak, ra)
5610 this%flwiter(n) = this%flwiter(n) + ra
5611 call this%lak_calculate_rainfall(n, hlak + delh, ra)
5612 this%precip1(n) = ra
5613 this%flwiter1(n) = this%flwiter1(n) + ra
5616 call this%lak_calculate_withdrawal(n, this%flwiter(n), wr)
5618 call this%lak_calculate_withdrawal(n, this%flwiter1(n), wr)
5622 call this%lak_calculate_evaporation(n, hlak, this%flwiter(n), ev)
5624 call this%lak_calculate_evaporation(n, hlak + delh, this%flwiter1(n), ev)
5628 call this%lak_calculate_outlet_outflow(n, hlak + delh, &
5631 call this%lak_calculate_outlet_outflow(n, hlak, this%flwiter(n), &
5635 call this%lak_calculate_outlet_inflow(n, this%surfin(n))
5639 if (this%iboundpak(n) > 0 .and. lupdate .eqv. .true.)
then
5642 hlak0 = this%xoldpak(n)
5643 hlak = this%xnewpak(n)
5644 call this%lak_calculate_vol(n, hlak0, v0)
5645 call this%lak_calculate_vol(n, hlak, v1)
5646 call this%lak_calculate_runoff(n, ro)
5647 call this%lak_calculate_inflow(n, qinf)
5648 call this%lak_calculate_external(n, ex)
5649 this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
5653 resid = this%precip(n) + this%evap(n) + this%withr(n) + ro + &
5654 qinf + ex + this%surfin(n) + &
5655 this%surfout(n) + this%seep(n)
5656 resid1 = this%precip1(n) + this%evap1(n) + this%withr1(n) + ro + &
5657 qinf + ex + this%surfin(n) + &
5658 this%surfout1(n) + this%seep1(n)
5661 hlak = this%xnewpak(n)
5662 if (this%gwfiss /= 1)
then
5663 call this%lak_calculate_vol(n, hlak, v1)
5664 resid = resid + (v0 - v1) /
delt
5665 call this%lak_calculate_vol(n, hlak + delh, v1)
5666 resid1 = resid1 + (v0 - v1) /
delt
5670 if (abs(resid1 - resid) >
dzero)
then
5671 derv = (resid1 - resid) / delh
5673 if (abs(derv) >
dprec)
then
5677 if (resid <
dzero)
then
5680 call this%lak_vol2stage(n, resid, dh)
5687 if (iter == 1) this%dh0(n) = dh
5689 adh0 = abs(this%dh0(n))
5690 if ((ts >= this%en2(n)) .or. (ts < this%en1(n)))
then
5693 if ((adh > adh0) .or. (ts - this%lakebot(n)) <
dprec)
then
5695 call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
5701 this%seep0(n) = this%seep(n)
5705 if (this%seep(n) * this%seep0(n) <
dprec)
then
5706 this%iseepc(n) = this%iseepc(n) + 1
5712 if (dh * this%dh0(n) <
dprec) idhp = 1
5715 if (adh > adh0) idhp = 1
5718 this%idhc(n) = this%idhc(n) + 1
5723 if (ibflg == 1)
then
5724 if (this%iseepc(n) > 7 .or. this%idhc(n) > 12)
then
5725 call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
5734 if (hlak < this%lakebot(n))
then
5735 hlak = this%lakebot(n)
5739 call this%lak_calculate_sarea(n, hlak, area)
5742 if (area >
dzero)
then
5743 qtolfact =
delt / area
5749 call this%lak_calculate_residual(n, hlak, resid)
5753 if (abs(dh) < delh .and. abs(resid) * qtolfact < this%dmaxchg)
then
5756 this%xnewpak(n) = hlak
5759 this%seep0(n) = this%seep(n)
5762 end subroutine lak_solve_single
5775 subroutine lak_estimate_seepage_single(this, n, ncnv)
5777 class(laktype),
intent(inout) :: this
5778 integer(I4B),
intent(in) :: n
5779 integer(I4B),
intent(in) :: ncnv
5783 integer(I4B) :: igwfnode
5784 integer(I4B) :: idry
5785 integer(I4B) :: idry1
5797 if (this%iboundpak(n) == 0)
return
5799 estseep:
do i = 1, 2
5801 if (this%gwfiss /= 0)
then
5802 this%xoldpak(n) = this%xnewpak(n)
5804 hlak = this%xnewpak(n)
5805 calcconnseep:
do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5806 igwfnode = this%cellid(j)
5807 head = this%xnew(igwfnode)
5808 if (this%ncncvr(n) /= 2)
then
5809 if (this%ibound(igwfnode) > 0)
then
5810 call this%lak_estimate_conn_exchange(i, n, j, idry, hlak, &
5814 call this%lak_estimate_conn_exchange(i, n, j, idry1, &
5815 hlak + delh, head, qlakgw1, &
5819 if (ncnv == 0 .and. i == 2)
then
5820 if (j == this%maxbound)
then
5824 this%hcof(j) = gwfhcof
5825 this%rhs(j) = gwfrhs
5827 this%hcof(j) =
dzero
5828 this%rhs(j) = qlakgw
5832 this%seep(n) = this%seep(n) + qlakgw
5833 this%seep1(n) = this%seep1(n) + qlakgw1
5839 end subroutine lak_estimate_seepage_single
5845 subroutine lak_bisection(this, n, ibflg, hlak, temporary_stage, dh, residual)
5847 class(laktype),
intent(inout) :: this
5848 integer(I4B),
intent(in) :: n
5849 integer(I4B),
intent(inout) :: ibflg
5850 real(DP),
intent(in) :: hlak
5851 real(DP),
intent(inout) :: temporary_stage
5852 real(DP),
intent(inout) :: dh
5853 real(DP),
intent(inout) :: residual
5856 real(DP) :: temporary_stage0
5857 real(DP) :: residuala
5858 real(DP) :: endpoint1
5859 real(DP) :: endpoint2
5862 temporary_stage0 = hlak
5863 endpoint1 = this%en1(n)
5864 endpoint2 = this%en2(n)
5865 call this%lak_calculate_residual(n, temporary_stage, residuala)
5866 if (hlak > endpoint1 .and. hlak < endpoint2)
then
5869 do i = 1, this%maxlakit
5870 temporary_stage =
dhalf * (endpoint1 + endpoint2)
5871 call this%lak_calculate_residual(n, temporary_stage, residual)
5872 if (abs(residual) ==
dzero .or. &
5873 abs(temporary_stage0 - temporary_stage) < this%dmaxchg)
then
5876 call this%lak_calculate_residual(n, endpoint1, residuala)
5879 if (sign(
done, residuala) == sign(
done, residual))
then
5880 endpoint1 = temporary_stage
5883 endpoint2 = temporary_stage
5885 temporary_stage0 = temporary_stage
5887 dh = hlak - temporary_stage
5888 end subroutine lak_bisection
5893 subroutine lak_calculate_available(this, n, hlak, avail, &
5894 ra, ro, qinf, ex, headp)
5898 class(laktype),
intent(inout) :: this
5899 integer(I4B),
intent(in) :: n
5900 real(DP),
intent(in) :: hlak
5901 real(DP),
intent(inout) :: avail
5902 real(DP),
intent(inout) :: ra
5903 real(DP),
intent(inout) :: ro
5904 real(DP),
intent(inout) :: qinf
5905 real(DP),
intent(inout) :: ex
5906 real(DP),
intent(in),
optional :: headp
5909 integer(I4B) :: idry
5910 integer(I4B) :: igwfnode
5917 if (
present(headp))
then
5927 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5928 igwfnode = this%cellid(j)
5929 if (this%ibound(igwfnode) == 0) cycle
5930 head = this%xnew(igwfnode) + hp
5931 call this%lak_estimate_conn_exchange(1, n, j, idry, hlak, head, qlakgw, &
5936 call this%lak_calculate_rainfall(n, hlak, ra)
5940 call this%lak_calculate_runoff(n, ro)
5944 call this%lak_calculate_inflow(n, qinf)
5945 avail = avail + qinf
5948 call this%lak_calculate_external(n, ex)
5952 call this%lak_calculate_vol(n, this%xoldpak(n), v0)
5953 avail = avail + v0 /
delt
5954 end subroutine lak_calculate_available
5958 subroutine lak_calculate_residual(this, n, hlak, resid, headp)
5962 class(laktype),
intent(inout) :: this
5963 integer(I4B),
intent(in) :: n
5964 real(DP),
intent(in) :: hlak
5965 real(DP),
intent(inout) :: resid
5966 real(DP),
intent(in),
optional :: headp
5969 integer(I4B) :: idry
5970 integer(I4B) :: igwfnode
5989 if (
present(headp))
then
6001 call this%lak_calculate_available(n, hlak, avail, &
6002 ra, ro, qinf, ex, hp)
6005 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
6006 igwfnode = this%cellid(j)
6007 if (this%ibound(igwfnode) == 0) cycle
6008 head = this%xnew(igwfnode) + hp
6009 call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, head, qlakgw, &
6011 seep = seep + qlakgw
6015 call this%lak_calculate_withdrawal(n, avail, wr)
6018 call this%lak_calculate_evaporation(n, hlak, avail, ev)
6021 call this%lak_calculate_outlet_outflow(n, hlak, avail, sout)
6024 call this%lak_calculate_outlet_inflow(n, sin)
6027 resid = ra + ev + wr + ro + qinf + ex + sin + sout + seep
6030 if (this%gwfiss /= 1)
then
6031 hlak0 = this%xoldpak(n)
6032 call this%lak_calculate_vol(n, hlak0, v0)
6033 call this%lak_calculate_vol(n, hlak, v1)
6034 resid = resid + (v0 - v1) /
delt
6036 end subroutine lak_calculate_residual
6040 subroutine lak_setup_budobj(this)
6044 class(laktype) :: this
6046 integer(I4B) :: nbudterm
6047 integer(I4B) :: nlen
6048 integer(I4B) :: j, n, n1, n2
6049 integer(I4B) :: maxlist, naux
6052 character(len=LENBUDTXT) :: text
6053 character(len=LENBUDTXT),
dimension(1) :: auxtxt
6059 do n = 1, this%noutlets
6060 if (this%lakein(n) > 0 .and. this%lakeout(n) > 0)
then
6064 if (nlen > 0) nbudterm = nbudterm + 1
6065 if (this%imover == 1) nbudterm = nbudterm + 2
6066 if (this%naux > 0) nbudterm = nbudterm + 1
6070 call this%budobj%budgetobject_df(this%nlakes, nbudterm, 0, 0, &
6071 ibudcsv=this%ibudcsv)
6077 text =
' FLOW-JA-FACE'
6081 call this%budobj%budterm(idx)%initialize(text, &
6086 maxlist, .false., .false., &
6087 naux, ordered_id1=.false.)
6090 call this%budobj%budterm(idx)%reset(2 * nlen)
6092 do n = 1, this%noutlets
6094 n2 = this%lakeout(n)
6095 if (n1 > 0 .and. n2 > 0)
then
6096 call this%budobj%budterm(idx)%update_term(n1, n2, q)
6097 call this%budobj%budterm(idx)%update_term(n2, n1, -q)
6105 maxlist = this%maxbound
6107 auxtxt(1) =
' FLOW-AREA'
6108 call this%budobj%budterm(idx)%initialize(text, &
6113 maxlist, .false., .true., &
6115 call this%budobj%budterm(idx)%reset(this%maxbound)
6117 do n = 1, this%nlakes
6118 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
6120 call this%budobj%budterm(idx)%update_term(n, n2, q)
6127 maxlist = this%nlakes
6129 call this%budobj%budterm(idx)%initialize(text, &
6134 maxlist, .false., .false., &
6138 text =
' EVAPORATION'
6140 maxlist = this%nlakes
6142 call this%budobj%budterm(idx)%initialize(text, &
6147 maxlist, .false., .false., &
6153 maxlist = this%nlakes
6155 call this%budobj%budterm(idx)%initialize(text, &
6160 maxlist, .false., .false., &
6164 text =
' EXT-INFLOW'
6166 maxlist = this%nlakes
6168 call this%budobj%budterm(idx)%initialize(text, &
6173 maxlist, .false., .false., &
6177 text =
' WITHDRAWAL'
6179 maxlist = this%nlakes
6181 call this%budobj%budterm(idx)%initialize(text, &
6186 maxlist, .false., .false., &
6190 text =
' EXT-OUTFLOW'
6192 maxlist = this%nlakes
6194 call this%budobj%budterm(idx)%initialize(text, &
6199 maxlist, .false., .false., &
6205 maxlist = this%nlakes
6207 auxtxt(1) =
' VOLUME'
6208 call this%budobj%budterm(idx)%initialize(text, &
6213 maxlist, .false., .false., &
6219 maxlist = this%nlakes
6221 call this%budobj%budterm(idx)%initialize(text, &
6226 maxlist, .false., .false., &
6230 if (this%imover == 1)
then
6235 maxlist = this%nlakes
6237 call this%budobj%budterm(idx)%initialize(text, &
6242 maxlist, .false., .false., &
6248 maxlist = this%noutlets
6250 call this%budobj%budterm(idx)%initialize(text, &
6255 maxlist, .false., .false., &
6256 naux, ordered_id1=.false.)
6259 call this%budobj%budterm(idx)%reset(this%noutlets)
6261 do n = 1, this%noutlets
6263 call this%budobj%budterm(idx)%update_term(n1, n1, q)
6274 maxlist = this%nlakes
6275 call this%budobj%budterm(idx)%initialize(text, &
6280 maxlist, .false., .false., &
6285 if (this%iprflow /= 0)
then
6286 call this%budobj%flowtable_df(this%iout)
6288 end subroutine lak_setup_budobj
6292 subroutine lak_fill_budobj(this)
6294 class(laktype) :: this
6296 integer(I4B) :: naux
6297 real(DP),
dimension(:),
allocatable :: auxvartmp
6306 integer(I4B) :: nlen
6309 real(DP) :: lkstg, gwhead, wa
6316 do n = 1, this%noutlets
6317 if (this%lakein(n) > 0 .and. this%lakeout(n) > 0)
then
6323 call this%budobj%budterm(idx)%reset(2 * nlen)
6324 do n = 1, this%noutlets
6326 n2 = this%lakeout(n)
6327 if (n1 > 0 .and. n2 > 0)
then
6328 q = this%simoutrate(n)
6329 if (this%imover == 1)
then
6330 q = q + this%pakmvrobj%get_qtomvr(n)
6332 call this%budobj%budterm(idx)%update_term(n1, n2, q)
6333 call this%budobj%budterm(idx)%update_term(n2, n1, -q)
6340 call this%budobj%budterm(idx)%reset(this%maxbound)
6341 do n = 1, this%nlakes
6342 do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
6345 lkstg = this%xnewpak(n)
6349 gwhead = this%xnew(n2)
6350 call this%lak_calculate_conn_warea(n, j, lkstg, gwhead, wa)
6354 if (this%belev(j) > lkstg) wa =
dzero
6355 this%qauxcbc(1) = wa
6356 call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
6362 call this%budobj%budterm(idx)%reset(this%nlakes)
6363 do n = 1, this%nlakes
6365 call this%budobj%budterm(idx)%update_term(n, n, q)
6370 call this%budobj%budterm(idx)%reset(this%nlakes)
6371 do n = 1, this%nlakes
6373 call this%budobj%budterm(idx)%update_term(n, n, q)
6378 call this%budobj%budterm(idx)%reset(this%nlakes)
6379 do n = 1, this%nlakes
6381 call this%budobj%budterm(idx)%update_term(n, n, q)
6386 call this%budobj%budterm(idx)%reset(this%nlakes)
6387 do n = 1, this%nlakes
6389 call this%budobj%budterm(idx)%update_term(n, n, q)
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
6404 call this%lak_get_external_outlet(n, q)
6406 call this%lak_get_external_mover(n, v)
6408 call this%budobj%budterm(idx)%update_term(n, n, q)
6413 call this%budobj%budterm(idx)%reset(this%nlakes)
6414 do n = 1, this%nlakes
6415 call this%lak_calculate_vol(n, this%xnewpak(n), v1)
6417 this%qauxcbc(1) = v1
6418 call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc)
6423 call this%budobj%budterm(idx)%reset(this%nlakes)
6424 do n = 1, this%nlakes
6426 call this%budobj%budterm(idx)%update_term(n, n, q)
6430 if (this%imover == 1)
then
6434 call this%budobj%budterm(idx)%reset(this%nlakes)
6435 do n = 1, this%nlakes
6436 q = this%pakmvrobj%get_qfrommvr(n)
6437 call this%budobj%budterm(idx)%update_term(n, n, q)
6442 call this%budobj%budterm(idx)%reset(this%noutlets)
6443 do n = 1, this%noutlets
6445 q = this%pakmvrobj%get_qtomvr(n)
6449 call this%budobj%budterm(idx)%update_term(n1, n1, q)
6458 allocate (auxvartmp(naux))
6459 call this%budobj%budterm(idx)%reset(this%nlakes)
6460 do n = 1, this%nlakes
6464 auxvartmp(jj) = this%lauxvar(jj, ii)
6466 call this%budobj%budterm(idx)%update_term(n, n, q, auxvartmp)
6468 deallocate (auxvartmp)
6472 call this%budobj%accumulate_terms()
6473 end subroutine lak_fill_budobj
6480 subroutine lak_setup_tableobj(this)
6484 class(laktype) :: this
6486 integer(I4B) :: nterms
6487 character(len=LINELENGTH) :: title
6488 character(len=LINELENGTH) :: text
6491 if (this%iprhed > 0)
then
6498 if (this%inamedbound == 1)
then
6503 title = trim(adjustl(this%text))//
' PACKAGE ('// &
6504 trim(adjustl(this%packName))//
') STAGES FOR EACH CONTROL VOLUME'
6507 call table_cr(this%stagetab, this%packName, title)
6508 call this%stagetab%table_df(this%nlakes, nterms, this%iout, &
6512 if (this%inamedbound == 1)
then
6514 call this%stagetab%initialize_column(text, 20, alignment=tableft)
6519 call this%stagetab%initialize_column(text, 10, alignment=tabcenter)
6523 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6526 text =
'SURFACE AREA'
6527 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6530 text =
'WETTED AREA'
6531 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6535 call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
6537 end subroutine lak_setup_tableobj
6541 subroutine lak_activate_density(this)
6543 class(laktype),
intent(inout) :: this
6545 integer(I4B) :: i, j
6549 call mem_reallocate(this%denseterms, 3, this%MAXBOUND,
'DENSETERMS', &
6551 do i = 1, this%maxbound
6553 this%denseterms(j, i) =
dzero
6556 write (this%iout,
'(/1x,a)')
'DENSITY TERMS HAVE BEEN ACTIVATED FOR LAKE &
6557 &PACKAGE: '//trim(adjustl(this%packName))
6558 end subroutine lak_activate_density
6564 subroutine lak_activate_viscosity(this)
6568 class(laktype),
intent(inout) :: this
6575 call mem_reallocate(this%viscratios, 2, this%MAXBOUND,
'VISCRATIOS', &
6577 do i = 1, this%maxbound
6579 this%viscratios(j, i) =
done
6582 write (this%iout,
'(/1x,a)')
'VISCOSITY HAS BEEN ACTIVATED FOR LAK &
6583 &PACKAGE: '//trim(adjustl(this%packName))
6584 end subroutine lak_activate_viscosity
6604 subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, &
6605 botl, flow, gwfhcof, gwfrhs)
6607 class(laktype),
intent(inout) :: this
6608 integer(I4B),
intent(in) :: iconn
6609 real(DP),
intent(in) :: stage
6610 real(DP),
intent(in) :: head
6611 real(DP),
intent(in) :: cond
6612 real(DP),
intent(in) :: botl
6613 real(DP),
intent(inout) :: flow
6614 real(DP),
intent(inout) :: gwfhcof
6615 real(DP),
intent(inout) :: gwfrhs
6620 real(DP) :: rdenselak
6621 real(DP) :: rdensegwf
6622 real(DP) :: rdenseavg
6628 logical(LGP) :: stage_below_bot
6629 logical(LGP) :: head_below_bot
6632 if (stage >= botl)
then
6634 stage_below_bot = .false.
6635 rdenselak = this%denseterms(1, iconn)
6638 stage_below_bot = .true.
6639 rdenselak = this%denseterms(2, iconn)
6643 if (head >= botl)
then
6645 head_below_bot = .false.
6646 rdensegwf = this%denseterms(2, iconn)
6649 head_below_bot = .true.
6650 rdensegwf = this%denseterms(1, iconn)
6654 if (rdensegwf ==
dzero)
return
6657 if (stage_below_bot .and. head_below_bot)
then
6664 rdenseavg =
dhalf * (rdenselak + rdensegwf)
6668 d1 = cond * (rdenseavg -
done)
6669 gwfhcof = gwfhcof - d1
6670 gwfrhs = gwfrhs - d1 * ss
6675 if (.not. stage_below_bot .and. .not. head_below_bot)
then
6679 elevgwf = this%denseterms(3, iconn)
6680 if (this%ictype(iconn) == 0 .or. this%ictype(iconn) == 3)
then
6687 elevavg =
dhalf * (elevlak + elevgwf)
6688 havg =
dhalf * (hh + ss)
6689 d2 = cond * (havg - elevavg) * (rdensegwf - rdenselak)
6690 gwfrhs = gwfrhs + d2
6694 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_estimate_seepage_single(this, n, ncnv)
Estimate the lakebed seepage for a single lake.
subroutine lak_solve(this, update, only_fallback)
Solve for lake stage.
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).