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)
call this%spdis_wa%destroy()
954 deallocate (this%spdis_wa)
960 if (this%intvk /= 0)
then
962 deallocate (this%tvk)
966 if (this%invsc /= 0)
then
1009 deallocate (this%aname)
1033 call this%NumericalPackageType%da()
1052 call this%NumericalPackageType%allocate_scalars()
1055 call mem_allocate(this%iname,
'INAME', this%memoryPath)
1056 call mem_allocate(this%ixt3d,
'IXT3D', this%memoryPath)
1057 call mem_allocate(this%ixt3drhs,
'IXT3DRHS', this%memoryPath)
1058 call mem_allocate(this%satomega,
'SATOMEGA', this%memoryPath)
1059 call mem_allocate(this%hnoflo,
'HNOFLO', this%memoryPath)
1061 call mem_allocate(this%icellavg,
'ICELLAVG', this%memoryPath)
1062 call mem_allocate(this%iavgkeff,
'IAVGKEFF', this%memoryPath)
1065 call mem_allocate(this%ik22overk,
'IK22OVERK', this%memoryPath)
1066 call mem_allocate(this%ik33overk,
'IK33OVERK', this%memoryPath)
1067 call mem_allocate(this%iperched,
'IPERCHED', this%memoryPath)
1068 call mem_allocate(this%ivarcv,
'IVARCV', this%memoryPath)
1069 call mem_allocate(this%idewatcv,
'IDEWATCV', this%memoryPath)
1070 call mem_allocate(this%ithickstrt,
'ITHICKSTRT', this%memoryPath)
1071 call mem_allocate(this%icalcspdis,
'ICALCSPDIS', this%memoryPath)
1072 call mem_allocate(this%isavspdis,
'ISAVSPDIS', this%memoryPath)
1073 call mem_allocate(this%isavsat,
'ISAVSAT', this%memoryPath)
1074 call mem_allocate(this%irewet,
'IREWET', this%memoryPath)
1075 call mem_allocate(this%wetfct,
'WETFCT', this%memoryPath)
1076 call mem_allocate(this%iwetit,
'IWETIT', this%memoryPath)
1077 call mem_allocate(this%ihdwet,
'IHDWET', this%memoryPath)
1078 call mem_allocate(this%iangle1,
'IANGLE1', this%memoryPath)
1079 call mem_allocate(this%iangle2,
'IANGLE2', this%memoryPath)
1080 call mem_allocate(this%iangle3,
'IANGLE3', this%memoryPath)
1081 call mem_allocate(this%iwetdry,
'IWETDRY', this%memoryPath)
1082 call mem_allocate(this%nedges,
'NEDGES', this%memoryPath)
1083 call mem_allocate(this%lastedge,
'LASTEDGE', this%memoryPath)
1084 call mem_allocate(this%intvk,
'INTVK', this%memoryPath)
1085 call mem_allocate(this%invsc,
'INVSC', this%memoryPath)
1086 call mem_allocate(this%kchangeper,
'KCHANGEPER', this%memoryPath)
1087 call mem_allocate(this%kchangestp,
'KCHANGESTP', this%memoryPath)
1090 call mem_setptr(this%igwfnewtonur,
'INEWTONUR', &
1097 this%satomega =
dzero
1129 this%iasym = this%inewton
1147 integer(I4B),
intent(in) :: ncells
1148 integer(I4B),
intent(in) :: njas
1154 this%k11input(n) = this%k11(n)
1155 this%k22input(n) = this%k22(n)
1156 this%k33input(n) = this%k33(n)
1165 integer(I4B),
intent(in) :: ncells
1166 integer(I4B),
intent(in) :: njas
1170 call mem_allocate(this%ithickstartflag, ncells,
'ITHICKSTARTFLAG', &
1172 call mem_allocate(this%icelltype, ncells,
'ICELLTYPE', this%memoryPath)
1173 call mem_allocate(this%k11, ncells,
'K11', this%memoryPath)
1174 call mem_allocate(this%sat, ncells,
'SAT', this%memoryPath)
1175 call mem_allocate(this%condsat, njas,
'CONDSAT', this%memoryPath)
1178 call mem_allocate(this%k22, ncells,
'K22', this%memoryPath)
1179 call mem_allocate(this%k33, ncells,
'K33', this%memoryPath)
1180 call mem_allocate(this%wetdry, ncells,
'WETDRY', this%memoryPath)
1181 call mem_allocate(this%angle1, ncells,
'ANGLE1', this%memoryPath)
1182 call mem_allocate(this%angle2, ncells,
'ANGLE2', this%memoryPath)
1183 call mem_allocate(this%angle3, ncells,
'ANGLE3', this%memoryPath)
1186 call mem_allocate(this%ibotnode, 0,
'IBOTNODE', this%memoryPath)
1187 call mem_allocate(this%nodedge, 0,
'NODEDGE', this%memoryPath)
1188 call mem_allocate(this%ihcedge, 0,
'IHCEDGE', this%memoryPath)
1189 call mem_allocate(this%propsedge, 0, 0,
'PROPSEDGE', this%memoryPath)
1190 call mem_allocate(this%iedge_ptr, 0,
'NREDGESNODE', this%memoryPath)
1191 call mem_allocate(this%edge_idxs, 0,
'EDGEIDXS', this%memoryPath)
1194 call mem_allocate(this%k11input, 0,
'K11INPUT', this%memoryPath)
1195 call mem_allocate(this%k22input, 0,
'K22INPUT', this%memoryPath)
1196 call mem_allocate(this%k33input, 0,
'K33INPUT', this%memoryPath)
1199 call mem_allocate(this%spdis, 3, 0,
'SPDIS', this%memoryPath)
1202 call mem_allocate(this%nodekchange, ncells,
'NODEKCHANGE', this%memoryPath)
1206 this%angle1(n) =
dzero
1207 this%angle2(n) =
dzero
1208 this%angle3(n) =
dzero
1209 this%wetdry(n) =
dzero
1210 this%nodekchange(n) =
dzero
1214 allocate (this%aname(this%iname))
1215 this%aname = [
' ICELLTYPE',
' K', &
1217 ' WETDRY',
' ANGLE1', &
1218 ' ANGLE2',
' ANGLE3']
1232 write (this%iout,
'(1x,a)')
'Setting NPF Options'
1233 if (found%iprflow) &
1234 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be printed &
1235 &to listing file whenever ICBCFL is not zero.'
1237 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be saved &
1238 &to binary file whenever ICBCFL is not zero.'
1239 if (found%cellavg) &
1240 write (this%iout,
'(4x,a,i0)')
'Alternative cell averaging [1=logarithmic, &
1241 &2=AMT-LMK, 3=AMT-HMK] set to: ', &
1243 if (found%ithickstrt) &
1244 write (this%iout,
'(4x,a)')
'THICKSTRT option has been activated.'
1245 if (found%iperched) &
1246 write (this%iout,
'(4x,a)')
'Vertical flow will be adjusted for perched &
1249 write (this%iout,
'(4x,a)')
'Vertical conductance varies with water table.'
1250 if (found%idewatcv) &
1251 write (this%iout,
'(4x,a)')
'Vertical conductance accounts for dewatered &
1252 &portion of an underlying cell.'
1253 if (found%ixt3d)
write (this%iout,
'(4x,a)')
'XT3D formulation is selected.'
1254 if (found%ixt3drhs) &
1255 write (this%iout,
'(4x,a)')
'XT3D RHS formulation is selected.'
1256 if (found%isavspdis) &
1257 write (this%iout,
'(4x,a)')
'Specific discharge will be calculated at cell &
1258 ¢ers and written to DATA-SPDIS in budget &
1259 &file when requested.'
1260 if (found%isavsat) &
1261 write (this%iout,
'(4x,a)')
'Saturation will be written to DATA-SAT in &
1262 &budget file when requested.'
1263 if (found%ik22overk) &
1264 write (this%iout,
'(4x,a)')
'Values specified for K22 are anisotropy &
1265 &ratios and will be multiplied by K before &
1266 &being used in calculations.'
1267 if (found%ik33overk) &
1268 write (this%iout,
'(4x,a)')
'Values specified for K33 are anisotropy &
1269 &ratios and will be multiplied by K before &
1270 &being used in calculations.'
1271 if (found%inewton) &
1272 write (this%iout,
'(4x,a)')
'NEWTON-RAPHSON method disabled for unconfined &
1274 if (found%satomega) &
1275 write (this%iout,
'(4x,a,1pg15.6)')
'Saturation omega: ', this%satomega
1277 write (this%iout,
'(4x,a)')
'Rewetting is active.'
1279 write (this%iout,
'(4x,a,1pg15.6)') &
1280 'Wetting factor (WETFCT) has been set to: ', this%wetfct
1282 write (this%iout,
'(4x,a,i5)') &
1283 'Wetting iteration interval (IWETIT) has been set to: ', this%iwetit
1285 write (this%iout,
'(4x,a,i5)') &
1286 'Head rewet equation (IHDWET) has been set to: ', this%ihdwet
1287 write (this%iout,
'(1x,a,/)')
'End Setting NPF Options'
1303 character(len=LENVARNAME),
dimension(3) :: cellavg_method = &
1304 &[character(len=LENVARNAME) ::
'LOGARITHMIC',
'AMT-LMK',
'AMT-HMK']
1306 character(len=LINELENGTH) :: tvk6_filename
1309 call mem_set_value(this%iprflow,
'IPRFLOW', this%input_mempath, found%iprflow)
1310 call mem_set_value(this%ipakcb,
'IPAKCB', this%input_mempath, found%ipakcb)
1311 call mem_set_value(this%icellavg,
'CELLAVG', this%input_mempath, &
1312 cellavg_method, found%cellavg)
1313 call mem_set_value(this%ithickstrt,
'ITHICKSTRT', this%input_mempath, &
1315 call mem_set_value(this%iperched,
'IPERCHED', this%input_mempath, &
1317 call mem_set_value(this%ivarcv,
'IVARCV', this%input_mempath, found%ivarcv)
1318 call mem_set_value(this%idewatcv,
'IDEWATCV', this%input_mempath, &
1320 call mem_set_value(this%ixt3d,
'IXT3D', this%input_mempath, found%ixt3d)
1321 call mem_set_value(this%ixt3drhs,
'IXT3DRHS', this%input_mempath, &
1323 call mem_set_value(this%isavspdis,
'ISAVSPDIS', this%input_mempath, &
1325 call mem_set_value(this%isavsat,
'ISAVSAT', this%input_mempath, found%isavsat)
1326 call mem_set_value(this%ik22overk,
'IK22OVERK', this%input_mempath, &
1328 call mem_set_value(this%ik33overk,
'IK33OVERK', this%input_mempath, &
1330 call mem_set_value(this%inewton,
'INEWTON', this%input_mempath, found%inewton)
1331 call mem_set_value(this%satomega,
'SATOMEGA', this%input_mempath, &
1333 call mem_set_value(this%irewet,
'IREWET', this%input_mempath, found%irewet)
1334 call mem_set_value(this%wetfct,
'WETFCT', this%input_mempath, found%wetfct)
1335 call mem_set_value(this%iwetit,
'IWETIT', this%input_mempath, found%iwetit)
1336 call mem_set_value(this%ihdwet,
'IHDWET', this%input_mempath, found%ihdwet)
1339 if (found%ipakcb) this%ipakcb = -1
1342 if (found%ixt3d .and. found%ixt3drhs) this%ixt3d = 2
1345 if (found%isavspdis) this%icalcspdis = this%isavspdis
1348 if (found%inewton)
then
1354 if (
filein_fname(tvk6_filename,
'TVK6_FILENAME', this%input_mempath, &
1355 this%input_fname))
then
1356 call openfile(this%intvk, this%iout, tvk6_filename,
'TVK')
1357 call tvk_cr(this%tvk, this%name_model, this%intvk, this%iout)
1361 if (this%iout > 0)
then
1362 call this%log_options(found)
1373 this%icellavg = options%icellavg
1374 this%ithickstrt = options%ithickstrt
1375 this%iperched = options%iperched
1376 this%ivarcv = options%ivarcv
1377 this%idewatcv = options%idewatcv
1378 this%irewet = options%irewet
1379 this%wetfct = options%wetfct
1380 this%iwetit = options%iwetit
1381 this%ihdwet = options%ihdwet
1394 if (this%inewton > 0)
then
1395 this%satomega = dem6
1398 if (this%inewton > 0)
then
1399 if (this%iperched > 0)
then
1400 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1401 'BE USED WITH PERCHED OPTION.'
1404 if (this%ivarcv > 0)
then
1405 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1406 'BE USED WITH VARIABLECV OPTION.'
1409 if (this%irewet > 0)
then
1410 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1411 'BE USED WITH REWET OPTION.'
1416 if (this%ixt3d /= 0)
then
1417 if (this%icellavg > 0)
then
1418 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. '// &
1419 'ALTERNATIVE_CELL_AVERAGING OPTION '// &
1420 'CANNOT BE USED WITH XT3D OPTION.'
1423 if (this%ithickstrt > 0)
then
1424 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. THICKSTRT OPTION '// &
1425 'CANNOT BE USED WITH XT3D OPTION.'
1428 if (this%iperched > 0)
then
1429 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. PERCHED OPTION '// &
1430 'CANNOT BE USED WITH XT3D OPTION.'
1433 if (this%ivarcv > 0)
then
1434 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. VARIABLECV OPTION '// &
1435 'CANNOT BE USED WITH XT3D OPTION.'
1455 write (this%iout,
'(1x,a)')
'Setting NPF Griddata'
1457 if (found%icelltype)
then
1458 write (this%iout,
'(4x,a)')
'ICELLTYPE set from input file'
1462 write (this%iout,
'(4x,a)')
'K set from input file'
1466 write (this%iout,
'(4x,a)')
'K33 set from input file'
1468 write (this%iout,
'(4x,a)')
'K33 not provided. Setting K33 = K.'
1472 write (this%iout,
'(4x,a)')
'K22 set from input file'
1474 write (this%iout,
'(4x,a)')
'K22 not provided. Setting K22 = K.'
1477 if (found%wetdry)
then
1478 write (this%iout,
'(4x,a)')
'WETDRY set from input file'
1481 if (found%angle1)
then
1482 write (this%iout,
'(4x,a)')
'ANGLE1 set from input file'
1485 if (found%angle2)
then
1486 write (this%iout,
'(4x,a)')
'ANGLE2 set from input file'
1489 if (found%angle3)
then
1490 write (this%iout,
'(4x,a)')
'ANGLE3 set from input file'
1493 write (this%iout,
'(1x,a,/)')
'End Setting NPF Griddata'
1507 character(len=LINELENGTH) :: errmsg
1509 logical,
dimension(2) :: afound
1510 integer(I4B),
dimension(:),
pointer,
contiguous :: map
1514 if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
1517 call mem_set_value(this%icelltype,
'ICELLTYPE', this%input_mempath, map, &
1519 call mem_set_value(this%k11,
'K', this%input_mempath, map, found%k)
1520 call mem_set_value(this%k33,
'K33', this%input_mempath, map, found%k33)
1521 call mem_set_value(this%k22,
'K22', this%input_mempath, map, found%k22)
1522 call mem_set_value(this%wetdry,
'WETDRY', this%input_mempath, map, &
1524 call mem_set_value(this%angle1,
'ANGLE1', this%input_mempath, map, &
1526 call mem_set_value(this%angle2,
'ANGLE2', this%input_mempath, map, &
1528 call mem_set_value(this%angle3,
'ANGLE3', this%input_mempath, map, &
1532 if (.not. found%icelltype)
then
1533 write (errmsg,
'(a)')
'Error in GRIDDATA block: ICELLTYPE not found.'
1538 if (.not. found%k)
then
1539 write (errmsg,
'(a)')
'Error in GRIDDATA block: K not found.'
1544 if (.not. found%k33 .and. this%ik33overk /= 0)
then
1545 write (errmsg,
'(a)')
'K33OVERK option specified but K33 not specified.'
1550 if (.not. found%k22 .and. this%ik22overk /= 0)
then
1551 write (errmsg,
'(a)')
'K22OVERK option specified but K22 not specified.'
1556 if (found%k33) this%ik33 = 1
1557 if (found%k22) this%ik22 = 1
1558 if (found%wetdry) this%iwetdry = 1
1559 if (found%angle1) this%iangle1 = 1
1560 if (found%angle2) this%iangle2 = 1
1561 if (found%angle3) this%iangle3 = 1
1564 if (.not. found%k33)
then
1565 call mem_set_value(this%k33,
'K', this%input_mempath, map, afound(1))
1567 if (.not. found%k22)
then
1568 call mem_set_value(this%k22,
'K', this%input_mempath, map, afound(2))
1570 if (.not. found%wetdry)
call mem_reallocate(this%wetdry, 1,
'WETDRY', &
1571 trim(this%memoryPath))
1572 if (.not. found%angle1 .and. this%ixt3d == 0) &
1573 call mem_reallocate(this%angle1, 0,
'ANGLE1', trim(this%memoryPath))
1574 if (.not. found%angle2 .and. this%ixt3d == 0) &
1575 call mem_reallocate(this%angle2, 0,
'ANGLE2', trim(this%memoryPath))
1576 if (.not. found%angle3 .and. this%ixt3d == 0) &
1577 call mem_reallocate(this%angle3, 0,
'ANGLE3', trim(this%memoryPath))
1580 if (this%iout > 0)
then
1581 call this%log_griddata(found)
1594 character(len=24),
dimension(:),
pointer :: aname
1595 character(len=LINELENGTH) :: cellstr, errmsg
1596 integer(I4B) :: nerr, n
1598 character(len=*),
parameter :: fmtkerr = &
1599 &
"(1x, 'Hydraulic property ',a,' is <= 0 for cell ',a, ' ', 1pg15.6)"
1600 character(len=*),
parameter :: fmtkerr2 = &
1601 &
"(1x, '... ', i0,' additional errors not shown for ',a)"
1608 do n = 1,
size(this%k11)
1609 if (this%k11(n) <= dzero)
then
1611 if (nerr <= 20)
then
1612 call this%dis%noder_to_string(n, cellstr)
1613 write (errmsg, fmtkerr) trim(adjustl(aname(2))), trim(cellstr), &
1620 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(2)))
1625 if (this%ik33 /= 0)
then
1629 do n = 1,
size(this%k33)
1630 if (this%ik33overk /= 0) this%k33(n) = this%k33(n) * this%k11(n)
1631 if (this%k33(n) <= dzero)
then
1633 if (nerr <= 20)
then
1634 call this%dis%noder_to_string(n, cellstr)
1635 write (errmsg, fmtkerr) trim(adjustl(aname(3))), trim(cellstr), &
1642 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(3)))
1648 if (this%ik22 /= 0)
then
1651 if (this%dis%con%ianglex == 0)
then
1652 write (errmsg,
'(a)')
'Error. ANGLDEGX not provided in '// &
1653 'discretization file, but K22 was specified. '
1659 do n = 1,
size(this%k22)
1660 if (this%ik22overk /= 0) this%k22(n) = this%k22(n) * this%k11(n)
1661 if (this%k22(n) <= dzero)
then
1663 if (nerr <= 20)
then
1664 call this%dis%noder_to_string(n, cellstr)
1665 write (errmsg, fmtkerr) trim(adjustl(aname(4))), trim(cellstr), &
1672 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(4)))
1678 if (this%irewet == 1)
then
1679 if (this%iwetdry == 0)
then
1680 write (errmsg,
'(a, a, a)')
'Error in GRIDDATA block: ', &
1681 trim(adjustl(aname(5))),
' not found.'
1687 if (this%iangle1 /= 0)
then
1688 do n = 1,
size(this%angle1)
1689 this%angle1(n) = this%angle1(n) *
dpio180
1692 if (this%ixt3d /= 0)
then
1694 write (this%iout,
'(a)')
'XT3D IN USE, BUT ANGLE1 NOT SPECIFIED. '// &
1695 'SETTING ANGLE1 TO ZERO.'
1696 do n = 1,
size(this%angle1)
1697 this%angle1(n) = dzero
1701 if (this%iangle2 /= 0)
then
1702 if (this%iangle1 == 0)
then
1703 write (errmsg,
'(a)')
'ANGLE2 SPECIFIED BUT NOT ANGLE1. '// &
1704 'ANGLE2 REQUIRES ANGLE1. '
1707 if (this%iangle3 == 0)
then
1708 write (errmsg,
'(a)')
'ANGLE2 SPECIFIED BUT NOT ANGLE3. '// &
1709 'SPECIFY BOTH OR NEITHER ONE. '
1712 do n = 1,
size(this%angle2)
1713 this%angle2(n) = this%angle2(n) *
dpio180
1716 if (this%iangle3 /= 0)
then
1717 if (this%iangle1 == 0)
then
1718 write (errmsg,
'(a)')
'ANGLE3 SPECIFIED BUT NOT ANGLE1. '// &
1719 'ANGLE3 REQUIRES ANGLE1. '
1722 if (this%iangle2 == 0)
then
1723 write (errmsg,
'(a)')
'ANGLE3 SPECIFIED BUT NOT ANGLE2. '// &
1724 'SPECIFY BOTH OR NEITHER ONE. '
1727 do n = 1,
size(this%angle3)
1728 this%angle3(n) = this%angle3(n) *
dpio180
1755 integer(I4B) :: n, m, ii, nn
1756 real(DP) :: hyn, hym
1757 real(DP) :: satn, topn, botn
1758 integer(I4B) :: nextn
1759 real(DP) :: minbot, botm
1761 character(len=LINELENGTH) :: cellstr, errmsg
1763 character(len=*),
parameter :: fmtcnv = &
1765 &' ELIMINATED BECAUSE ALL HYDRAULIC CONDUCTIVITIES TO NODE ARE 0.')"
1766 character(len=*),
parameter :: fmtnct = &
1767 &
"(1X,'Negative cell thickness at cell ', A)"
1768 character(len=*),
parameter :: fmtihbe = &
1769 &
"(1X,'Initial head, bottom elevation:',1P,2G13.5)"
1770 character(len=*),
parameter :: fmttebe = &
1771 &
"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
1773 do n = 1, this%dis%nodes
1774 this%ithickstartflag(n) = 0
1780 nodeloop:
do n = 1, this%dis%nodes
1783 if (this%ibound(n) == 0)
then
1784 if (this%irewet /= 0)
then
1785 if (this%wetdry(n) == dzero) cycle nodeloop
1792 if (this%k11(n) /= dzero) cycle nodeloop
1796 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1797 m = this%dis%con%ja(ii)
1798 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0)
then
1800 if (this%ik33 /= 0) hyn = this%k33(n)
1801 if (hyn /= dzero)
then
1803 if (this%ik33 /= 0) hym = this%k33(m)
1804 if (hym /= dzero) cycle
1812 this%hnew(n) = this%hnoflo
1813 if (this%irewet /= 0) this%wetdry(n) = dzero
1814 call this%dis%noder_to_string(n, cellstr)
1815 write (this%iout, fmtcnv) trim(adjustl(cellstr))
1820 if (this%inewton == 0)
then
1823 call this%wd(0, this%hnew)
1829 if (this%ivarcv == 1)
then
1830 do n = 1, this%dis%nodes
1831 if (this%hnew(n) < this%dis%bot(n))
then
1832 this%hnew(n) = this%dis%bot(n) + dem6
1840 if (this%ithickstrt == 0)
then
1841 do n = 1, this%dis%nodes
1842 if (this%icelltype(n) < 0)
then
1843 this%icelltype(n) = 1
1853 do n = 1, this%dis%nodes
1854 if (this%ibound(n) == 0)
then
1856 if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0)
then
1857 this%ithickstartflag(n) = 1
1858 this%icelltype(n) = 0
1861 topn = this%dis%top(n)
1862 botn = this%dis%bot(n)
1863 if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0)
then
1864 call this%thksat(n, this%ic%strt(n), satn)
1865 if (botn > this%ic%strt(n))
then
1866 call this%dis%noder_to_string(n, cellstr)
1867 write (errmsg, fmtnct) trim(adjustl(cellstr))
1869 write (errmsg, fmtihbe) this%ic%strt(n), botn
1872 this%ithickstartflag(n) = 1
1873 this%icelltype(n) = 0
1876 if (botn > topn)
then
1877 call this%dis%noder_to_string(n, cellstr)
1878 write (errmsg, fmtnct) trim(adjustl(cellstr))
1880 write (errmsg, fmttebe) topn, botn
1893 if (this%ixt3d == 0)
then
1898 do n = 1, this%dis%nodes
1899 call this%calc_condsat(n, .true.)
1905 if (this%igwfnewtonur /= 0)
then
1907 trim(this%memoryPath))
1908 do n = 1, this%dis%nodes
1910 minbot = this%dis%bot(n)
1913 do while (.not. finished)
1917 do ii = this%dis%con%ia(nn) + 1, this%dis%con%ia(nn + 1) - 1
1920 m = this%dis%con%ja(ii)
1921 botm = this%dis%bot(m)
1924 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0)
then
1925 if (m > nn .and. botm < minbot)
then
1937 this%ibotnode(n) = nn
1942 this%igwfnewtonur => null()
1953 integer(I4B),
intent(in) :: node
1954 logical,
intent(in) :: upperOnly
1956 integer(I4B) :: ii, m, n, ihc, jj
1957 real(DP) :: topm, topn, topnode, botm, botn, botnode, satm, satn, satnode
1958 real(DP) :: hyn, hym, hn, hm, fawidth, csat
1960 satnode = this%calc_initial_sat(node)
1962 topnode = this%dis%top(node)
1963 botnode = this%dis%bot(node)
1966 do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
1970 m = this%dis%con%ja(ii)
1971 jj = this%dis%con%jas(ii)
1973 if (upperonly) cycle
1980 topn = this%dis%top(n)
1981 botn = this%dis%bot(n)
1982 satn = this%calc_initial_sat(n)
1989 topm = this%dis%top(m)
1990 botm = this%dis%bot(m)
1991 satm = this%calc_initial_sat(m)
1994 ihc = this%dis%con%ihc(jj)
1995 hyn = this%hy_eff(n, m, ihc, ipos=ii)
1996 hym = this%hy_eff(m, n, ihc, ipos=ii)
1997 if (this%ithickstartflag(n) == 0)
then
2000 hn = this%ic%strt(n)
2002 if (this%ithickstartflag(m) == 0)
then
2005 hm = this%ic%strt(m)
2013 csat =
vcond(1, 1, 1, 1, 0, 1, 1,
done, &
2019 this%dis%con%hwva(jj))
2023 fawidth = this%dis%con%hwva(jj)
2024 csat =
hcond(1, 1, 1, 1, 0, &
2028 hn, hm, satn, satm, hyn, hym, &
2031 this%dis%con%cl1(jj), &
2032 this%dis%con%cl2(jj), &
2035 this%condsat(jj) = csat
2049 integer(I4B),
intent(in) :: n
2054 if (this%ibound(n) /= 0 .and. this%ithickstartflag(n) /= 0)
then
2055 call this%thksat(n, this%ic%strt(n), satn)
2068 integer(I4B),
intent(in) :: kiter
2069 real(DP),
intent(inout),
dimension(:) :: hnew
2071 integer(I4B) :: n, m, ii, ihc
2072 real(DP) :: ttop, bbot, thick
2073 integer(I4B) :: ncnvrt, ihdcnv
2074 character(len=30),
dimension(5) :: nodcnvrt
2075 character(len=30) :: nodestr
2076 character(len=3),
dimension(5) :: acnvrt
2077 character(len=LINELENGTH) :: errmsg
2078 integer(I4B) :: irewet
2080 character(len=*),
parameter :: fmtnct = &
2081 "(1X,/1X,'Negative cell thickness at (layer,row,col)', &
2083 character(len=*),
parameter :: fmttopbot = &
2084 &
"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
2085 character(len=*),
parameter :: fmttopbotthk = &
2086 &
"(1X,'Top elevation, bottom elevation, thickness:',1P,3G13.5)"
2087 character(len=*),
parameter :: fmtdrychd = &
2088 &
"(1X,/1X,'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED')"
2089 character(len=*),
parameter :: fmtni = &
2090 &
"(1X,'CELLID=',a,' ITERATION=',I0,' TIME STEP=',I0,' STRESS PERIOD=',I0)"
2097 do n = 1, this%dis%nodes
2098 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2099 m = this%dis%con%ja(ii)
2100 ihc = this%dis%con%ihc(this%dis%con%jas(ii))
2101 call this%rewet_check(kiter, n, hnew(m), this%ibound(m), ihc, hnew, &
2103 if (irewet == 1)
then
2104 call this%wdmsg(2, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2110 do n = 1, this%dis%nodes
2113 if (this%ibound(n) == 0) cycle
2114 if (this%icelltype(n) == 0) cycle
2117 bbot = this%dis%bot(n)
2118 ttop = this%dis%top(n)
2119 if (bbot > ttop)
then
2120 write (errmsg, fmtnct) n
2122 write (errmsg, fmttopbot) ttop, bbot
2128 if (this%icelltype(n) /= 0)
then
2129 if (hnew(n) < ttop) ttop = hnew(n)
2134 if (thick <= dzero)
then
2135 call this%wdmsg(1, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2137 if (this%ibound(n) < 0)
then
2138 errmsg =
'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED'
2140 write (errmsg, fmttopbotthk) ttop, bbot, thick
2142 call this%dis%noder_to_string(n, nodestr)
2143 write (errmsg, fmtni) trim(adjustl(nodestr)), kiter,
kstp,
kper
2152 call this%wdmsg(0, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2155 do n = 1, this%dis%nodes
2156 if (this%ibound(n) == 30000) this%ibound(n) = 1
2167 subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet)
2170 integer(I4B),
intent(in) :: kiter
2171 integer(I4B),
intent(in) :: node
2172 real(DP),
intent(in) :: hm
2173 integer(I4B),
intent(in) :: ibdm
2174 integer(I4B),
intent(in) :: ihc
2175 real(DP),
intent(inout),
dimension(:) :: hnew
2176 integer(I4B),
intent(out) :: irewet
2178 integer(I4B) :: itflg
2179 real(DP) :: wd, awd, turnon, bbot
2184 if (this%irewet > 0)
then
2185 itflg = mod(kiter, this%iwetit)
2186 if (itflg == 0)
then
2187 if (this%ibound(node) == 0 .and. this%wetdry(node) /=
dzero)
then
2190 bbot = this%dis%bot(node)
2191 wd = this%wetdry(node)
2193 if (wd < 0) awd = -wd
2201 if (ibdm > 0 .and. hm >= turnon) irewet = 1
2203 if (wd >
dzero)
then
2206 if (ibdm > 0 .and. hm >= turnon) irewet = 1
2210 if (irewet == 1)
then
2213 if (this%ihdwet == 0)
then
2214 hnew(node) = bbot + this%wetfct * (hm - bbot)
2216 hnew(node) = bbot + this%wetfct * awd
2218 this%ibound(node) = 30000
2233 integer(I4B),
intent(in) :: icode
2234 integer(I4B),
intent(inout) :: ncnvrt
2235 character(len=30),
dimension(5),
intent(inout) :: nodcnvrt
2236 character(len=3),
dimension(5),
intent(inout) :: acnvrt
2237 integer(I4B),
intent(inout) :: ihdcnv
2238 integer(I4B),
intent(in) :: kiter
2239 integer(I4B),
intent(in) :: n
2243 character(len=*),
parameter :: fmtcnvtn = &
2244 "(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I0, &
2245 &' STEP=',I0,' PERIOD=',I0,' (NODE or LRC)')"
2246 character(len=*),
parameter :: fmtnode =
"(1X,3X,5(A4, A20))"
2251 call this%dis%noder_to_string(n, nodcnvrt(ncnvrt))
2252 if (icode == 1)
then
2253 acnvrt(ncnvrt) =
'DRY'
2255 acnvrt(ncnvrt) =
'WET'
2261 if (ncnvrt == 5 .or. (icode == 0 .and. ncnvrt > 0))
then
2262 if (ihdcnv == 0)
write (this%iout, fmtcnvtn) kiter,
kstp,
kper
2264 write (this%iout, fmtnode) &
2265 (acnvrt(l), trim(adjustl(nodcnvrt(l))), l=1, ncnvrt)
2280 function hy_eff(this, n, m, ihc, ipos, vg)
result(hy)
2285 integer(I4B),
intent(in) :: n
2286 integer(I4B),
intent(in) :: m
2287 integer(I4B),
intent(in) :: ihc
2288 integer(I4B),
intent(in),
optional :: ipos
2289 real(dp),
dimension(3),
intent(in),
optional :: vg
2291 integer(I4B) :: iipos
2292 real(dp) :: hy11, hy22, hy33
2293 real(dp) :: ang1, ang2, ang3
2294 real(dp) :: vg1, vg2, vg3
2298 if (
present(ipos)) iipos = ipos
2312 if (this%iangle2 > 0)
then
2313 if (
present(vg))
then
2318 call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2320 ang1 = this%angle1(n)
2321 ang2 = this%angle2(n)
2323 if (this%iangle3 > 0) ang3 = this%angle3(n)
2324 hy =
hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2332 if (this%ik22 > 0)
then
2333 if (
present(vg))
then
2338 call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2343 if (this%iangle1 > 0)
then
2344 ang1 = this%angle1(n)
2345 if (this%iangle2 > 0)
then
2346 ang2 = this%angle2(n)
2347 if (this%iangle3 > 0) ang3 = this%angle3(n)
2350 hy =
hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2364 real(DP),
intent(in),
dimension(:) :: flowja
2368 integer(I4B) :: ipos
2369 integer(I4B) :: iedge
2370 integer(I4B) :: isympos
2398 logical :: nozee = .true.
2402 if (this%icalcspdis /= 0 .and. this%dis%con%ianglex == 0)
then
2403 call store_error(
'Error. ANGLDEGX not provided in '// &
2404 'discretization file. ANGLDEGX required for '// &
2405 'calculation of specific discharge.', terminate=.true.)
2408 swa => this%spdis_wa
2409 if (.not. swa%is_created())
then
2411 call this%spdis_wa%create(this%calc_max_conns())
2414 if (this%nedges > 0)
call this%prepare_edge_lookup()
2418 do n = 1, this%dis%nodes
2428 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2429 m = this%dis%con%ja(ipos)
2430 isympos = this%dis%con%jas(ipos)
2431 ihc = this%dis%con%ihc(isympos)
2432 area = this%dis%con%hwva(isympos)
2438 call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2439 ihc, xc, yc, zc, dltot)
2440 cl1 = this%dis%con%cl1(isympos)
2441 cl2 = this%dis%con%cl2(isympos)
2443 cl1 = this%dis%con%cl2(isympos)
2444 cl2 = this%dis%con%cl1(isympos)
2446 ooclsum =
done / (cl1 + cl2)
2447 swa%diz(iz) = dltot * cl1 * ooclsum
2450 swa%viz(iz) = qz / area
2455 dz =
thksatnm(this%ibound(n), this%ibound(m), &
2456 this%icelltype(n), this%icelltype(m), &
2457 this%inewton, ihc, &
2458 this%hnew(n), this%hnew(m), this%sat(n), this%sat(m), &
2459 this%dis%top(n), this%dis%top(m), this%dis%bot(n), &
2462 call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
2463 call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2464 ihc, xc, yc, zc, dltot)
2465 cl1 = this%dis%con%cl1(isympos)
2466 cl2 = this%dis%con%cl2(isympos)
2468 cl1 = this%dis%con%cl2(isympos)
2469 cl2 = this%dis%con%cl1(isympos)
2471 ooclsum =
done / (cl1 + cl2)
2474 swa%di(ic) = dltot * cl1 * ooclsum
2475 if (area >
dzero)
then
2476 swa%vi(ic) = flowja(ipos) / area
2484 if (this%nedges > 0)
then
2485 do ipos = this%iedge_ptr(n), this%iedge_ptr(n + 1) - 1
2486 iedge = this%edge_idxs(ipos)
2489 ihc = this%ihcedge(iedge)
2490 area = this%propsedge(2, iedge)
2493 swa%viz(iz) = this%propsedge(1, iedge) / area
2494 swa%diz(iz) = this%propsedge(5, iedge)
2497 swa%nix(ic) = -this%propsedge(3, iedge)
2498 swa%niy(ic) = -this%propsedge(4, iedge)
2499 swa%di(ic) = this%propsedge(5, iedge)
2500 if (area >
dzero)
then
2501 swa%vi(ic) = this%propsedge(1, iedge) / area
2519 dsumz = dsumz + swa%diz(iz)
2521 denom = (ncz -
done)
2523 dsumz = dsumz +
dem10 * dsumz
2525 if (dsumz >
dzero) swa%wiz(iz) =
done - swa%diz(iz) / dsumz
2527 swa%wiz(iz) = swa%wiz(iz) / denom
2535 vz = vz + swa%wiz(iz) * swa%viz(iz)
2544 swa%wix(ic) = swa%di(ic) * abs(swa%nix(ic))
2545 swa%wiy(ic) = swa%di(ic) * abs(swa%niy(ic))
2546 dsumx = dsumx + swa%wix(ic)
2547 dsumy = dsumy + swa%wiy(ic)
2554 dsumx = dsumx +
dem10 * dsumx
2555 dsumy = dsumy +
dem10 * dsumy
2557 swa%wix(ic) = (dsumx - swa%wix(ic)) * abs(swa%nix(ic))
2558 swa%wiy(ic) = (dsumy - swa%wiy(ic)) * abs(swa%niy(ic))
2565 swa%bix(ic) = swa%wix(ic) * sign(
done, swa%nix(ic))
2566 swa%biy(ic) = swa%wiy(ic) * sign(
done, swa%niy(ic))
2567 dsumx = dsumx + swa%wix(ic) * abs(swa%nix(ic))
2568 dsumy = dsumy + swa%wiy(ic) * abs(swa%niy(ic))
2575 swa%bix(ic) = swa%bix(ic) * dsumx
2576 swa%biy(ic) = swa%biy(ic) * dsumy
2577 axy = axy + swa%bix(ic) * swa%niy(ic)
2578 ayx = ayx + swa%biy(ic) * swa%nix(ic)
2591 vx = vx + (swa%bix(ic) - axy * swa%biy(ic)) * swa%vi(ic)
2592 vy = vy + (swa%biy(ic) - ayx * swa%bix(ic)) * swa%vi(ic)
2594 denom =
done - axy * ayx
2595 if (denom /=
dzero)
then
2600 this%spdis(1, n) = vx
2601 this%spdis(2, n) = vy
2602 this%spdis(3, n) = vz
2613 integer(I4B),
intent(in) :: ibinun
2615 character(len=16) :: text
2616 character(len=16),
dimension(3) :: auxtxt
2618 integer(I4B) :: naux
2621 text =
' DATA-SPDIS'
2623 auxtxt(:) = [
' qx',
' qy',
' qz']
2624 call this%dis%record_srcdst_list_header(text, this%name_model, &
2625 this%packName, this%name_model, &
2626 this%packName, naux, auxtxt, ibinun, &
2627 this%dis%nodes, this%iout)
2630 do n = 1, this%dis%nodes
2631 call this%dis%record_mf6_list_entry(ibinun, n, n,
dzero, naux, &
2641 integer(I4B),
intent(in) :: ibinun
2643 character(len=16) :: text
2644 character(len=16),
dimension(1) :: auxtxt
2645 real(DP),
dimension(1) :: a
2647 integer(I4B) :: naux
2652 auxtxt(:) = [
' sat']
2653 call this%dis%record_srcdst_list_header(text, this%name_model, &
2654 this%packName, this%name_model, &
2655 this%packName, naux, auxtxt, ibinun, &
2656 this%dis%nodes, this%iout)
2659 do n = 1, this%dis%nodes
2661 call this%dis%record_mf6_list_entry(ibinun, n, n,
dzero, naux, a)
2673 integer(I4B),
intent(in) :: nedges
2675 this%nedges = this%nedges + nedges
2682 integer(I4B) :: max_conns
2684 integer(I4B) :: n, m, ic
2687 do n = 1, this%dis%nodes
2690 ic = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
2693 do m = 1, this%nedges
2694 if (this%nodedge(m) == n)
then
2700 if (ic > max_conns) max_conns = ic
2711 integer(I4B),
intent(in) :: nodedge
2712 integer(I4B),
intent(in) :: ihcedge
2713 real(DP),
intent(in) :: q
2714 real(DP),
intent(in) :: area
2715 real(DP),
intent(in) :: nx
2716 real(DP),
intent(in) :: ny
2717 real(DP),
intent(in) :: distance
2719 integer(I4B) :: lastedge
2721 this%lastedge = this%lastedge + 1
2722 lastedge = this%lastedge
2723 this%nodedge(lastedge) = nodedge
2724 this%ihcedge(lastedge) = ihcedge
2725 this%propsedge(1, lastedge) = q
2726 this%propsedge(2, lastedge) = area
2727 this%propsedge(3, lastedge) = nx
2728 this%propsedge(4, lastedge) = ny
2729 this%propsedge(5, lastedge) = distance
2733 if (this%lastedge == this%nedges) this%lastedge = 0
2739 integer(I4B) :: i, inode, iedge
2740 integer(I4B) :: n, start, end
2741 integer(I4B) :: prev_cnt, strt_idx, ipos
2743 do i = 1,
size(this%iedge_ptr)
2744 this%iedge_ptr(i) = 0
2746 do i = 1,
size(this%edge_idxs)
2747 this%edge_idxs(i) = 0
2751 do iedge = 1, this%nedges
2752 n = this%nodedge(iedge)
2753 this%iedge_ptr(n) = this%iedge_ptr(n) + 1
2757 prev_cnt = this%iedge_ptr(1)
2758 this%iedge_ptr(1) = 1
2759 do inode = 2, this%dis%nodes + 1
2760 strt_idx = this%iedge_ptr(inode - 1) + prev_cnt
2761 prev_cnt = this%iedge_ptr(inode)
2762 this%iedge_ptr(inode) = strt_idx
2766 do iedge = 1, this%nedges
2767 n = this%nodedge(iedge)
2768 start = this%iedge_ptr(n)
2769 end = this%iedge_ptr(n + 1) - 1
2770 do ipos = start,
end
2771 if (this%edge_idxs(ipos) > 0) cycle
2772 this%edge_idxs(ipos) = iedge
2788 real(dp) :: satthickness
2790 satthickness = thksatnm(this%ibound(n), &
2792 this%icelltype(n), &
2793 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.