40 integer(I4B),
pointer :: iname => null()
41 character(len=24),
dimension(:),
pointer :: aname => null()
42 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
43 real(dp),
dimension(:),
pointer,
contiguous :: hnew => null()
44 integer(I4B),
pointer :: ixt3d => null()
45 integer(I4B),
pointer :: ixt3drhs => null()
46 integer(I4B),
pointer :: iperched => null()
47 integer(I4B),
pointer :: ivarcv => null()
48 integer(I4B),
pointer :: idewatcv => null()
49 integer(I4B),
pointer :: ithickstrt => null()
50 integer(I4B),
pointer :: igwfnewtonur => null()
51 integer(I4B),
pointer :: icalcspdis => null()
52 integer(I4B),
pointer :: isavspdis => null()
53 integer(I4B),
pointer :: isavsat => null()
54 real(dp),
pointer :: hnoflo => null()
55 real(dp),
pointer :: satomega => null()
56 integer(I4B),
pointer :: irewet => null()
57 integer(I4B),
pointer :: iwetit => null()
58 integer(I4B),
pointer :: ihdwet => null()
59 integer(I4B),
pointer :: icellavg => null()
60 real(dp),
pointer :: wetfct => null()
61 real(dp),
pointer :: hdry => null()
62 integer(I4B),
dimension(:),
pointer,
contiguous :: icelltype => null()
63 integer(I4B),
dimension(:),
pointer,
contiguous :: ithickstartflag => null()
66 real(dp),
dimension(:),
pointer,
contiguous :: k11 => null()
67 real(dp),
dimension(:),
pointer,
contiguous :: k22 => null()
68 real(dp),
dimension(:),
pointer,
contiguous :: k33 => null()
69 real(dp),
dimension(:),
pointer,
contiguous :: k11input => null()
70 real(dp),
dimension(:),
pointer,
contiguous :: k22input => null()
71 real(dp),
dimension(:),
pointer,
contiguous :: k33input => null()
72 integer(I4B),
pointer :: iavgkeff => null()
73 integer(I4B),
pointer :: ik22 => null()
74 integer(I4B),
pointer :: ik33 => null()
75 integer(I4B),
pointer :: ik22overk => null()
76 integer(I4B),
pointer :: ik33overk => null()
77 integer(I4B),
pointer :: iangle1 => null()
78 integer(I4B),
pointer :: iangle2 => null()
79 integer(I4B),
pointer :: iangle3 => null()
80 real(dp),
dimension(:),
pointer,
contiguous :: angle1 => null()
81 real(dp),
dimension(:),
pointer,
contiguous :: angle2 => null()
82 real(dp),
dimension(:),
pointer,
contiguous :: angle3 => null()
84 integer(I4B),
pointer :: iwetdry => null()
85 real(dp),
dimension(:),
pointer,
contiguous :: wetdry => null()
86 real(dp),
dimension(:),
pointer,
contiguous :: sat => null()
87 real(dp),
dimension(:),
pointer,
contiguous :: condsat => null()
88 integer(I4B),
dimension(:),
pointer,
contiguous :: ibotnode => null()
90 real(dp),
dimension(:, :),
pointer,
contiguous :: spdis => null()
91 integer(I4B),
pointer :: nedges => null()
92 integer(I4B),
pointer :: lastedge => null()
93 integer(I4B),
dimension(:),
pointer,
contiguous :: nodedge => null()
94 integer(I4B),
dimension(:),
pointer,
contiguous :: ihcedge => null()
95 real(dp),
dimension(:, :),
pointer,
contiguous :: propsedge => null()
96 integer(I4B),
dimension(:),
pointer,
contiguous :: iedge_ptr => null()
97 integer(I4B),
dimension(:),
pointer,
contiguous :: edge_idxs => null()
100 integer(I4B),
pointer :: intvk => null()
101 integer(I4B),
pointer :: invsc => null()
103 integer(I4B),
pointer :: kchangeper => null()
104 integer(I4B),
pointer :: kchangestp => null()
105 integer(I4B),
dimension(:),
pointer,
contiguous :: nodekchange => null()
157 subroutine npf_cr(npfobj, name_model, input_mempath, inunit, iout)
163 character(len=*),
intent(in) :: name_model
164 character(len=*),
intent(in) :: input_mempath
165 integer(I4B),
intent(in) :: inunit
166 integer(I4B),
intent(in) :: iout
168 character(len=*),
parameter :: fmtheader = &
169 "(1x, /1x, 'NPF -- NODE PROPERTY FLOW PACKAGE, VERSION 1, 3/30/2015', &
170 &' INPUT READ FROM MEMPATH: ', A, /)"
176 call npfobj%set_names(1, name_model,
'NPF',
'NPF', input_mempath)
179 call npfobj%allocate_scalars()
182 npfobj%inunit = inunit
189 write (iout, fmtheader) input_mempath
193 allocate (npfobj%spdis_wa)
204 subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options)
212 integer(I4B),
intent(in) :: ingnc
213 integer(I4B),
intent(in) :: invsc
220 if (invsc > 0) this%invsc = invsc
222 if (.not.
present(npf_options))
then
225 call this%source_options()
228 call this%allocate_arrays(this%dis%nodes, this%dis%njas)
231 call this%source_griddata()
232 call this%prepcheck()
234 call this%set_options(npf_options)
237 call this%allocate_arrays(this%dis%nodes, this%dis%njas)
240 call this%check_options()
244 if (this%ixt3d /= 0) xt3d%ixt3d = this%ixt3d
245 call this%xt3d%xt3d_df(dis)
248 if (this%ixt3d /= 0 .and. ingnc > 0)
then
249 call store_error(
'Error in model '//trim(this%name_model)// &
250 '. The XT3D option cannot be used with the GNC &
251 &Package.', terminate=.true.)
262 integer(I4B),
intent(in) :: moffset
266 if (this%ixt3d /= 0)
call this%xt3d%xt3d_ac(moffset, sparse)
271 subroutine npf_mc(this, moffset, matrix_sln)
274 integer(I4B),
intent(in) :: moffset
277 if (this%ixt3d /= 0)
call this%xt3d%xt3d_mc(moffset, matrix_sln)
285 subroutine npf_ar(this, ic, vsc, ibound, hnew)
290 type(
gwfictype),
pointer,
intent(in) :: ic
292 integer(I4B),
dimension(:),
pointer,
contiguous,
intent(inout) :: ibound
293 real(DP),
dimension(:),
pointer,
contiguous,
intent(inout) :: hnew
299 this%ibound => ibound
302 if (this%icalcspdis == 1)
then
303 call mem_reallocate(this%spdis, 3, this%dis%nodes,
'SPDIS', this%memoryPath)
304 call mem_reallocate(this%nodedge, this%nedges,
'NODEDGE', this%memoryPath)
305 call mem_reallocate(this%ihcedge, this%nedges,
'IHCEDGE', this%memoryPath)
306 call mem_reallocate(this%propsedge, 5, this%nedges,
'PROPSEDGE', &
309 'NREDGESNODE', this%memoryPath)
311 'EDGEIDXS', this%memoryPath)
313 do n = 1, this%nedges
314 this%edge_idxs(n) = 0
316 do n = 1, this%dis%nodes
317 this%iedge_ptr(n) = 0
318 this%spdis(:, n) =
dzero
323 if (this%invsc /= 0)
then
328 if (this%invsc > 0)
then
340 call this%store_original_k_arrays(this%dis%nodes, this%dis%njas)
344 call this%preprocess_input()
347 if (this%ixt3d /= 0)
then
348 call this%xt3d%xt3d_ar(ibound, this%k11, this%ik33, this%k33, &
349 this%sat, this%ik22, this%k22, &
350 this%iangle1, this%iangle2, this%iangle3, &
351 this%angle1, this%angle2, this%angle3, &
352 this%inewton, this%icelltype)
356 if (this%intvk /= 0)
then
357 call this%tvk%ar(this%dis)
371 if (this%intvk /= 0)
then
380 subroutine npf_ad(this, nodes, hold, hnew, irestore)
387 integer(I4B),
intent(in) :: nodes
388 real(DP),
dimension(nodes),
intent(inout) :: hold
389 real(DP),
dimension(nodes),
intent(inout) :: hnew
390 integer(I4B),
intent(in) :: irestore
395 if (this%irewet > 0)
then
396 do n = 1, this%dis%nodes
397 if (this%wetdry(n) ==
dzero) cycle
398 if (this%ibound(n) /= 0) cycle
399 hold(n) = this%dis%bot(n)
403 do n = 1, this%dis%nodes
404 if (this%wetdry(n) ==
dzero) cycle
405 if (this%ibound(n) /= 0) cycle
411 if (this%intvk /= 0)
then
417 if (this%invsc /= 0)
then
418 call this%vsc%update_k_with_vsc()
422 if (this%kchangeper ==
kper .and. this%kchangestp ==
kstp)
then
423 if (this%ixt3d == 0)
then
427 do n = 1, this%dis%nodes
428 if (this%nodekchange(n) == 1)
then
429 call this%calc_condsat(n, .false.)
435 if (this%xt3d%lamatsaved .and. .not. this%xt3d%ldispersion)
then
436 call this%xt3d%xt3d_fcpc(this%dis%nodes, .true.)
444 subroutine npf_cf(this, kiter, nodes, hnew)
447 integer(I4B) :: kiter
448 integer(I4B),
intent(in) :: nodes
449 real(DP),
intent(inout),
dimension(nodes) :: hnew
455 if (this%inewton /= 1)
then
456 call this%wd(kiter, hnew)
460 do n = 1, this%dis%nodes
461 if (this%icelltype(n) /= 0)
then
462 if (this%ibound(n) == 0)
then
465 call this%thksat(n, hnew(n), satn)
474 subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
479 integer(I4B) :: kiter
481 integer(I4B),
intent(in),
dimension(:) :: idxglo
482 real(DP),
intent(inout),
dimension(:) :: rhs
483 real(DP),
intent(inout),
dimension(:) :: hnew
485 integer(I4B) :: n, m, ii, idiag, ihc
486 integer(I4B) :: isymcon, idiagm
492 if (this%ixt3d /= 0)
then
493 call this%xt3d%xt3d_fc(kiter, matrix_sln, idxglo, rhs, hnew)
495 do n = 1, this%dis%nodes
496 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
497 if (this%dis%con%mask(ii) == 0) cycle
499 m = this%dis%con%ja(ii)
504 ihc = this%dis%con%ihc(this%dis%con%jas(ii))
505 hyn = this%hy_eff(n, m, ihc, ipos=ii)
506 hym = this%hy_eff(m, n, ihc, ipos=ii)
509 if (ihc == c3d_vertical)
then
512 cond =
vcond(this%ibound(n), this%ibound(m), &
513 this%icelltype(n), this%icelltype(m), this%inewton, &
514 this%ivarcv, this%idewatcv, &
515 this%condsat(this%dis%con%jas(ii)), hnew(n), hnew(m), &
517 this%sat(n), this%sat(m), &
518 this%dis%top(n), this%dis%top(m), &
519 this%dis%bot(n), this%dis%bot(m), &
520 this%dis%con%hwva(this%dis%con%jas(ii)))
523 if (this%iperched /= 0)
then
524 if (this%icelltype(m) /= 0)
then
525 if (hnew(m) < this%dis%top(m))
then
528 idiag = this%dis%con%ia(n)
529 rhs(n) = rhs(n) - cond * this%dis%bot(n)
530 call matrix_sln%add_value_pos(idxglo(idiag), -cond)
533 isymcon = this%dis%con%isym(ii)
534 call matrix_sln%add_value_pos(idxglo(isymcon), cond)
535 rhs(m) = rhs(m) + cond * this%dis%bot(n)
546 cond =
hcond(this%ibound(n), this%ibound(m), &
547 this%icelltype(n), this%icelltype(m), &
549 this%dis%con%ihc(this%dis%con%jas(ii)), &
551 this%condsat(this%dis%con%jas(ii)), &
552 hnew(n), hnew(m), this%sat(n), this%sat(m), hyn, hym, &
553 this%dis%top(n), this%dis%top(m), &
554 this%dis%bot(n), this%dis%bot(m), &
555 this%dis%con%cl1(this%dis%con%jas(ii)), &
556 this%dis%con%cl2(this%dis%con%jas(ii)), &
557 this%dis%con%hwva(this%dis%con%jas(ii)))
561 idiag = this%dis%con%ia(n)
562 call matrix_sln%add_value_pos(idxglo(ii), cond)
563 call matrix_sln%add_value_pos(idxglo(idiag), -cond)
566 isymcon = this%dis%con%isym(ii)
567 idiagm = this%dis%con%ia(m)
568 call matrix_sln%add_value_pos(idxglo(isymcon), cond)
569 call matrix_sln%add_value_pos(idxglo(idiagm), -cond)
578 subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew)
581 integer(I4B) :: kiter
583 integer(I4B),
intent(in),
dimension(:) :: idxglo
584 real(DP),
intent(inout),
dimension(:) :: rhs
585 real(DP),
intent(inout),
dimension(:) :: hnew
587 integer(I4B) :: nodes, nja
588 integer(I4B) :: n, m, ii, idiag
589 integer(I4B) :: isymcon, idiagm
594 real(DP) :: filledterm
602 nodes = this%dis%nodes
603 nja = this%dis%con%nja
604 if (this%ixt3d /= 0)
then
605 call this%xt3d%xt3d_fn(kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew)
609 idiag = this%dis%con%ia(n)
610 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
611 if (this%dis%con%mask(ii) == 0) cycle
613 m = this%dis%con%ja(ii)
614 isymcon = this%dis%con%isym(ii)
617 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0 .and. &
618 this%ivarcv == 0)
then
624 if (hnew(m) < hnew(n)) iups = n
626 if (iups == n) idn = m
629 if (this%icelltype(iups) == 0) cycle
633 topup = this%dis%top(iups)
634 botup = this%dis%bot(iups)
635 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 2)
then
636 topup = min(this%dis%top(n), this%dis%top(m))
637 botup = max(this%dis%bot(n), this%dis%bot(m))
641 cond = this%condsat(this%dis%con%jas(ii))
644 consterm = -cond * (hnew(iups) - hnew(idn))
646 filledterm = matrix_sln%get_value_pos(idxglo(ii))
649 idiagm = this%dis%con%ia(m)
654 term = consterm * derv
655 rhs(n) = rhs(n) + term * hnew(n)
656 rhs(m) = rhs(m) - term * hnew(n)
658 call matrix_sln%add_value_pos(idxglo(idiag), term)
660 if (this%ibound(n) > 0)
then
661 filledterm = matrix_sln%get_value_pos(idxglo(ii))
662 call matrix_sln%set_value_pos(idxglo(ii), filledterm)
665 filledterm = matrix_sln%get_value_pos(idxglo(idiagm))
666 call matrix_sln%set_value_pos(idxglo(idiagm), filledterm)
668 if (this%ibound(m) > 0)
then
669 call matrix_sln%add_value_pos(idxglo(isymcon), -term)
674 term = -consterm * derv
675 rhs(n) = rhs(n) + term * hnew(m)
676 rhs(m) = rhs(m) - term * hnew(m)
678 filledterm = matrix_sln%get_value_pos(idxglo(idiag))
679 call matrix_sln%set_value_pos(idxglo(idiag), filledterm)
681 if (this%ibound(n) > 0)
then
682 call matrix_sln%add_value_pos(idxglo(ii), term)
685 call matrix_sln%add_value_pos(idxglo(idiagm), -term)
687 if (this%ibound(m) > 0)
then
688 filledterm = matrix_sln%get_value_pos(idxglo(isymcon))
689 call matrix_sln%set_value_pos(idxglo(isymcon), filledterm)
705 subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
708 integer(I4B),
intent(in) :: neqmod
709 real(DP),
dimension(neqmod),
intent(inout) :: x
710 real(DP),
dimension(neqmod),
intent(in) :: xtemp
711 real(DP),
dimension(neqmod),
intent(inout) :: dx
712 integer(I4B),
intent(inout) :: inewtonur
713 real(DP),
intent(inout) :: dxmax
714 integer(I4B),
intent(inout) :: locmax
722 do n = 1, this%dis%nodes
723 if (this%ibound(n) < 1) cycle
724 if (this%icelltype(n) > 0)
then
725 botm = this%dis%bot(this%ibotnode(n))
728 if (x(n) < botm)
then
732 if (abs(dxx) > abs(dxmax))
then
748 real(DP),
intent(inout),
dimension(:) :: hnew
749 real(DP),
intent(inout),
dimension(:) :: flowja
751 integer(I4B) :: n, ipos, m
756 if (this%ixt3d /= 0)
then
757 call this%xt3d%xt3d_flowja(hnew, flowja)
760 do n = 1, this%dis%nodes
761 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
762 m = this%dis%con%ja(ipos)
764 call this%qcalc(n, m, hnew(n), hnew(m), ipos, qnm)
766 flowja(this%dis%con%isym(ipos)) = -qnm
778 integer(I4B),
intent(in) :: n
779 real(DP),
intent(in) :: hn
780 real(DP),
intent(inout) :: thksat
783 if (hn >= this%dis%top(n))
then
786 thksat = (hn - this%dis%bot(n)) / (this%dis%top(n) - this%dis%bot(n))
790 if (this%inewton /= 0)
then
801 integer(I4B),
intent(in) :: n
802 integer(I4B),
intent(in) :: m
803 real(DP),
intent(in) :: hn
804 real(DP),
intent(in) :: hm
805 integer(I4B),
intent(in) :: icon
806 real(DP),
intent(inout) :: qnm
810 real(DP) :: hntemp, hmtemp
814 ihc = this%dis%con%ihc(this%dis%con%jas(icon))
815 hyn = this%hy_eff(n, m, ihc, ipos=icon)
816 hym = this%hy_eff(m, n, ihc, ipos=icon)
820 condnm =
vcond(this%ibound(n), this%ibound(m), &
821 this%icelltype(n), this%icelltype(m), this%inewton, &
822 this%ivarcv, this%idewatcv, &
823 this%condsat(this%dis%con%jas(icon)), hn, hm, &
825 this%sat(n), this%sat(m), &
826 this%dis%top(n), this%dis%top(m), &
827 this%dis%bot(n), this%dis%bot(m), &
828 this%dis%con%hwva(this%dis%con%jas(icon)))
830 condnm =
hcond(this%ibound(n), this%ibound(m), &
831 this%icelltype(n), this%icelltype(m), &
833 this%dis%con%ihc(this%dis%con%jas(icon)), &
835 this%condsat(this%dis%con%jas(icon)), &
836 hn, hm, this%sat(n), this%sat(m), hyn, hym, &
837 this%dis%top(n), this%dis%top(m), &
838 this%dis%bot(n), this%dis%bot(m), &
839 this%dis%con%cl1(this%dis%con%jas(icon)), &
840 this%dis%con%cl2(this%dis%con%jas(icon)), &
841 this%dis%con%hwva(this%dis%con%jas(icon)))
849 if (this%iperched /= 0)
then
850 if (this%dis%con%ihc(this%dis%con%jas(icon)) == 0)
then
852 if (this%icelltype(n) /= 0)
then
853 if (hn < this%dis%top(n)) hntemp = this%dis%bot(m)
856 if (this%icelltype(m) /= 0)
then
857 if (hm < this%dis%top(m)) hmtemp = this%dis%bot(n)
864 qnm = condnm * (hmtemp - hntemp)
872 real(DP),
dimension(:),
intent(in) :: flowja
873 integer(I4B),
intent(in) :: icbcfl
874 integer(I4B),
intent(in) :: icbcun
876 integer(I4B) :: ibinun
879 if (this%ipakcb < 0)
then
881 elseif (this%ipakcb == 0)
then
886 if (icbcfl == 0) ibinun = 0
889 if (ibinun /= 0)
then
890 call this%dis%record_connection_array(flowja, ibinun, this%iout)
894 if (this%isavspdis /= 0)
then
895 if (ibinun /= 0)
call this%sav_spdis(ibinun)
899 if (this%isavsat /= 0)
then
900 if (ibinun /= 0)
call this%sav_sat(ibinun)
912 integer(I4B),
intent(in) :: ibudfl
913 real(DP),
intent(inout),
dimension(:) :: flowja
915 character(len=LENBIGLINE) :: line
916 character(len=30) :: tempstr
917 integer(I4B) :: n, ipos, m
920 character(len=*),
parameter :: fmtiprflow = &
921 &
"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)"
924 if (ibudfl /= 0 .and. this%iprflow > 0)
then
925 write (this%iout, fmtiprflow)
kper,
kstp
926 do n = 1, this%dis%nodes
928 call this%dis%noder_to_string(n, tempstr)
929 line = trim(tempstr)//
':'
930 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
931 m = this%dis%con%ja(ipos)
932 call this%dis%noder_to_string(m, tempstr)
933 line = trim(line)//
' '//trim(tempstr)
935 write (tempstr,
'(1pg15.6)') qnm
936 line = trim(line)//
' '//trim(adjustl(tempstr))
938 write (this%iout,
'(a)') trim(line)
953 if (this%icalcspdis == 1 .and. this%spdis_wa%is_created()) &
954 call this%spdis_wa%destroy()
955 deallocate (this%spdis_wa)
961 if (this%intvk /= 0)
then
963 deallocate (this%tvk)
967 if (this%invsc /= 0)
then
1010 deallocate (this%aname)
1034 call this%NumericalPackageType%da()
1053 call this%NumericalPackageType%allocate_scalars()
1056 call mem_allocate(this%iname,
'INAME', this%memoryPath)
1057 call mem_allocate(this%ixt3d,
'IXT3D', this%memoryPath)
1058 call mem_allocate(this%ixt3drhs,
'IXT3DRHS', this%memoryPath)
1059 call mem_allocate(this%satomega,
'SATOMEGA', this%memoryPath)
1060 call mem_allocate(this%hnoflo,
'HNOFLO', this%memoryPath)
1062 call mem_allocate(this%icellavg,
'ICELLAVG', this%memoryPath)
1063 call mem_allocate(this%iavgkeff,
'IAVGKEFF', this%memoryPath)
1066 call mem_allocate(this%ik22overk,
'IK22OVERK', this%memoryPath)
1067 call mem_allocate(this%ik33overk,
'IK33OVERK', this%memoryPath)
1068 call mem_allocate(this%iperched,
'IPERCHED', this%memoryPath)
1069 call mem_allocate(this%ivarcv,
'IVARCV', this%memoryPath)
1070 call mem_allocate(this%idewatcv,
'IDEWATCV', this%memoryPath)
1071 call mem_allocate(this%ithickstrt,
'ITHICKSTRT', this%memoryPath)
1072 call mem_allocate(this%icalcspdis,
'ICALCSPDIS', this%memoryPath)
1073 call mem_allocate(this%isavspdis,
'ISAVSPDIS', this%memoryPath)
1074 call mem_allocate(this%isavsat,
'ISAVSAT', this%memoryPath)
1075 call mem_allocate(this%irewet,
'IREWET', this%memoryPath)
1076 call mem_allocate(this%wetfct,
'WETFCT', this%memoryPath)
1077 call mem_allocate(this%iwetit,
'IWETIT', this%memoryPath)
1078 call mem_allocate(this%ihdwet,
'IHDWET', this%memoryPath)
1079 call mem_allocate(this%iangle1,
'IANGLE1', this%memoryPath)
1080 call mem_allocate(this%iangle2,
'IANGLE2', this%memoryPath)
1081 call mem_allocate(this%iangle3,
'IANGLE3', this%memoryPath)
1082 call mem_allocate(this%iwetdry,
'IWETDRY', this%memoryPath)
1083 call mem_allocate(this%nedges,
'NEDGES', this%memoryPath)
1084 call mem_allocate(this%lastedge,
'LASTEDGE', this%memoryPath)
1085 call mem_allocate(this%intvk,
'INTVK', this%memoryPath)
1086 call mem_allocate(this%invsc,
'INVSC', this%memoryPath)
1087 call mem_allocate(this%kchangeper,
'KCHANGEPER', this%memoryPath)
1088 call mem_allocate(this%kchangestp,
'KCHANGESTP', this%memoryPath)
1091 call mem_setptr(this%igwfnewtonur,
'INEWTONUR', &
1098 this%satomega =
dzero
1130 this%iasym = this%inewton
1148 integer(I4B),
intent(in) :: ncells
1149 integer(I4B),
intent(in) :: njas
1155 this%k11input(n) = this%k11(n)
1156 this%k22input(n) = this%k22(n)
1157 this%k33input(n) = this%k33(n)
1166 integer(I4B),
intent(in) :: ncells
1167 integer(I4B),
intent(in) :: njas
1171 call mem_allocate(this%ithickstartflag, ncells,
'ITHICKSTARTFLAG', &
1173 call mem_allocate(this%icelltype, ncells,
'ICELLTYPE', this%memoryPath)
1174 call mem_allocate(this%k11, ncells,
'K11', this%memoryPath)
1175 call mem_allocate(this%sat, ncells,
'SAT', this%memoryPath)
1176 call mem_allocate(this%condsat, njas,
'CONDSAT', this%memoryPath)
1179 call mem_allocate(this%k22, ncells,
'K22', this%memoryPath)
1180 call mem_allocate(this%k33, ncells,
'K33', this%memoryPath)
1181 call mem_allocate(this%wetdry, ncells,
'WETDRY', this%memoryPath)
1182 call mem_allocate(this%angle1, ncells,
'ANGLE1', this%memoryPath)
1183 call mem_allocate(this%angle2, ncells,
'ANGLE2', this%memoryPath)
1184 call mem_allocate(this%angle3, ncells,
'ANGLE3', this%memoryPath)
1187 call mem_allocate(this%ibotnode, 0,
'IBOTNODE', this%memoryPath)
1188 call mem_allocate(this%nodedge, 0,
'NODEDGE', this%memoryPath)
1189 call mem_allocate(this%ihcedge, 0,
'IHCEDGE', this%memoryPath)
1190 call mem_allocate(this%propsedge, 0, 0,
'PROPSEDGE', this%memoryPath)
1191 call mem_allocate(this%iedge_ptr, 0,
'NREDGESNODE', this%memoryPath)
1192 call mem_allocate(this%edge_idxs, 0,
'EDGEIDXS', this%memoryPath)
1195 call mem_allocate(this%k11input, 0,
'K11INPUT', this%memoryPath)
1196 call mem_allocate(this%k22input, 0,
'K22INPUT', this%memoryPath)
1197 call mem_allocate(this%k33input, 0,
'K33INPUT', this%memoryPath)
1200 call mem_allocate(this%spdis, 3, 0,
'SPDIS', this%memoryPath)
1203 call mem_allocate(this%nodekchange, ncells,
'NODEKCHANGE', this%memoryPath)
1207 this%angle1(n) =
dzero
1208 this%angle2(n) =
dzero
1209 this%angle3(n) =
dzero
1210 this%wetdry(n) =
dzero
1211 this%nodekchange(n) =
dzero
1215 allocate (this%aname(this%iname))
1216 this%aname = [
' ICELLTYPE',
' K', &
1218 ' WETDRY',
' ANGLE1', &
1219 ' ANGLE2',
' ANGLE3']
1233 write (this%iout,
'(1x,a)')
'Setting NPF Options'
1234 if (found%iprflow) &
1235 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be printed &
1236 &to listing file whenever ICBCFL is not zero.'
1238 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be saved &
1239 &to binary file whenever ICBCFL is not zero.'
1240 if (found%cellavg) &
1241 write (this%iout,
'(4x,a,i0)')
'Alternative cell averaging [1=logarithmic, &
1242 &2=AMT-LMK, 3=AMT-HMK] set to: ', &
1244 if (found%ithickstrt) &
1245 write (this%iout,
'(4x,a)')
'THICKSTRT option has been activated.'
1246 if (found%iperched) &
1247 write (this%iout,
'(4x,a)')
'Vertical flow will be adjusted for perched &
1250 write (this%iout,
'(4x,a)')
'Vertical conductance varies with water table.'
1251 if (found%idewatcv) &
1252 write (this%iout,
'(4x,a)')
'Vertical conductance is calculated using &
1253 &only the saturated thickness and properties &
1254 &of the overlying cell if the head in the &
1255 &underlying cell is below its top.'
1256 if (found%ixt3d)
write (this%iout,
'(4x,a)')
'XT3D formulation is selected.'
1257 if (found%ixt3drhs) &
1258 write (this%iout,
'(4x,a)')
'XT3D RHS formulation is selected.'
1259 if (found%isavspdis) &
1260 write (this%iout,
'(4x,a)')
'Specific discharge will be calculated at cell &
1261 ¢ers and written to DATA-SPDIS in budget &
1262 &file when requested.'
1263 if (found%isavsat) &
1264 write (this%iout,
'(4x,a)')
'Saturation will be written to DATA-SAT in &
1265 &budget file when requested.'
1266 if (found%ik22overk) &
1267 write (this%iout,
'(4x,a)')
'Values specified for K22 are anisotropy &
1268 &ratios and will be multiplied by K before &
1269 &being used in calculations.'
1270 if (found%ik33overk) &
1271 write (this%iout,
'(4x,a)')
'Values specified for K33 are anisotropy &
1272 &ratios and will be multiplied by K before &
1273 &being used in calculations.'
1274 if (found%inewton) &
1275 write (this%iout,
'(4x,a)')
'NEWTON-RAPHSON method disabled for unconfined &
1277 if (found%satomega) &
1278 write (this%iout,
'(4x,a,1pg15.6)')
'Saturation omega: ', this%satomega
1280 write (this%iout,
'(4x,a)')
'Rewetting is active.'
1282 write (this%iout,
'(4x,a,1pg15.6)') &
1283 'Wetting factor (WETFCT) has been set to: ', this%wetfct
1285 write (this%iout,
'(4x,a,i5)') &
1286 'Wetting iteration interval (IWETIT) has been set to: ', this%iwetit
1288 write (this%iout,
'(4x,a,i5)') &
1289 'Head rewet equation (IHDWET) has been set to: ', this%ihdwet
1290 write (this%iout,
'(1x,a,/)')
'End Setting NPF Options'
1306 character(len=LENVARNAME),
dimension(3) :: cellavg_method = &
1307 &[character(len=LENVARNAME) ::
'LOGARITHMIC',
'AMT-LMK',
'AMT-HMK']
1309 character(len=LINELENGTH) :: tvk6_filename
1312 call mem_set_value(this%iprflow,
'IPRFLOW', this%input_mempath, found%iprflow)
1313 call mem_set_value(this%ipakcb,
'IPAKCB', this%input_mempath, found%ipakcb)
1314 call mem_set_value(this%icellavg,
'CELLAVG', this%input_mempath, &
1315 cellavg_method, found%cellavg)
1316 call mem_set_value(this%ithickstrt,
'ITHICKSTRT', this%input_mempath, &
1318 call mem_set_value(this%iperched,
'IPERCHED', this%input_mempath, &
1320 call mem_set_value(this%ivarcv,
'IVARCV', this%input_mempath, found%ivarcv)
1321 call mem_set_value(this%idewatcv,
'IDEWATCV', this%input_mempath, &
1323 call mem_set_value(this%ixt3d,
'IXT3D', this%input_mempath, found%ixt3d)
1324 call mem_set_value(this%ixt3drhs,
'IXT3DRHS', this%input_mempath, &
1326 call mem_set_value(this%isavspdis,
'ISAVSPDIS', this%input_mempath, &
1328 call mem_set_value(this%isavsat,
'ISAVSAT', this%input_mempath, found%isavsat)
1329 call mem_set_value(this%ik22overk,
'IK22OVERK', this%input_mempath, &
1331 call mem_set_value(this%ik33overk,
'IK33OVERK', this%input_mempath, &
1333 call mem_set_value(this%inewton,
'INEWTON', this%input_mempath, found%inewton)
1334 call mem_set_value(this%satomega,
'SATOMEGA', this%input_mempath, &
1336 call mem_set_value(this%irewet,
'IREWET', this%input_mempath, found%irewet)
1337 call mem_set_value(this%wetfct,
'WETFCT', this%input_mempath, found%wetfct)
1338 call mem_set_value(this%iwetit,
'IWETIT', this%input_mempath, found%iwetit)
1339 call mem_set_value(this%ihdwet,
'IHDWET', this%input_mempath, found%ihdwet)
1342 if (found%ipakcb) this%ipakcb = -1
1345 if (found%ixt3d .and. found%ixt3drhs) this%ixt3d = 2
1348 if (found%isavspdis) this%icalcspdis = this%isavspdis
1351 if (found%inewton)
then
1357 if (
filein_fname(tvk6_filename,
'TVK6_FILENAME', this%input_mempath, &
1358 this%input_fname))
then
1359 call openfile(this%intvk, this%iout, tvk6_filename,
'TVK')
1360 call tvk_cr(this%tvk, this%name_model, this%intvk, this%iout)
1364 if (this%iout > 0)
then
1365 call this%log_options(found)
1376 this%icellavg = options%icellavg
1377 this%ithickstrt = options%ithickstrt
1378 this%iperched = options%iperched
1379 this%ivarcv = options%ivarcv
1380 this%idewatcv = options%idewatcv
1381 this%irewet = options%irewet
1382 this%wetfct = options%wetfct
1383 this%iwetit = options%iwetit
1384 this%ihdwet = options%ihdwet
1397 if (this%inewton > 0)
then
1398 this%satomega = dem6
1401 if (this%inewton > 0)
then
1402 if (this%iperched > 0)
then
1403 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1404 'BE USED WITH PERCHED OPTION.'
1407 if (this%ivarcv > 0)
then
1408 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1409 'BE USED WITH VARIABLECV OPTION.'
1412 if (this%irewet > 0)
then
1413 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1414 'BE USED WITH REWET OPTION.'
1419 if (this%ixt3d /= 0)
then
1420 if (this%icellavg > 0)
then
1421 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. '// &
1422 'ALTERNATIVE_CELL_AVERAGING OPTION '// &
1423 'CANNOT BE USED WITH XT3D OPTION.'
1426 if (this%ithickstrt > 0)
then
1427 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. THICKSTRT OPTION '// &
1428 'CANNOT BE USED WITH XT3D OPTION.'
1431 if (this%iperched > 0)
then
1432 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. PERCHED OPTION '// &
1433 'CANNOT BE USED WITH XT3D OPTION.'
1436 if (this%ivarcv > 0)
then
1437 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. VARIABLECV OPTION '// &
1438 'CANNOT BE USED WITH XT3D OPTION.'
1458 write (this%iout,
'(1x,a)')
'Setting NPF Griddata'
1460 if (found%icelltype)
then
1461 write (this%iout,
'(4x,a)')
'ICELLTYPE set from input file'
1465 write (this%iout,
'(4x,a)')
'K set from input file'
1469 write (this%iout,
'(4x,a)')
'K33 set from input file'
1471 write (this%iout,
'(4x,a)')
'K33 not provided. Setting K33 = K.'
1475 write (this%iout,
'(4x,a)')
'K22 set from input file'
1477 write (this%iout,
'(4x,a)')
'K22 not provided. Setting K22 = K.'
1480 if (found%wetdry)
then
1481 write (this%iout,
'(4x,a)')
'WETDRY set from input file'
1484 if (found%angle1)
then
1485 write (this%iout,
'(4x,a)')
'ANGLE1 set from input file'
1488 if (found%angle2)
then
1489 write (this%iout,
'(4x,a)')
'ANGLE2 set from input file'
1492 if (found%angle3)
then
1493 write (this%iout,
'(4x,a)')
'ANGLE3 set from input file'
1496 write (this%iout,
'(1x,a,/)')
'End Setting NPF Griddata'
1510 character(len=LINELENGTH) :: errmsg
1512 logical,
dimension(2) :: afound
1513 integer(I4B),
dimension(:),
pointer,
contiguous :: map
1517 if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
1520 call mem_set_value(this%icelltype,
'ICELLTYPE', this%input_mempath, map, &
1522 call mem_set_value(this%k11,
'K', this%input_mempath, map, found%k)
1523 call mem_set_value(this%k33,
'K33', this%input_mempath, map, found%k33)
1524 call mem_set_value(this%k22,
'K22', this%input_mempath, map, found%k22)
1525 call mem_set_value(this%wetdry,
'WETDRY', this%input_mempath, map, &
1527 call mem_set_value(this%angle1,
'ANGLE1', this%input_mempath, map, &
1529 call mem_set_value(this%angle2,
'ANGLE2', this%input_mempath, map, &
1531 call mem_set_value(this%angle3,
'ANGLE3', this%input_mempath, map, &
1535 if (.not. found%icelltype)
then
1536 write (errmsg,
'(a)')
'Error in GRIDDATA block: ICELLTYPE not found.'
1541 if (.not. found%k)
then
1542 write (errmsg,
'(a)')
'Error in GRIDDATA block: K not found.'
1547 if (.not. found%k33 .and. this%ik33overk /= 0)
then
1548 write (errmsg,
'(a)')
'K33OVERK option specified but K33 not specified.'
1553 if (.not. found%k22 .and. this%ik22overk /= 0)
then
1554 write (errmsg,
'(a)')
'K22OVERK option specified but K22 not specified.'
1559 if (found%k33) this%ik33 = 1
1560 if (found%k22) this%ik22 = 1
1561 if (found%wetdry) this%iwetdry = 1
1562 if (found%angle1) this%iangle1 = 1
1563 if (found%angle2) this%iangle2 = 1
1564 if (found%angle3) this%iangle3 = 1
1567 if (.not. found%k33)
then
1568 call mem_set_value(this%k33,
'K', this%input_mempath, map, afound(1))
1570 if (.not. found%k22)
then
1571 call mem_set_value(this%k22,
'K', this%input_mempath, map, afound(2))
1573 if (.not. found%wetdry)
call mem_reallocate(this%wetdry, 1,
'WETDRY', &
1574 trim(this%memoryPath))
1575 if (.not. found%angle1 .and. this%ixt3d == 0) &
1576 call mem_reallocate(this%angle1, 0,
'ANGLE1', trim(this%memoryPath))
1577 if (.not. found%angle2 .and. this%ixt3d == 0) &
1578 call mem_reallocate(this%angle2, 0,
'ANGLE2', trim(this%memoryPath))
1579 if (.not. found%angle3 .and. this%ixt3d == 0) &
1580 call mem_reallocate(this%angle3, 0,
'ANGLE3', trim(this%memoryPath))
1583 if (this%iout > 0)
then
1584 call this%log_griddata(found)
1597 character(len=24),
dimension(:),
pointer :: aname
1598 character(len=LINELENGTH) :: cellstr, errmsg
1599 integer(I4B) :: nerr, n
1601 character(len=*),
parameter :: fmtkerr = &
1602 &
"(1x, 'Hydraulic property ',a,' is <= 0 for cell ',a, ' ', 1pg15.6)"
1603 character(len=*),
parameter :: fmtkerr2 = &
1604 &
"(1x, '... ', i0,' additional errors not shown for ',a)"
1611 do n = 1,
size(this%k11)
1612 if (this%k11(n) <= dzero)
then
1614 if (nerr <= 20)
then
1615 call this%dis%noder_to_string(n, cellstr)
1616 write (errmsg, fmtkerr) trim(adjustl(aname(2))), trim(cellstr), &
1623 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(2)))
1628 if (this%ik33 /= 0)
then
1632 do n = 1,
size(this%k33)
1633 if (this%ik33overk /= 0) this%k33(n) = this%k33(n) * this%k11(n)
1634 if (this%k33(n) <= dzero)
then
1636 if (nerr <= 20)
then
1637 call this%dis%noder_to_string(n, cellstr)
1638 write (errmsg, fmtkerr) trim(adjustl(aname(3))), trim(cellstr), &
1645 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(3)))
1651 if (this%ik22 /= 0)
then
1654 if (this%dis%con%ianglex == 0)
then
1655 write (errmsg,
'(a)')
'Error. ANGLDEGX not provided in '// &
1656 'discretization file, but K22 was specified. '
1662 do n = 1,
size(this%k22)
1663 if (this%ik22overk /= 0) this%k22(n) = this%k22(n) * this%k11(n)
1664 if (this%k22(n) <= dzero)
then
1666 if (nerr <= 20)
then
1667 call this%dis%noder_to_string(n, cellstr)
1668 write (errmsg, fmtkerr) trim(adjustl(aname(4))), trim(cellstr), &
1675 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(4)))
1681 if (this%irewet == 1)
then
1682 if (this%iwetdry == 0)
then
1683 write (errmsg,
'(a, a, a)')
'Error in GRIDDATA block: ', &
1684 trim(adjustl(aname(5))),
' not found.'
1690 if (this%iangle1 /= 0)
then
1691 do n = 1,
size(this%angle1)
1692 this%angle1(n) = this%angle1(n) *
dpio180
1695 if (this%ixt3d /= 0)
then
1697 write (this%iout,
'(a)')
'XT3D IN USE, BUT ANGLE1 NOT SPECIFIED. '// &
1698 'SETTING ANGLE1 TO ZERO.'
1699 do n = 1,
size(this%angle1)
1700 this%angle1(n) = dzero
1704 if (this%iangle2 /= 0)
then
1705 if (this%iangle1 == 0)
then
1706 write (errmsg,
'(a)')
'ANGLE2 SPECIFIED BUT NOT ANGLE1. '// &
1707 'ANGLE2 REQUIRES ANGLE1. '
1710 if (this%iangle3 == 0)
then
1711 write (errmsg,
'(a)')
'ANGLE2 SPECIFIED BUT NOT ANGLE3. '// &
1712 'SPECIFY BOTH OR NEITHER ONE. '
1715 do n = 1,
size(this%angle2)
1716 this%angle2(n) = this%angle2(n) *
dpio180
1719 if (this%iangle3 /= 0)
then
1720 if (this%iangle1 == 0)
then
1721 write (errmsg,
'(a)')
'ANGLE3 SPECIFIED BUT NOT ANGLE1. '// &
1722 'ANGLE3 REQUIRES ANGLE1. '
1725 if (this%iangle2 == 0)
then
1726 write (errmsg,
'(a)')
'ANGLE3 SPECIFIED BUT NOT ANGLE2. '// &
1727 'SPECIFY BOTH OR NEITHER ONE. '
1730 do n = 1,
size(this%angle3)
1731 this%angle3(n) = this%angle3(n) *
dpio180
1758 integer(I4B) :: n, m, ii, nn
1759 real(DP) :: hyn, hym
1760 real(DP) :: satn, topn, botn
1761 integer(I4B) :: nextn
1762 real(DP) :: minbot, botm
1764 character(len=LINELENGTH) :: cellstr, errmsg
1766 character(len=*),
parameter :: fmtcnv = &
1768 &' ELIMINATED BECAUSE ALL HYDRAULIC CONDUCTIVITIES TO NODE ARE 0.')"
1769 character(len=*),
parameter :: fmtnct = &
1770 &
"(1X,'Negative cell thickness at cell ', A)"
1771 character(len=*),
parameter :: fmtihbe = &
1772 &
"(1X,'Initial head, bottom elevation:',1P,2G13.5)"
1773 character(len=*),
parameter :: fmttebe = &
1774 &
"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
1776 do n = 1, this%dis%nodes
1777 this%ithickstartflag(n) = 0
1783 nodeloop:
do n = 1, this%dis%nodes
1786 if (this%ibound(n) == 0)
then
1787 if (this%irewet /= 0)
then
1788 if (this%wetdry(n) == dzero) cycle nodeloop
1795 if (this%k11(n) /= dzero) cycle nodeloop
1799 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1800 m = this%dis%con%ja(ii)
1801 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0)
then
1803 if (this%ik33 /= 0) hyn = this%k33(n)
1804 if (hyn /= dzero)
then
1806 if (this%ik33 /= 0) hym = this%k33(m)
1807 if (hym /= dzero) cycle
1815 this%hnew(n) = this%hnoflo
1816 if (this%irewet /= 0) this%wetdry(n) = dzero
1817 call this%dis%noder_to_string(n, cellstr)
1818 write (this%iout, fmtcnv) trim(adjustl(cellstr))
1823 if (this%inewton == 0)
then
1826 call this%wd(0, this%hnew)
1832 if (this%ivarcv == 1)
then
1833 do n = 1, this%dis%nodes
1834 if (this%hnew(n) < this%dis%bot(n))
then
1835 this%hnew(n) = this%dis%bot(n) + dem6
1843 if (this%ithickstrt == 0)
then
1844 do n = 1, this%dis%nodes
1845 if (this%icelltype(n) < 0)
then
1846 this%icelltype(n) = 1
1856 do n = 1, this%dis%nodes
1857 if (this%ibound(n) == 0)
then
1859 if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0)
then
1860 this%ithickstartflag(n) = 1
1861 this%icelltype(n) = 0
1864 topn = this%dis%top(n)
1865 botn = this%dis%bot(n)
1866 if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0)
then
1867 call this%thksat(n, this%ic%strt(n), satn)
1868 if (botn > this%ic%strt(n))
then
1869 call this%dis%noder_to_string(n, cellstr)
1870 write (errmsg, fmtnct) trim(adjustl(cellstr))
1872 write (errmsg, fmtihbe) this%ic%strt(n), botn
1875 this%ithickstartflag(n) = 1
1876 this%icelltype(n) = 0
1879 if (botn > topn)
then
1880 call this%dis%noder_to_string(n, cellstr)
1881 write (errmsg, fmtnct) trim(adjustl(cellstr))
1883 write (errmsg, fmttebe) topn, botn
1896 if (this%ixt3d == 0)
then
1901 do n = 1, this%dis%nodes
1902 call this%calc_condsat(n, .true.)
1908 if (this%igwfnewtonur /= 0)
then
1910 trim(this%memoryPath))
1911 do n = 1, this%dis%nodes
1913 minbot = this%dis%bot(n)
1916 do while (.not. finished)
1920 do ii = this%dis%con%ia(nn) + 1, this%dis%con%ia(nn + 1) - 1
1923 m = this%dis%con%ja(ii)
1924 botm = this%dis%bot(m)
1927 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0)
then
1928 if (m > nn .and. botm < minbot)
then
1940 this%ibotnode(n) = nn
1945 this%igwfnewtonur => null()
1956 integer(I4B),
intent(in) :: node
1957 logical,
intent(in) :: upperOnly
1959 integer(I4B) :: ii, m, n, ihc, jj
1960 real(DP) :: topm, topn, topnode, botm, botn, botnode, satm, satn, satnode
1961 real(DP) :: hyn, hym, hn, hm, fawidth, csat
1963 satnode = this%calc_initial_sat(node)
1965 topnode = this%dis%top(node)
1966 botnode = this%dis%bot(node)
1969 do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
1973 m = this%dis%con%ja(ii)
1974 jj = this%dis%con%jas(ii)
1976 if (upperonly) cycle
1983 topn = this%dis%top(n)
1984 botn = this%dis%bot(n)
1985 satn = this%calc_initial_sat(n)
1992 topm = this%dis%top(m)
1993 botm = this%dis%bot(m)
1994 satm = this%calc_initial_sat(m)
1997 ihc = this%dis%con%ihc(jj)
1998 hyn = this%hy_eff(n, m, ihc, ipos=ii)
1999 hym = this%hy_eff(m, n, ihc, ipos=ii)
2000 if (this%ithickstartflag(n) == 0)
then
2003 hn = this%ic%strt(n)
2005 if (this%ithickstartflag(m) == 0)
then
2008 hm = this%ic%strt(m)
2016 csat =
vcond(1, 1, 1, 1, 0, 1, 1,
done, &
2022 this%dis%con%hwva(jj))
2026 fawidth = this%dis%con%hwva(jj)
2027 csat =
hcond(1, 1, 1, 1, 0, &
2031 hn, hm, satn, satm, hyn, hym, &
2034 this%dis%con%cl1(jj), &
2035 this%dis%con%cl2(jj), &
2038 this%condsat(jj) = csat
2052 integer(I4B),
intent(in) :: n
2057 if (this%ibound(n) /= 0 .and. this%ithickstartflag(n) /= 0)
then
2058 call this%thksat(n, this%ic%strt(n), satn)
2071 integer(I4B),
intent(in) :: kiter
2072 real(DP),
intent(inout),
dimension(:) :: hnew
2074 integer(I4B) :: n, m, ii, ihc
2075 real(DP) :: ttop, bbot, thick
2076 integer(I4B) :: ncnvrt, ihdcnv
2077 character(len=30),
dimension(5) :: nodcnvrt
2078 character(len=30) :: nodestr
2079 character(len=3),
dimension(5) :: acnvrt
2080 character(len=LINELENGTH) :: errmsg
2081 integer(I4B) :: irewet
2083 character(len=*),
parameter :: fmtnct = &
2084 "(1X,/1X,'Negative cell thickness at (layer,row,col)', &
2086 character(len=*),
parameter :: fmttopbot = &
2087 &
"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
2088 character(len=*),
parameter :: fmttopbotthk = &
2089 &
"(1X,'Top elevation, bottom elevation, thickness:',1P,3G13.5)"
2090 character(len=*),
parameter :: fmtdrychd = &
2091 &
"(1X,/1X,'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED')"
2092 character(len=*),
parameter :: fmtni = &
2093 &
"(1X,'CELLID=',a,' ITERATION=',I0,' TIME STEP=',I0,' STRESS PERIOD=',I0)"
2100 do n = 1, this%dis%nodes
2101 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2102 m = this%dis%con%ja(ii)
2103 ihc = this%dis%con%ihc(this%dis%con%jas(ii))
2104 call this%rewet_check(kiter, n, hnew(m), this%ibound(m), ihc, hnew, &
2106 if (irewet == 1)
then
2107 call this%wdmsg(2, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2113 do n = 1, this%dis%nodes
2116 if (this%ibound(n) == 0) cycle
2117 if (this%icelltype(n) == 0) cycle
2120 bbot = this%dis%bot(n)
2121 ttop = this%dis%top(n)
2122 if (bbot > ttop)
then
2123 write (errmsg, fmtnct) n
2125 write (errmsg, fmttopbot) ttop, bbot
2131 if (this%icelltype(n) /= 0)
then
2132 if (hnew(n) < ttop) ttop = hnew(n)
2137 if (thick <= dzero)
then
2138 call this%wdmsg(1, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2140 if (this%ibound(n) < 0)
then
2141 errmsg =
'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED'
2143 write (errmsg, fmttopbotthk) ttop, bbot, thick
2145 call this%dis%noder_to_string(n, nodestr)
2146 write (errmsg, fmtni) trim(adjustl(nodestr)), kiter,
kstp,
kper
2155 call this%wdmsg(0, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2158 do n = 1, this%dis%nodes
2159 if (this%ibound(n) == 30000) this%ibound(n) = 1
2170 subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet)
2173 integer(I4B),
intent(in) :: kiter
2174 integer(I4B),
intent(in) :: node
2175 real(DP),
intent(in) :: hm
2176 integer(I4B),
intent(in) :: ibdm
2177 integer(I4B),
intent(in) :: ihc
2178 real(DP),
intent(inout),
dimension(:) :: hnew
2179 integer(I4B),
intent(out) :: irewet
2181 integer(I4B) :: itflg
2182 real(DP) :: wd, awd, turnon, bbot
2187 if (this%irewet > 0)
then
2188 itflg = mod(kiter, this%iwetit)
2189 if (itflg == 0)
then
2190 if (this%ibound(node) == 0 .and. this%wetdry(node) /=
dzero)
then
2193 bbot = this%dis%bot(node)
2194 wd = this%wetdry(node)
2196 if (wd < 0) awd = -wd
2204 if (ibdm > 0 .and. hm >= turnon) irewet = 1
2206 if (wd >
dzero)
then
2209 if (ibdm > 0 .and. hm >= turnon) irewet = 1
2213 if (irewet == 1)
then
2216 if (this%ihdwet == 0)
then
2217 hnew(node) = bbot + this%wetfct * (hm - bbot)
2219 hnew(node) = bbot + this%wetfct * awd
2221 this%ibound(node) = 30000
2236 integer(I4B),
intent(in) :: icode
2237 integer(I4B),
intent(inout) :: ncnvrt
2238 character(len=30),
dimension(5),
intent(inout) :: nodcnvrt
2239 character(len=3),
dimension(5),
intent(inout) :: acnvrt
2240 integer(I4B),
intent(inout) :: ihdcnv
2241 integer(I4B),
intent(in) :: kiter
2242 integer(I4B),
intent(in) :: n
2246 character(len=*),
parameter :: fmtcnvtn = &
2247 "(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I0, &
2248 &' STEP=',I0,' PERIOD=',I0,' (NODE or LRC)')"
2249 character(len=*),
parameter :: fmtnode =
"(1X,3X,5(A4, A20))"
2254 call this%dis%noder_to_string(n, nodcnvrt(ncnvrt))
2255 if (icode == 1)
then
2256 acnvrt(ncnvrt) =
'DRY'
2258 acnvrt(ncnvrt) =
'WET'
2264 if (ncnvrt == 5 .or. (icode == 0 .and. ncnvrt > 0))
then
2265 if (ihdcnv == 0)
write (this%iout, fmtcnvtn) kiter,
kstp,
kper
2267 write (this%iout, fmtnode) &
2268 (acnvrt(l), trim(adjustl(nodcnvrt(l))), l=1, ncnvrt)
2283 function hy_eff(this, n, m, ihc, ipos, vg)
result(hy)
2288 integer(I4B),
intent(in) :: n
2289 integer(I4B),
intent(in) :: m
2290 integer(I4B),
intent(in) :: ihc
2291 integer(I4B),
intent(in),
optional :: ipos
2292 real(dp),
dimension(3),
intent(in),
optional :: vg
2294 integer(I4B) :: iipos
2295 real(dp) :: hy11, hy22, hy33
2296 real(dp) :: ang1, ang2, ang3
2297 real(dp) :: vg1, vg2, vg3
2301 if (
present(ipos)) iipos = ipos
2315 if (this%iangle2 > 0)
then
2316 if (
present(vg))
then
2321 call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2323 ang1 = this%angle1(n)
2324 ang2 = this%angle2(n)
2326 if (this%iangle3 > 0) ang3 = this%angle3(n)
2327 hy =
hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2335 if (this%ik22 > 0)
then
2336 if (
present(vg))
then
2341 call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2346 if (this%iangle1 > 0)
then
2347 ang1 = this%angle1(n)
2348 if (this%iangle2 > 0)
then
2349 ang2 = this%angle2(n)
2350 if (this%iangle3 > 0) ang3 = this%angle3(n)
2353 hy =
hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2367 real(DP),
intent(in),
dimension(:) :: flowja
2371 integer(I4B) :: ipos
2372 integer(I4B) :: iedge
2373 integer(I4B) :: isympos
2401 logical :: nozee = .true.
2405 if (this%icalcspdis /= 0 .and. this%dis%con%ianglex == 0)
then
2406 call store_error(
'Error. ANGLDEGX not provided in '// &
2407 'discretization file. ANGLDEGX required for '// &
2408 'calculation of specific discharge.', terminate=.true.)
2411 swa => this%spdis_wa
2412 if (.not. swa%is_created())
then
2414 call this%spdis_wa%create(this%calc_max_conns())
2417 if (this%nedges > 0)
call this%prepare_edge_lookup()
2421 do n = 1, this%dis%nodes
2431 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2432 m = this%dis%con%ja(ipos)
2433 isympos = this%dis%con%jas(ipos)
2434 ihc = this%dis%con%ihc(isympos)
2435 area = this%dis%con%hwva(isympos)
2441 call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2442 ihc, xc, yc, zc, dltot)
2443 cl1 = this%dis%con%cl1(isympos)
2444 cl2 = this%dis%con%cl2(isympos)
2446 cl1 = this%dis%con%cl2(isympos)
2447 cl2 = this%dis%con%cl1(isympos)
2449 ooclsum =
done / (cl1 + cl2)
2450 swa%diz(iz) = dltot * cl1 * ooclsum
2453 swa%viz(iz) = qz / area
2458 dz =
thksatnm(this%ibound(n), this%ibound(m), &
2459 this%icelltype(n), this%icelltype(m), &
2460 this%inewton, ihc, &
2461 this%hnew(n), this%hnew(m), this%sat(n), this%sat(m), &
2462 this%dis%top(n), this%dis%top(m), this%dis%bot(n), &
2465 call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
2466 call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2467 ihc, xc, yc, zc, dltot)
2468 cl1 = this%dis%con%cl1(isympos)
2469 cl2 = this%dis%con%cl2(isympos)
2471 cl1 = this%dis%con%cl2(isympos)
2472 cl2 = this%dis%con%cl1(isympos)
2474 ooclsum =
done / (cl1 + cl2)
2477 swa%di(ic) = dltot * cl1 * ooclsum
2478 if (area >
dzero)
then
2479 swa%vi(ic) = flowja(ipos) / area
2487 if (this%nedges > 0)
then
2488 do ipos = this%iedge_ptr(n), this%iedge_ptr(n + 1) - 1
2489 iedge = this%edge_idxs(ipos)
2492 ihc = this%ihcedge(iedge)
2493 area = this%propsedge(2, iedge)
2496 swa%viz(iz) = this%propsedge(1, iedge) / area
2497 swa%diz(iz) = this%propsedge(5, iedge)
2500 swa%nix(ic) = -this%propsedge(3, iedge)
2501 swa%niy(ic) = -this%propsedge(4, iedge)
2502 swa%di(ic) = this%propsedge(5, iedge)
2503 if (area >
dzero)
then
2504 swa%vi(ic) = this%propsedge(1, iedge) / area
2522 dsumz = dsumz + swa%diz(iz)
2524 denom = (ncz -
done)
2526 dsumz = dsumz +
dem10 * dsumz
2528 if (dsumz >
dzero) swa%wiz(iz) =
done - swa%diz(iz) / dsumz
2530 swa%wiz(iz) = swa%wiz(iz) / denom
2538 vz = vz + swa%wiz(iz) * swa%viz(iz)
2547 swa%wix(ic) = swa%di(ic) * abs(swa%nix(ic))
2548 swa%wiy(ic) = swa%di(ic) * abs(swa%niy(ic))
2549 dsumx = dsumx + swa%wix(ic)
2550 dsumy = dsumy + swa%wiy(ic)
2557 dsumx = dsumx +
dem10 * dsumx
2558 dsumy = dsumy +
dem10 * dsumy
2560 swa%wix(ic) = (dsumx - swa%wix(ic)) * abs(swa%nix(ic))
2561 swa%wiy(ic) = (dsumy - swa%wiy(ic)) * abs(swa%niy(ic))
2568 swa%bix(ic) = swa%wix(ic) * sign(
done, swa%nix(ic))
2569 swa%biy(ic) = swa%wiy(ic) * sign(
done, swa%niy(ic))
2570 dsumx = dsumx + swa%wix(ic) * abs(swa%nix(ic))
2571 dsumy = dsumy + swa%wiy(ic) * abs(swa%niy(ic))
2578 swa%bix(ic) = swa%bix(ic) * dsumx
2579 swa%biy(ic) = swa%biy(ic) * dsumy
2580 axy = axy + swa%bix(ic) * swa%niy(ic)
2581 ayx = ayx + swa%biy(ic) * swa%nix(ic)
2594 vx = vx + (swa%bix(ic) - axy * swa%biy(ic)) * swa%vi(ic)
2595 vy = vy + (swa%biy(ic) - ayx * swa%bix(ic)) * swa%vi(ic)
2597 denom =
done - axy * ayx
2598 if (denom /=
dzero)
then
2603 this%spdis(1, n) = vx
2604 this%spdis(2, n) = vy
2605 this%spdis(3, n) = vz
2616 integer(I4B),
intent(in) :: ibinun
2618 character(len=16) :: text
2619 character(len=16),
dimension(3) :: auxtxt
2621 integer(I4B) :: naux
2624 text =
' DATA-SPDIS'
2626 auxtxt(:) = [
' qx',
' qy',
' qz']
2627 call this%dis%record_srcdst_list_header(text, this%name_model, &
2628 this%packName, this%name_model, &
2629 this%packName, naux, auxtxt, ibinun, &
2630 this%dis%nodes, this%iout)
2633 do n = 1, this%dis%nodes
2634 call this%dis%record_mf6_list_entry(ibinun, n, n,
dzero, naux, &
2644 integer(I4B),
intent(in) :: ibinun
2646 character(len=16) :: text
2647 character(len=16),
dimension(1) :: auxtxt
2648 real(DP),
dimension(1) :: a
2650 integer(I4B) :: naux
2655 auxtxt(:) = [
' sat']
2656 call this%dis%record_srcdst_list_header(text, this%name_model, &
2657 this%packName, this%name_model, &
2658 this%packName, naux, auxtxt, ibinun, &
2659 this%dis%nodes, this%iout)
2662 do n = 1, this%dis%nodes
2664 call this%dis%record_mf6_list_entry(ibinun, n, n,
dzero, naux, a)
2676 integer(I4B),
intent(in) :: nedges
2678 this%nedges = this%nedges + nedges
2685 integer(I4B) :: max_conns
2687 integer(I4B) :: n, m, ic
2690 do n = 1, this%dis%nodes
2693 ic = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
2696 do m = 1, this%nedges
2697 if (this%nodedge(m) == n)
then
2703 if (ic > max_conns) max_conns = ic
2714 integer(I4B),
intent(in) :: nodedge
2715 integer(I4B),
intent(in) :: ihcedge
2716 real(DP),
intent(in) :: q
2717 real(DP),
intent(in) :: area
2718 real(DP),
intent(in) :: nx
2719 real(DP),
intent(in) :: ny
2720 real(DP),
intent(in) :: distance
2722 integer(I4B) :: lastedge
2724 this%lastedge = this%lastedge + 1
2725 lastedge = this%lastedge
2726 this%nodedge(lastedge) = nodedge
2727 this%ihcedge(lastedge) = ihcedge
2728 this%propsedge(1, lastedge) = q
2729 this%propsedge(2, lastedge) = area
2730 this%propsedge(3, lastedge) = nx
2731 this%propsedge(4, lastedge) = ny
2732 this%propsedge(5, lastedge) = distance
2736 if (this%lastedge == this%nedges) this%lastedge = 0
2742 integer(I4B) :: i, inode, iedge
2743 integer(I4B) :: n, start, end
2744 integer(I4B) :: prev_cnt, strt_idx, ipos
2746 do i = 1,
size(this%iedge_ptr)
2747 this%iedge_ptr(i) = 0
2749 do i = 1,
size(this%edge_idxs)
2750 this%edge_idxs(i) = 0
2754 do iedge = 1, this%nedges
2755 n = this%nodedge(iedge)
2756 this%iedge_ptr(n) = this%iedge_ptr(n) + 1
2760 prev_cnt = this%iedge_ptr(1)
2761 this%iedge_ptr(1) = 1
2762 do inode = 2, this%dis%nodes + 1
2763 strt_idx = this%iedge_ptr(inode - 1) + prev_cnt
2764 prev_cnt = this%iedge_ptr(inode)
2765 this%iedge_ptr(inode) = strt_idx
2769 do iedge = 1, this%nedges
2770 n = this%nodedge(iedge)
2771 start = this%iedge_ptr(n)
2772 end = this%iedge_ptr(n + 1) - 1
2773 do ipos = start,
end
2774 if (this%edge_idxs(ipos) > 0) cycle
2775 this%edge_idxs(ipos) = iedge
2791 real(dp) :: satthickness
2793 satthickness = thksatnm(this%ibound(n), &
2795 this%icelltype(n), &
2796 this%icelltype(m), &
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
@ c3d_vertical
vertical connection
real(dp), parameter dhdry
real dry cell constant
real(dp), parameter dp9
real constant 9/10
real(dp), parameter dem10
real constant 1e-10
real(dp), parameter dem7
real constant 1e-7
real(dp), parameter dem8
real constant 1e-8
integer(i4b), parameter lenbigline
maximum length of a big line
real(dp), parameter dhnoflo
real no flow constant
integer(i4b), parameter lenvarname
maximum length of a variable name
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dpio180
real constant
real(dp), parameter dem6
real constant 1e-6
real(dp), parameter dzero
real constant zero
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 lenmempath
maximum length of the memory path
real(dp), parameter done
real constant 1
This module contains stateless conductance functions.
real(dp) function, public thksatnm(ibdn, ibdm, ictn, ictm, iupstream, ihc, hn, hm, satn, satm, topn, topm, botn, botm)
Calculate wetted cell thickness at interface between two cells.
real(dp) function, public condmean(k1, k2, thick1, thick2, cl1, cl2, width, iavgmeth)
Calculate the conductance between two cells.
real(dp) function, public hcond(ibdn, ibdm, ictn, ictm, iupstream, ihc, icellavg, condsat, hn, hm, satn, satm, hkn, hkm, topn, topm, botn, botm, cln, clm, fawidth)
Horizontal conductance between two cells.
@, public ccond_hmean
Harmonic mean.
real(dp) function, public vcond(ibdn, ibdm, ictn, ictm, inewton, ivarcv, idewatcv, condsat, hn, hm, vkn, vkm, satn, satm, topn, topm, botn, botm, flowarea)
Vertical conductance between two cells.
subroutine set_options(this, options)
Set options in the NPF object.
subroutine npf_save_model_flows(this, flowja, icbcfl, icbcun)
Record flowja and calculate specific discharge if requested.
subroutine source_options(this)
Update simulation options from input mempath.
real(dp) function calc_initial_sat(this, n)
Calculate initial saturation for the given node.
subroutine calc_condsat(this, node, upperOnly)
Calculate CONDSAT array entries for the given node.
subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet)
Determine if a cell should rewet.
integer(i4b) function calc_max_conns(this)
Calculate the maximum number of connections for any cell.
subroutine npf_mc(this, moffset, matrix_sln)
Map connections and construct iax, jax, and idxglox.
subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
Formulate coefficients.
subroutine npf_ac(this, moffset, sparse)
Add connections for extended neighbors to the sparse matrix.
subroutine sgwf_npf_wetdry(this, kiter, hnew)
Perform wetting and drying.
subroutine source_griddata(this)
Update simulation griddata from input mempath.
subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
Under-relaxation.
subroutine sav_spdis(this, ibinun)
Save specific discharge in binary format to ibinun.
subroutine sgwf_npf_thksat(this, n, hn, thksat)
Fractional cell saturation.
subroutine preprocess_input(this)
preprocess the NPF input data
subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, distance)
Provide the npf package with edge properties.
subroutine prepare_edge_lookup(this)
subroutine npf_da(this)
Deallocate variables.
subroutine log_griddata(this, found)
Write dimensions to list file.
real(dp) function hy_eff(this, n, m, ihc, ipos, vg)
Calculate the effective hydraulic conductivity for the n-m connection.
subroutine calc_spdis(this, flowja)
Calculate the 3 components of specific discharge at the cell center.
subroutine, public npf_cr(npfobj, name_model, input_mempath, inunit, iout)
Create a new NPF object. Pass a inunit value of 0 if npf data will initialized from memory.
subroutine increase_edge_count(this, nedges)
Reserve space for nedges cells that have an edge on them.
subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew)
Fill newton terms.
subroutine log_options(this, found)
Log npf options sourced from the input mempath.
subroutine npf_print_model_flows(this, ibudfl, flowja)
Print budget.
subroutine allocate_arrays(this, ncells, njas)
Allocate npf arrays.
subroutine sav_sat(this, ibinun)
Save saturation in binary format to ibinun.
subroutine prepcheck(this)
Initialize and check NPF data.
subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options)
Define the NPF package instance.
subroutine npf_ad(this, nodes, hold, hnew, irestore)
Advance.
subroutine sgwf_npf_wdmsg(this, icode, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
Print wet/dry message.
subroutine sgwf_npf_qcalc(this, n, m, hn, hm, icon, qnm)
Flow between two cells.
subroutine npf_rp(this)
Read and prepare method for package.
real(dp) function calcsatthickness(this, n, m, ihc)
Calculate saturated thickness between cell n and m.
subroutine store_original_k_arrays(this, ncells, njas)
@ brief Store backup copy of hydraulic conductivity when the VSC package is activate
subroutine allocate_scalars(this)
@ brief Allocate scalars
subroutine npf_cq(this, hnew, flowja)
Calculate flowja.
subroutine npf_ar(this, ic, vsc, ibound, hnew)
Allocate and read this NPF instance.
subroutine check_options(this)
Check for conflicting NPF options.
subroutine npf_cf(this, kiter, nodes, hnew)
Routines associated fill coefficients.
General-purpose hydrogeologic functions.
real(dp) function, public hyeff(k11, k22, k33, ang1, ang2, ang3, vg1, vg2, vg3, iavgmeth)
Calculate the effective horizontal hydraulic conductivity from an ellipse using a specified direction...
This module defines variable data types.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
subroutine, public memorystore_remove(component, subcomponent, context)
subroutine, public get_isize(name, mem_path, isize)
@ brief Get the number of elements for this variable
This module contains the base numerical package type.
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
This module contains simulation variables.
character(len=maxcharlen) errmsg
error message string
character(len=linelength) idm_context
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
This module contains the SourceCommonModule.
logical(lgp) function, public filein_fname(filename, tagname, input_mempath, input_fname)
enforce and set a single input filename provided via FILEIN keyword
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number
This module contains time-varying conductivity package methods.
subroutine, public tvk_cr(tvk, name_model, inunit, iout)
Create a new TvkType object.
subroutine, public xt3d_cr(xt3dobj, name_model, inunit, iout, ldispopt)
Create a new xt3d object.
This class is used to store a single deferred-length character string. It was designed to work in an ...
Data structure and helper methods for passing NPF options into npf_df, as an alternative to reading t...
Helper class with work arrays for the SPDIS calculation in NPF.