39 integer(I4B),
pointer :: iname => null()
40 character(len=24),
dimension(:),
pointer :: aname => null()
41 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
42 real(dp),
dimension(:),
pointer,
contiguous :: hnew => null()
43 integer(I4B),
pointer :: ixt3d => null()
44 integer(I4B),
pointer :: ixt3drhs => null()
45 integer(I4B),
pointer :: iperched => null()
46 integer(I4B),
pointer :: ivarcv => null()
47 integer(I4B),
pointer :: idewatcv => null()
48 integer(I4B),
pointer :: ithickstrt => null()
49 integer(I4B),
pointer :: igwfnewtonur => null()
50 integer(I4B),
pointer :: icalcspdis => null()
51 integer(I4B),
pointer :: isavspdis => null()
52 integer(I4B),
pointer :: isavsat => null()
53 real(dp),
pointer :: hnoflo => null()
54 real(dp),
pointer :: satomega => null()
55 integer(I4B),
pointer :: irewet => null()
56 integer(I4B),
pointer :: iwetit => null()
57 integer(I4B),
pointer :: ihdwet => null()
58 integer(I4B),
pointer :: icellavg => null()
59 real(dp),
pointer :: wetfct => null()
60 real(dp),
pointer :: hdry => null()
61 integer(I4B),
dimension(:),
pointer,
contiguous :: icelltype => null()
62 integer(I4B),
dimension(:),
pointer,
contiguous :: ithickstartflag => null()
65 real(dp),
dimension(:),
pointer,
contiguous :: k11 => null()
66 real(dp),
dimension(:),
pointer,
contiguous :: k22 => null()
67 real(dp),
dimension(:),
pointer,
contiguous :: k33 => null()
68 real(dp),
dimension(:),
pointer,
contiguous :: k11input => null()
69 real(dp),
dimension(:),
pointer,
contiguous :: k22input => null()
70 real(dp),
dimension(:),
pointer,
contiguous :: k33input => null()
71 integer(I4B),
pointer :: iavgkeff => null()
72 integer(I4B),
pointer :: ik22 => null()
73 integer(I4B),
pointer :: ik33 => null()
74 integer(I4B),
pointer :: ik22overk => null()
75 integer(I4B),
pointer :: ik33overk => null()
76 integer(I4B),
pointer :: iangle1 => null()
77 integer(I4B),
pointer :: iangle2 => null()
78 integer(I4B),
pointer :: iangle3 => null()
79 real(dp),
dimension(:),
pointer,
contiguous :: angle1 => null()
80 real(dp),
dimension(:),
pointer,
contiguous :: angle2 => null()
81 real(dp),
dimension(:),
pointer,
contiguous :: angle3 => null()
83 integer(I4B),
pointer :: iwetdry => null()
84 real(dp),
dimension(:),
pointer,
contiguous :: wetdry => null()
85 real(dp),
dimension(:),
pointer,
contiguous :: sat => null()
86 real(dp),
dimension(:),
pointer,
contiguous :: condsat => null()
87 integer(I4B),
dimension(:),
pointer,
contiguous :: ibotnode => null()
89 real(dp),
dimension(:, :),
pointer,
contiguous :: spdis => null()
90 integer(I4B),
pointer :: nedges => null()
91 integer(I4B),
pointer :: lastedge => null()
92 integer(I4B),
dimension(:),
pointer,
contiguous :: nodedge => null()
93 integer(I4B),
dimension(:),
pointer,
contiguous :: ihcedge => null()
94 real(dp),
dimension(:, :),
pointer,
contiguous :: propsedge => null()
96 integer(I4B),
pointer :: intvk => null()
97 integer(I4B),
pointer :: invsc => null()
99 integer(I4B),
pointer :: kchangeper => null()
100 integer(I4B),
pointer :: kchangestp => null()
101 integer(I4B),
dimension(:),
pointer,
contiguous :: nodekchange => null()
152 subroutine npf_cr(npfobj, name_model, input_mempath, inunit, iout)
158 character(len=*),
intent(in) :: name_model
159 character(len=*),
intent(in) :: input_mempath
160 integer(I4B),
intent(in) :: inunit
161 integer(I4B),
intent(in) :: iout
163 character(len=*),
parameter :: fmtheader = &
164 "(1x, /1x, 'NPF -- NODE PROPERTY FLOW PACKAGE, VERSION 1, 3/30/2015', &
165 &' INPUT READ FROM MEMPATH: ', A, /)"
171 call npfobj%set_names(1, name_model,
'NPF',
'NPF', input_mempath)
174 call npfobj%allocate_scalars()
177 npfobj%inunit = inunit
184 write (iout, fmtheader) input_mempath
198 subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options)
206 integer(I4B),
intent(in) :: ingnc
207 integer(I4B),
intent(in) :: invsc
214 if (invsc > 0) this%invsc = invsc
216 if (.not.
present(npf_options))
then
219 call this%source_options()
222 call this%allocate_arrays(this%dis%nodes, this%dis%njas)
225 call this%source_griddata()
226 call this%prepcheck()
228 call this%set_options(npf_options)
231 call this%allocate_arrays(this%dis%nodes, this%dis%njas)
234 call this%check_options()
238 if (this%ixt3d /= 0) xt3d%ixt3d = this%ixt3d
239 call this%xt3d%xt3d_df(dis)
242 if (this%ixt3d /= 0 .and. ingnc > 0)
then
243 call store_error(
'Error in model '//trim(this%name_model)// &
244 '. The XT3D option cannot be used with the GNC &
245 &Package.', terminate=.true.)
259 integer(I4B),
intent(in) :: moffset
263 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)
288 subroutine npf_ar(this, ic, vsc, ibound, hnew)
293 type(
gwfictype),
pointer,
intent(in) :: ic
295 integer(I4B),
dimension(:),
pointer,
contiguous,
intent(inout) :: ibound
296 real(DP),
dimension(:),
pointer,
contiguous,
intent(inout) :: hnew
302 this%ibound => ibound
305 if (this%icalcspdis == 1)
then
306 call mem_reallocate(this%spdis, 3, this%dis%nodes,
'SPDIS', this%memoryPath)
307 call mem_reallocate(this%nodedge, this%nedges,
'NODEDGE', this%memoryPath)
308 call mem_reallocate(this%ihcedge, this%nedges,
'IHCEDGE', this%memoryPath)
309 call mem_reallocate(this%propsedge, 5, this%nedges,
'PROPSEDGE', &
311 do n = 1, this%dis%nodes
312 this%spdis(:, n) =
dzero
317 if (this%invsc /= 0)
then
322 if (this%invsc > 0)
then
334 call this%store_original_k_arrays(this%dis%nodes, this%dis%njas)
338 call this%preprocess_input()
341 if (this%ixt3d /= 0)
then
342 call this%xt3d%xt3d_ar(ibound, this%k11, this%ik33, this%k33, &
343 this%sat, this%ik22, this%k22, &
344 this%iangle1, this%iangle2, this%iangle3, &
345 this%angle1, this%angle2, this%angle3, &
346 this%inewton, this%icelltype)
350 if (this%intvk /= 0)
then
351 call this%tvk%ar(this%dis)
368 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.)
447 subroutine npf_cf(this, kiter, nodes, hnew)
450 integer(I4B) :: kiter
451 integer(I4B),
intent(in) :: nodes
452 real(DP),
intent(inout),
dimension(nodes) :: hnew
458 if (this%inewton /= 1)
then
459 call this%wd(kiter, hnew)
463 do n = 1, this%dis%nodes
464 if (this%icelltype(n) /= 0)
then
465 if (this%ibound(n) == 0)
then
468 call this%thksat(n, hnew(n), satn)
480 subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
485 integer(I4B) :: kiter
487 integer(I4B),
intent(in),
dimension(:) :: idxglo
488 real(DP),
intent(inout),
dimension(:) :: rhs
489 real(DP),
intent(inout),
dimension(:) :: hnew
491 integer(I4B) :: n, m, ii, idiag, ihc
492 integer(I4B) :: isymcon, idiagm
498 if (this%ixt3d /= 0)
then
499 call this%xt3d%xt3d_fc(kiter, matrix_sln, idxglo, rhs, hnew)
501 do n = 1, this%dis%nodes
502 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
503 if (this%dis%con%mask(ii) == 0) cycle
505 m = this%dis%con%ja(ii)
510 ihc = this%dis%con%ihc(this%dis%con%jas(ii))
511 hyn = this%hy_eff(n, m, ihc, ipos=ii)
512 hym = this%hy_eff(m, n, ihc, ipos=ii)
515 if (ihc == c3d_vertical)
then
518 cond =
vcond(this%ibound(n), this%ibound(m), &
519 this%icelltype(n), this%icelltype(m), this%inewton, &
520 this%ivarcv, this%idewatcv, &
521 this%condsat(this%dis%con%jas(ii)), hnew(n), hnew(m), &
523 this%sat(n), this%sat(m), &
524 this%dis%top(n), this%dis%top(m), &
525 this%dis%bot(n), this%dis%bot(m), &
526 this%dis%con%hwva(this%dis%con%jas(ii)))
529 if (this%iperched /= 0)
then
530 if (this%icelltype(m) /= 0)
then
531 if (hnew(m) < this%dis%top(m))
then
534 idiag = this%dis%con%ia(n)
535 rhs(n) = rhs(n) - cond * this%dis%bot(n)
536 call matrix_sln%add_value_pos(idxglo(idiag), -cond)
539 isymcon = this%dis%con%isym(ii)
540 call matrix_sln%add_value_pos(idxglo(isymcon), cond)
541 rhs(m) = rhs(m) + cond * this%dis%bot(n)
552 cond =
hcond(this%ibound(n), this%ibound(m), &
553 this%icelltype(n), this%icelltype(m), &
555 this%dis%con%ihc(this%dis%con%jas(ii)), &
557 this%condsat(this%dis%con%jas(ii)), &
558 hnew(n), hnew(m), this%sat(n), this%sat(m), hyn, hym, &
559 this%dis%top(n), this%dis%top(m), &
560 this%dis%bot(n), this%dis%bot(m), &
561 this%dis%con%cl1(this%dis%con%jas(ii)), &
562 this%dis%con%cl2(this%dis%con%jas(ii)), &
563 this%dis%con%hwva(this%dis%con%jas(ii)))
567 idiag = this%dis%con%ia(n)
568 call matrix_sln%add_value_pos(idxglo(ii), cond)
569 call matrix_sln%add_value_pos(idxglo(idiag), -cond)
572 isymcon = this%dis%con%isym(ii)
573 idiagm = this%dis%con%ia(m)
574 call matrix_sln%add_value_pos(idxglo(isymcon), cond)
575 call matrix_sln%add_value_pos(idxglo(idiagm), -cond)
587 subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew)
590 integer(I4B) :: kiter
592 integer(I4B),
intent(in),
dimension(:) :: idxglo
593 real(DP),
intent(inout),
dimension(:) :: rhs
594 real(DP),
intent(inout),
dimension(:) :: hnew
596 integer(I4B) :: nodes, nja
597 integer(I4B) :: n, m, ii, idiag
598 integer(I4B) :: isymcon, idiagm
603 real(DP) :: filledterm
611 nodes = this%dis%nodes
612 nja = this%dis%con%nja
613 if (this%ixt3d /= 0)
then
614 call this%xt3d%xt3d_fn(kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew)
618 idiag = this%dis%con%ia(n)
619 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
620 if (this%dis%con%mask(ii) == 0) cycle
622 m = this%dis%con%ja(ii)
623 isymcon = this%dis%con%isym(ii)
626 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0 .and. &
627 this%ivarcv == 0)
then
633 if (hnew(m) < hnew(n)) iups = n
635 if (iups == n) idn = m
638 if (this%icelltype(iups) == 0) cycle
642 topup = this%dis%top(iups)
643 botup = this%dis%bot(iups)
644 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 2)
then
645 topup = min(this%dis%top(n), this%dis%top(m))
646 botup = max(this%dis%bot(n), this%dis%bot(m))
650 cond = this%condsat(this%dis%con%jas(ii))
653 consterm = -cond * (hnew(iups) - hnew(idn))
655 filledterm = matrix_sln%get_value_pos(idxglo(ii))
658 idiagm = this%dis%con%ia(m)
663 term = consterm * derv
664 rhs(n) = rhs(n) + term * hnew(n)
665 rhs(m) = rhs(m) - term * hnew(n)
667 call matrix_sln%add_value_pos(idxglo(idiag), term)
669 if (this%ibound(n) > 0)
then
670 filledterm = matrix_sln%get_value_pos(idxglo(ii))
671 call matrix_sln%set_value_pos(idxglo(ii), filledterm)
674 filledterm = matrix_sln%get_value_pos(idxglo(idiagm))
675 call matrix_sln%set_value_pos(idxglo(idiagm), filledterm)
677 if (this%ibound(m) > 0)
then
678 call matrix_sln%add_value_pos(idxglo(isymcon), -term)
683 term = -consterm * derv
684 rhs(n) = rhs(n) + term * hnew(m)
685 rhs(m) = rhs(m) - term * hnew(m)
687 filledterm = matrix_sln%get_value_pos(idxglo(idiag))
688 call matrix_sln%set_value_pos(idxglo(idiag), filledterm)
690 if (this%ibound(n) > 0)
then
691 call matrix_sln%add_value_pos(idxglo(ii), term)
694 call matrix_sln%add_value_pos(idxglo(idiagm), -term)
696 if (this%ibound(m) > 0)
then
697 filledterm = matrix_sln%get_value_pos(idxglo(isymcon))
698 call matrix_sln%set_value_pos(idxglo(isymcon), filledterm)
717 subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
720 integer(I4B),
intent(in) :: neqmod
721 real(DP),
dimension(neqmod),
intent(inout) :: x
722 real(DP),
dimension(neqmod),
intent(in) :: xtemp
723 real(DP),
dimension(neqmod),
intent(inout) :: dx
724 integer(I4B),
intent(inout) :: inewtonur
725 real(DP),
intent(inout) :: dxmax
726 integer(I4B),
intent(inout) :: locmax
734 do n = 1, this%dis%nodes
735 if (this%ibound(n) < 1) cycle
736 if (this%icelltype(n) > 0)
then
737 botm = this%dis%bot(this%ibotnode(n))
740 if (x(n) < botm)
then
744 if (abs(dxx) > abs(dxmax))
then
763 real(DP),
intent(inout),
dimension(:) :: hnew
764 real(DP),
intent(inout),
dimension(:) :: flowja
766 integer(I4B) :: n, ipos, m
771 if (this%ixt3d /= 0)
then
772 call this%xt3d%xt3d_flowja(hnew, flowja)
775 do n = 1, this%dis%nodes
776 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
777 m = this%dis%con%ja(ipos)
779 call this%qcalc(n, m, hnew(n), hnew(m), ipos, qnm)
781 flowja(this%dis%con%isym(ipos)) = -qnm
796 integer(I4B),
intent(in) :: n
797 real(DP),
intent(in) :: hn
798 real(DP),
intent(inout) :: thksat
801 if (hn >= this%dis%top(n))
then
804 thksat = (hn - this%dis%bot(n)) / (this%dis%top(n) - this%dis%bot(n))
808 if (this%inewton /= 0)
then
822 integer(I4B),
intent(in) :: n
823 integer(I4B),
intent(in) :: m
824 real(DP),
intent(in) :: hn
825 real(DP),
intent(in) :: hm
826 integer(I4B),
intent(in) :: icon
827 real(DP),
intent(inout) :: qnm
831 real(DP) :: hntemp, hmtemp
835 ihc = this%dis%con%ihc(this%dis%con%jas(icon))
836 hyn = this%hy_eff(n, m, ihc, ipos=icon)
837 hym = this%hy_eff(m, n, ihc, ipos=icon)
841 condnm =
vcond(this%ibound(n), this%ibound(m), &
842 this%icelltype(n), this%icelltype(m), this%inewton, &
843 this%ivarcv, this%idewatcv, &
844 this%condsat(this%dis%con%jas(icon)), hn, hm, &
846 this%sat(n), this%sat(m), &
847 this%dis%top(n), this%dis%top(m), &
848 this%dis%bot(n), this%dis%bot(m), &
849 this%dis%con%hwva(this%dis%con%jas(icon)))
851 condnm =
hcond(this%ibound(n), this%ibound(m), &
852 this%icelltype(n), this%icelltype(m), &
854 this%dis%con%ihc(this%dis%con%jas(icon)), &
856 this%condsat(this%dis%con%jas(icon)), &
857 hn, hm, this%sat(n), this%sat(m), hyn, hym, &
858 this%dis%top(n), this%dis%top(m), &
859 this%dis%bot(n), this%dis%bot(m), &
860 this%dis%con%cl1(this%dis%con%jas(icon)), &
861 this%dis%con%cl2(this%dis%con%jas(icon)), &
862 this%dis%con%hwva(this%dis%con%jas(icon)))
870 if (this%iperched /= 0)
then
871 if (this%dis%con%ihc(this%dis%con%jas(icon)) == 0)
then
873 if (this%icelltype(n) /= 0)
then
874 if (hn < this%dis%top(n)) hntemp = this%dis%bot(m)
877 if (this%icelltype(m) /= 0)
then
878 if (hm < this%dis%top(m)) hmtemp = this%dis%bot(n)
885 qnm = condnm * (hmtemp - hntemp)
896 real(DP),
dimension(:),
intent(in) :: flowja
897 integer(I4B),
intent(in) :: icbcfl
898 integer(I4B),
intent(in) :: icbcun
900 integer(I4B) :: ibinun
903 if (this%ipakcb < 0)
then
905 elseif (this%ipakcb == 0)
then
910 if (icbcfl == 0) ibinun = 0
913 if (ibinun /= 0)
then
914 call this%dis%record_connection_array(flowja, ibinun, this%iout)
918 if (this%isavspdis /= 0)
then
919 if (ibinun /= 0)
call this%sav_spdis(ibinun)
923 if (this%isavsat /= 0)
then
924 if (ibinun /= 0)
call this%sav_sat(ibinun)
939 integer(I4B),
intent(in) :: ibudfl
940 real(DP),
intent(inout),
dimension(:) :: flowja
942 character(len=LENBIGLINE) :: line
943 character(len=30) :: tempstr
944 integer(I4B) :: n, ipos, m
947 character(len=*),
parameter :: fmtiprflow = &
948 &
"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)"
951 if (ibudfl /= 0 .and. this%iprflow > 0)
then
952 write (this%iout, fmtiprflow)
kper,
kstp
953 do n = 1, this%dis%nodes
955 call this%dis%noder_to_string(n, tempstr)
956 line = trim(tempstr)//
':'
957 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
958 m = this%dis%con%ja(ipos)
959 call this%dis%noder_to_string(m, tempstr)
960 line = trim(line)//
' '//trim(tempstr)
962 write (tempstr,
'(1pg15.6)') qnm
963 line = trim(line)//
' '//trim(adjustl(tempstr))
965 write (this%iout,
'(a)') trim(line)
986 if (this%intvk /= 0)
then
988 deallocate (this%tvk)
992 if (this%invsc /= 0)
then
1035 deallocate (this%aname)
1057 call this%NumericalPackageType%da()
1076 call this%NumericalPackageType%allocate_scalars()
1079 call mem_allocate(this%iname,
'INAME', this%memoryPath)
1080 call mem_allocate(this%ixt3d,
'IXT3D', this%memoryPath)
1081 call mem_allocate(this%ixt3drhs,
'IXT3DRHS', this%memoryPath)
1082 call mem_allocate(this%satomega,
'SATOMEGA', this%memoryPath)
1083 call mem_allocate(this%hnoflo,
'HNOFLO', this%memoryPath)
1085 call mem_allocate(this%icellavg,
'ICELLAVG', this%memoryPath)
1086 call mem_allocate(this%iavgkeff,
'IAVGKEFF', this%memoryPath)
1089 call mem_allocate(this%ik22overk,
'IK22OVERK', this%memoryPath)
1090 call mem_allocate(this%ik33overk,
'IK33OVERK', this%memoryPath)
1091 call mem_allocate(this%iperched,
'IPERCHED', this%memoryPath)
1092 call mem_allocate(this%ivarcv,
'IVARCV', this%memoryPath)
1093 call mem_allocate(this%idewatcv,
'IDEWATCV', this%memoryPath)
1094 call mem_allocate(this%ithickstrt,
'ITHICKSTRT', this%memoryPath)
1095 call mem_allocate(this%icalcspdis,
'ICALCSPDIS', this%memoryPath)
1096 call mem_allocate(this%isavspdis,
'ISAVSPDIS', this%memoryPath)
1097 call mem_allocate(this%isavsat,
'ISAVSAT', this%memoryPath)
1098 call mem_allocate(this%irewet,
'IREWET', this%memoryPath)
1099 call mem_allocate(this%wetfct,
'WETFCT', this%memoryPath)
1100 call mem_allocate(this%iwetit,
'IWETIT', this%memoryPath)
1101 call mem_allocate(this%ihdwet,
'IHDWET', this%memoryPath)
1102 call mem_allocate(this%iangle1,
'IANGLE1', this%memoryPath)
1103 call mem_allocate(this%iangle2,
'IANGLE2', this%memoryPath)
1104 call mem_allocate(this%iangle3,
'IANGLE3', this%memoryPath)
1105 call mem_allocate(this%iwetdry,
'IWETDRY', this%memoryPath)
1106 call mem_allocate(this%nedges,
'NEDGES', this%memoryPath)
1107 call mem_allocate(this%lastedge,
'LASTEDGE', this%memoryPath)
1108 call mem_allocate(this%intvk,
'INTVK', this%memoryPath)
1109 call mem_allocate(this%invsc,
'INVSC', this%memoryPath)
1110 call mem_allocate(this%kchangeper,
'KCHANGEPER', this%memoryPath)
1111 call mem_allocate(this%kchangestp,
'KCHANGESTP', this%memoryPath)
1114 call mem_setptr(this%igwfnewtonur,
'INEWTONUR', &
1121 this%satomega =
dzero
1153 this%iasym = this%inewton
1174 integer(I4B),
intent(in) :: ncells
1175 integer(I4B),
intent(in) :: njas
1181 this%k11input(n) = this%k11(n)
1182 this%k22input(n) = this%k22(n)
1183 this%k33input(n) = this%k33(n)
1195 integer(I4B),
intent(in) :: ncells
1196 integer(I4B),
intent(in) :: njas
1200 call mem_allocate(this%ithickstartflag, ncells,
'ITHICKSTARTFLAG', &
1202 call mem_allocate(this%icelltype, ncells,
'ICELLTYPE', this%memoryPath)
1203 call mem_allocate(this%k11, ncells,
'K11', this%memoryPath)
1204 call mem_allocate(this%sat, ncells,
'SAT', this%memoryPath)
1205 call mem_allocate(this%condsat, njas,
'CONDSAT', this%memoryPath)
1208 call mem_allocate(this%k22, ncells,
'K22', this%memoryPath)
1209 call mem_allocate(this%k33, ncells,
'K33', this%memoryPath)
1210 call mem_allocate(this%wetdry, ncells,
'WETDRY', this%memoryPath)
1211 call mem_allocate(this%angle1, ncells,
'ANGLE1', this%memoryPath)
1212 call mem_allocate(this%angle2, ncells,
'ANGLE2', this%memoryPath)
1213 call mem_allocate(this%angle3, ncells,
'ANGLE3', this%memoryPath)
1216 call mem_allocate(this%ibotnode, 0,
'IBOTNODE', this%memoryPath)
1217 call mem_allocate(this%nodedge, 0,
'NODEDGE', this%memoryPath)
1218 call mem_allocate(this%ihcedge, 0,
'IHCEDGE', this%memoryPath)
1219 call mem_allocate(this%propsedge, 0, 0,
'PROPSEDGE', this%memoryPath)
1222 call mem_allocate(this%k11input, 0,
'K11INPUT', this%memoryPath)
1223 call mem_allocate(this%k22input, 0,
'K22INPUT', this%memoryPath)
1224 call mem_allocate(this%k33input, 0,
'K33INPUT', this%memoryPath)
1227 call mem_allocate(this%spdis, 3, 0,
'SPDIS', this%memoryPath)
1230 call mem_allocate(this%nodekchange, ncells,
'NODEKCHANGE', this%memoryPath)
1234 this%angle1(n) =
dzero
1235 this%angle2(n) =
dzero
1236 this%angle3(n) =
dzero
1237 this%wetdry(n) =
dzero
1238 this%nodekchange(n) =
dzero
1242 allocate (this%aname(this%iname))
1243 this%aname = [
' ICELLTYPE',
' K', &
1245 ' WETDRY',
' ANGLE1', &
1246 ' ANGLE2',
' ANGLE3']
1263 write (this%iout,
'(1x,a)')
'Setting NPF Options'
1264 if (found%iprflow) &
1265 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be printed &
1266 &to listing file whenever ICBCFL is not zero.'
1268 write (this%iout,
'(4x,a)')
'Cell-by-cell flow information will be saved &
1269 &to binary file whenever ICBCFL is not zero.'
1270 if (found%cellavg) &
1271 write (this%iout,
'(4x,a,i0)')
'Alternative cell averaging [1=logarithmic, &
1272 &2=AMT-LMK, 3=AMT-HMK] set to: ', &
1274 if (found%ithickstrt) &
1275 write (this%iout,
'(4x,a)')
'THICKSTRT option has been activated.'
1276 if (found%iperched) &
1277 write (this%iout,
'(4x,a)')
'Vertical flow will be adjusted for perched &
1280 write (this%iout,
'(4x,a)')
'Vertical conductance varies with water table.'
1281 if (found%idewatcv) &
1282 write (this%iout,
'(4x,a)')
'Vertical conductance accounts for dewatered &
1283 &portion of an underlying cell.'
1284 if (found%ixt3d)
write (this%iout,
'(4x,a)')
'XT3D formulation is selected.'
1285 if (found%ixt3drhs) &
1286 write (this%iout,
'(4x,a)')
'XT3D RHS formulation is selected.'
1287 if (found%isavspdis) &
1288 write (this%iout,
'(4x,a)')
'Specific discharge will be calculated at cell &
1289 ¢ers and written to DATA-SPDIS in budget &
1290 &file when requested.'
1291 if (found%isavsat) &
1292 write (this%iout,
'(4x,a)')
'Saturation will be written to DATA-SAT in &
1293 &budget file when requested.'
1294 if (found%ik22overk) &
1295 write (this%iout,
'(4x,a)')
'Values specified for K22 are anisotropy &
1296 &ratios and will be multiplied by K before &
1297 &being used in calculations.'
1298 if (found%ik33overk) &
1299 write (this%iout,
'(4x,a)')
'Values specified for K33 are anisotropy &
1300 &ratios and will be multiplied by K before &
1301 &being used in calculations.'
1302 if (found%inewton) &
1303 write (this%iout,
'(4x,a)')
'NEWTON-RAPHSON method disabled for unconfined &
1305 if (found%satomega) &
1306 write (this%iout,
'(4x,a,1pg15.6)')
'Saturation omega: ', this%satomega
1308 write (this%iout,
'(4x,a)')
'Rewetting is active.'
1310 write (this%iout,
'(4x,a,1pg15.6)') &
1311 'Wetting factor (WETFCT) has been set to: ', this%wetfct
1313 write (this%iout,
'(4x,a,i5)') &
1314 'Wetting iteration interval (IWETIT) has been set to: ', this%iwetit
1316 write (this%iout,
'(4x,a,i5)') &
1317 'Head rewet equation (IHDWET) has been set to: ', this%ihdwet
1318 write (this%iout,
'(1x,a,/)')
'End Setting NPF Options'
1337 character(len=LENVARNAME),
dimension(3) :: cellavg_method = &
1338 &[character(len=LENVARNAME) ::
'LOGARITHMIC',
'AMT-LMK',
'AMT-HMK']
1340 character(len=LINELENGTH) :: tvk6_filename
1343 call mem_set_value(this%iprflow,
'IPRFLOW', this%input_mempath, found%iprflow)
1344 call mem_set_value(this%ipakcb,
'IPAKCB', this%input_mempath, found%ipakcb)
1345 call mem_set_value(this%icellavg,
'CELLAVG', this%input_mempath, &
1346 cellavg_method, found%cellavg)
1347 call mem_set_value(this%ithickstrt,
'ITHICKSTRT', this%input_mempath, &
1349 call mem_set_value(this%iperched,
'IPERCHED', this%input_mempath, &
1351 call mem_set_value(this%ivarcv,
'IVARCV', this%input_mempath, found%ivarcv)
1352 call mem_set_value(this%idewatcv,
'IDEWATCV', this%input_mempath, &
1354 call mem_set_value(this%ixt3d,
'IXT3D', this%input_mempath, found%ixt3d)
1355 call mem_set_value(this%ixt3drhs,
'IXT3DRHS', this%input_mempath, &
1357 call mem_set_value(this%isavspdis,
'ISAVSPDIS', this%input_mempath, &
1359 call mem_set_value(this%isavsat,
'ISAVSAT', this%input_mempath, found%isavsat)
1360 call mem_set_value(this%ik22overk,
'IK22OVERK', this%input_mempath, &
1362 call mem_set_value(this%ik33overk,
'IK33OVERK', this%input_mempath, &
1364 call mem_set_value(this%inewton,
'INEWTON', this%input_mempath, found%inewton)
1365 call mem_set_value(this%satomega,
'SATOMEGA', this%input_mempath, &
1367 call mem_set_value(this%irewet,
'IREWET', this%input_mempath, found%irewet)
1368 call mem_set_value(this%wetfct,
'WETFCT', this%input_mempath, found%wetfct)
1369 call mem_set_value(this%iwetit,
'IWETIT', this%input_mempath, found%iwetit)
1370 call mem_set_value(this%ihdwet,
'IHDWET', this%input_mempath, found%ihdwet)
1373 if (found%ipakcb) this%ipakcb = -1
1376 if (found%ixt3d .and. found%ixt3drhs) this%ixt3d = 2
1379 if (found%isavspdis) this%icalcspdis = this%isavspdis
1382 if (found%inewton)
then
1388 if (
filein_fname(tvk6_filename,
'TVK6_FILENAME', this%input_mempath, &
1389 this%input_fname))
then
1390 call openfile(this%intvk, this%iout, tvk6_filename,
'TVK')
1391 call tvk_cr(this%tvk, this%name_model, this%intvk, this%iout)
1395 if (this%iout > 0)
then
1396 call this%log_options(found)
1410 this%icellavg = options%icellavg
1411 this%ithickstrt = options%ithickstrt
1412 this%iperched = options%iperched
1413 this%ivarcv = options%ivarcv
1414 this%idewatcv = options%idewatcv
1415 this%irewet = options%irewet
1416 this%wetfct = options%wetfct
1417 this%iwetit = options%iwetit
1418 this%ihdwet = options%ihdwet
1434 if (this%inewton > 0)
then
1435 this%satomega = dem6
1438 if (this%inewton > 0)
then
1439 if (this%iperched > 0)
then
1440 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1441 'BE USED WITH PERCHED OPTION.'
1444 if (this%ivarcv > 0)
then
1445 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1446 'BE USED WITH VARIABLECV OPTION.'
1449 if (this%irewet > 0)
then
1450 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1451 'BE USED WITH REWET OPTION.'
1456 if (this%ixt3d /= 0)
then
1457 if (this%icellavg > 0)
then
1458 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. '// &
1459 'ALTERNATIVE_CELL_AVERAGING OPTION '// &
1460 'CANNOT BE USED WITH XT3D OPTION.'
1463 if (this%ithickstrt > 0)
then
1464 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. THICKSTRT OPTION '// &
1465 'CANNOT BE USED WITH XT3D OPTION.'
1468 if (this%iperched > 0)
then
1469 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. PERCHED OPTION '// &
1470 'CANNOT BE USED WITH XT3D OPTION.'
1473 if (this%ivarcv > 0)
then
1474 write (
errmsg,
'(a)')
'ERROR IN NPF OPTIONS. VARIABLECV OPTION '// &
1475 'CANNOT BE USED WITH XT3D OPTION.'
1498 write (this%iout,
'(1x,a)')
'Setting NPF Griddata'
1500 if (found%icelltype)
then
1501 write (this%iout,
'(4x,a)')
'ICELLTYPE set from input file'
1505 write (this%iout,
'(4x,a)')
'K set from input file'
1509 write (this%iout,
'(4x,a)')
'K33 set from input file'
1511 write (this%iout,
'(4x,a)')
'K33 not provided. Setting K33 = K.'
1515 write (this%iout,
'(4x,a)')
'K22 set from input file'
1517 write (this%iout,
'(4x,a)')
'K22 not provided. Setting K22 = K.'
1520 if (found%wetdry)
then
1521 write (this%iout,
'(4x,a)')
'WETDRY set from input file'
1524 if (found%angle1)
then
1525 write (this%iout,
'(4x,a)')
'ANGLE1 set from input file'
1528 if (found%angle2)
then
1529 write (this%iout,
'(4x,a)')
'ANGLE2 set from input file'
1532 if (found%angle3)
then
1533 write (this%iout,
'(4x,a)')
'ANGLE3 set from input file'
1536 write (this%iout,
'(1x,a,/)')
'End Setting NPF Griddata'
1553 character(len=LINELENGTH) :: errmsg
1555 logical,
dimension(2) :: afound
1556 integer(I4B),
dimension(:),
pointer,
contiguous :: map
1560 if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
1563 call mem_set_value(this%icelltype,
'ICELLTYPE', this%input_mempath, map, &
1565 call mem_set_value(this%k11,
'K', this%input_mempath, map, found%k)
1566 call mem_set_value(this%k33,
'K33', this%input_mempath, map, found%k33)
1567 call mem_set_value(this%k22,
'K22', this%input_mempath, map, found%k22)
1568 call mem_set_value(this%wetdry,
'WETDRY', this%input_mempath, map, &
1570 call mem_set_value(this%angle1,
'ANGLE1', this%input_mempath, map, &
1572 call mem_set_value(this%angle2,
'ANGLE2', this%input_mempath, map, &
1574 call mem_set_value(this%angle3,
'ANGLE3', this%input_mempath, map, &
1578 if (.not. found%icelltype)
then
1579 write (errmsg,
'(a)')
'Error in GRIDDATA block: ICELLTYPE not found.'
1584 if (.not. found%k)
then
1585 write (errmsg,
'(a)')
'Error in GRIDDATA block: K not found.'
1590 if (.not. found%k33 .and. this%ik33overk /= 0)
then
1591 write (errmsg,
'(a)')
'K33OVERK option specified but K33 not specified.'
1596 if (.not. found%k22 .and. this%ik22overk /= 0)
then
1597 write (errmsg,
'(a)')
'K22OVERK option specified but K22 not specified.'
1602 if (found%k33) this%ik33 = 1
1603 if (found%k22) this%ik22 = 1
1604 if (found%wetdry) this%iwetdry = 1
1605 if (found%angle1) this%iangle1 = 1
1606 if (found%angle2) this%iangle2 = 1
1607 if (found%angle3) this%iangle3 = 1
1610 if (.not. found%k33)
then
1611 call mem_set_value(this%k33,
'K', this%input_mempath, map, afound(1))
1613 if (.not. found%k22)
then
1614 call mem_set_value(this%k22,
'K', this%input_mempath, map, afound(2))
1616 if (.not. found%wetdry)
call mem_reallocate(this%wetdry, 1,
'WETDRY', &
1617 trim(this%memoryPath))
1618 if (.not. found%angle1 .and. this%ixt3d == 0) &
1619 call mem_reallocate(this%angle1, 0,
'ANGLE1', trim(this%memoryPath))
1620 if (.not. found%angle2 .and. this%ixt3d == 0) &
1621 call mem_reallocate(this%angle2, 0,
'ANGLE2', trim(this%memoryPath))
1622 if (.not. found%angle3 .and. this%ixt3d == 0) &
1623 call mem_reallocate(this%angle3, 0,
'ANGLE3', trim(this%memoryPath))
1626 if (this%iout > 0)
then
1627 call this%log_griddata(found)
1643 character(len=24),
dimension(:),
pointer :: aname
1644 character(len=LINELENGTH) :: cellstr, errmsg
1645 integer(I4B) :: nerr, n
1647 character(len=*),
parameter :: fmtkerr = &
1648 &
"(1x, 'Hydraulic property ',a,' is <= 0 for cell ',a, ' ', 1pg15.6)"
1649 character(len=*),
parameter :: fmtkerr2 = &
1650 &
"(1x, '... ', i0,' additional errors not shown for ',a)"
1657 do n = 1,
size(this%k11)
1658 if (this%k11(n) <= dzero)
then
1660 if (nerr <= 20)
then
1661 call this%dis%noder_to_string(n, cellstr)
1662 write (errmsg, fmtkerr) trim(adjustl(aname(2))), trim(cellstr), &
1669 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(2)))
1674 if (this%ik33 /= 0)
then
1678 do n = 1,
size(this%k33)
1679 if (this%ik33overk /= 0) this%k33(n) = this%k33(n) * this%k11(n)
1680 if (this%k33(n) <= dzero)
then
1682 if (nerr <= 20)
then
1683 call this%dis%noder_to_string(n, cellstr)
1684 write (errmsg, fmtkerr) trim(adjustl(aname(3))), trim(cellstr), &
1691 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(3)))
1697 if (this%ik22 /= 0)
then
1700 if (this%dis%con%ianglex == 0)
then
1701 write (errmsg,
'(a)')
'Error. ANGLDEGX not provided in '// &
1702 'discretization file, but K22 was specified. '
1708 do n = 1,
size(this%k22)
1709 if (this%ik22overk /= 0) this%k22(n) = this%k22(n) * this%k11(n)
1710 if (this%k22(n) <= dzero)
then
1712 if (nerr <= 20)
then
1713 call this%dis%noder_to_string(n, cellstr)
1714 write (errmsg, fmtkerr) trim(adjustl(aname(4))), trim(cellstr), &
1721 write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(4)))
1727 if (this%irewet == 1)
then
1728 if (this%iwetdry == 0)
then
1729 write (errmsg,
'(a, a, a)')
'Error in GRIDDATA block: ', &
1730 trim(adjustl(aname(5))),
' not found.'
1736 if (this%iangle1 /= 0)
then
1737 do n = 1,
size(this%angle1)
1738 this%angle1(n) = this%angle1(n) *
dpio180
1741 if (this%ixt3d /= 0)
then
1743 write (this%iout,
'(a)')
'XT3D IN USE, BUT ANGLE1 NOT SPECIFIED. '// &
1744 'SETTING ANGLE1 TO ZERO.'
1745 do n = 1,
size(this%angle1)
1746 this%angle1(n) = dzero
1750 if (this%iangle2 /= 0)
then
1751 if (this%iangle1 == 0)
then
1752 write (errmsg,
'(a)')
'ANGLE2 SPECIFIED BUT NOT ANGLE1. '// &
1753 'ANGLE2 REQUIRES ANGLE1. '
1756 if (this%iangle3 == 0)
then
1757 write (errmsg,
'(a)')
'ANGLE2 SPECIFIED BUT NOT ANGLE3. '// &
1758 'SPECIFY BOTH OR NEITHER ONE. '
1761 do n = 1,
size(this%angle2)
1762 this%angle2(n) = this%angle2(n) *
dpio180
1765 if (this%iangle3 /= 0)
then
1766 if (this%iangle1 == 0)
then
1767 write (errmsg,
'(a)')
'ANGLE3 SPECIFIED BUT NOT ANGLE1. '// &
1768 'ANGLE3 REQUIRES ANGLE1. '
1771 if (this%iangle2 == 0)
then
1772 write (errmsg,
'(a)')
'ANGLE3 SPECIFIED BUT NOT ANGLE2. '// &
1773 'SPECIFY BOTH OR NEITHER ONE. '
1776 do n = 1,
size(this%angle3)
1777 this%angle3(n) = this%angle3(n) *
dpio180
1807 integer(I4B) :: n, m, ii, nn
1808 real(DP) :: hyn, hym
1809 real(DP) :: satn, topn, botn
1810 integer(I4B) :: nextn
1811 real(DP) :: minbot, botm
1813 character(len=LINELENGTH) :: cellstr, errmsg
1815 character(len=*),
parameter :: fmtcnv = &
1817 &' ELIMINATED BECAUSE ALL HYDRAULIC CONDUCTIVITIES TO NODE ARE 0.')"
1818 character(len=*),
parameter :: fmtnct = &
1819 &
"(1X,'Negative cell thickness at cell ', A)"
1820 character(len=*),
parameter :: fmtihbe = &
1821 &
"(1X,'Initial head, bottom elevation:',1P,2G13.5)"
1822 character(len=*),
parameter :: fmttebe = &
1823 &
"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
1825 do n = 1, this%dis%nodes
1826 this%ithickstartflag(n) = 0
1832 nodeloop:
do n = 1, this%dis%nodes
1835 if (this%ibound(n) == 0)
then
1836 if (this%irewet /= 0)
then
1837 if (this%wetdry(n) == dzero) cycle nodeloop
1844 if (this%k11(n) /= dzero) cycle nodeloop
1848 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1849 m = this%dis%con%ja(ii)
1850 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0)
then
1852 if (this%ik33 /= 0) hyn = this%k33(n)
1853 if (hyn /= dzero)
then
1855 if (this%ik33 /= 0) hym = this%k33(m)
1856 if (hym /= dzero) cycle
1864 this%hnew(n) = this%hnoflo
1865 if (this%irewet /= 0) this%wetdry(n) = dzero
1866 call this%dis%noder_to_string(n, cellstr)
1867 write (this%iout, fmtcnv) trim(adjustl(cellstr))
1872 if (this%inewton == 0)
then
1875 call this%wd(0, this%hnew)
1881 if (this%ivarcv == 1)
then
1882 do n = 1, this%dis%nodes
1883 if (this%hnew(n) < this%dis%bot(n))
then
1884 this%hnew(n) = this%dis%bot(n) + dem6
1892 if (this%ithickstrt == 0)
then
1893 do n = 1, this%dis%nodes
1894 if (this%icelltype(n) < 0)
then
1895 this%icelltype(n) = 1
1905 do n = 1, this%dis%nodes
1906 if (this%ibound(n) == 0)
then
1908 if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0)
then
1909 this%ithickstartflag(n) = 1
1910 this%icelltype(n) = 0
1913 topn = this%dis%top(n)
1914 botn = this%dis%bot(n)
1915 if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0)
then
1916 call this%thksat(n, this%ic%strt(n), satn)
1917 if (botn > this%ic%strt(n))
then
1918 call this%dis%noder_to_string(n, cellstr)
1919 write (errmsg, fmtnct) trim(adjustl(cellstr))
1921 write (errmsg, fmtihbe) this%ic%strt(n), botn
1924 this%ithickstartflag(n) = 1
1925 this%icelltype(n) = 0
1928 if (botn > topn)
then
1929 call this%dis%noder_to_string(n, cellstr)
1930 write (errmsg, fmtnct) trim(adjustl(cellstr))
1932 write (errmsg, fmttebe) topn, botn
1945 if (this%ixt3d == 0)
then
1950 do n = 1, this%dis%nodes
1951 call this%calc_condsat(n, .true.)
1957 if (this%igwfnewtonur /= 0)
then
1959 trim(this%memoryPath))
1960 do n = 1, this%dis%nodes
1962 minbot = this%dis%bot(n)
1965 do while (.not. finished)
1969 do ii = this%dis%con%ia(nn) + 1, this%dis%con%ia(nn + 1) - 1
1972 m = this%dis%con%ja(ii)
1973 botm = this%dis%bot(m)
1976 if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0)
then
1977 if (m > nn .and. botm < minbot)
then
1989 this%ibotnode(n) = nn
1994 this%igwfnewtonur => null()
2008 integer(I4B),
intent(in) :: node
2009 logical,
intent(in) :: upperOnly
2011 integer(I4B) :: ii, m, n, ihc, jj
2012 real(DP) :: topm, topn, topnode, botm, botn, botnode, satm, satn, satnode
2013 real(DP) :: hyn, hym, hn, hm, fawidth, csat
2015 satnode = this%calc_initial_sat(node)
2017 topnode = this%dis%top(node)
2018 botnode = this%dis%bot(node)
2021 do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
2025 m = this%dis%con%ja(ii)
2026 jj = this%dis%con%jas(ii)
2028 if (upperonly) cycle
2035 topn = this%dis%top(n)
2036 botn = this%dis%bot(n)
2037 satn = this%calc_initial_sat(n)
2044 topm = this%dis%top(m)
2045 botm = this%dis%bot(m)
2046 satm = this%calc_initial_sat(m)
2049 ihc = this%dis%con%ihc(jj)
2050 hyn = this%hy_eff(n, m, ihc, ipos=ii)
2051 hym = this%hy_eff(m, n, ihc, ipos=ii)
2052 if (this%ithickstartflag(n) == 0)
then
2055 hn = this%ic%strt(n)
2057 if (this%ithickstartflag(m) == 0)
then
2060 hm = this%ic%strt(m)
2068 csat =
vcond(1, 1, 1, 1, 0, 1, 1,
done, &
2074 this%dis%con%hwva(jj))
2078 fawidth = this%dis%con%hwva(jj)
2079 csat =
hcond(1, 1, 1, 1, 0, &
2083 hn, hm, satn, satm, hyn, hym, &
2086 this%dis%con%cl1(jj), &
2087 this%dis%con%cl2(jj), &
2090 this%condsat(jj) = csat
2107 integer(I4B),
intent(in) :: n
2112 if (this%ibound(n) /= 0 .and. this%ithickstartflag(n) /= 0)
then
2113 call this%thksat(n, this%ic%strt(n), satn)
2128 integer(I4B),
intent(in) :: kiter
2129 real(DP),
intent(inout),
dimension(:) :: hnew
2131 integer(I4B) :: n, m, ii, ihc
2132 real(DP) :: ttop, bbot, thick
2133 integer(I4B) :: ncnvrt, ihdcnv
2134 character(len=30),
dimension(5) :: nodcnvrt
2135 character(len=30) :: nodestr
2136 character(len=3),
dimension(5) :: acnvrt
2137 character(len=LINELENGTH) :: errmsg
2138 integer(I4B) :: irewet
2140 character(len=*),
parameter :: fmtnct = &
2141 "(1X,/1X,'Negative cell thickness at (layer,row,col)', &
2143 character(len=*),
parameter :: fmttopbot = &
2144 &
"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
2145 character(len=*),
parameter :: fmttopbotthk = &
2146 &
"(1X,'Top elevation, bottom elevation, thickness:',1P,3G13.5)"
2147 character(len=*),
parameter :: fmtdrychd = &
2148 &
"(1X,/1X,'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED')"
2149 character(len=*),
parameter :: fmtni = &
2150 &
"(1X,'CELLID=',a,' ITERATION=',I0,' TIME STEP=',I0,' STRESS PERIOD=',I0)"
2157 do n = 1, this%dis%nodes
2158 do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2159 m = this%dis%con%ja(ii)
2160 ihc = this%dis%con%ihc(this%dis%con%jas(ii))
2161 call this%rewet_check(kiter, n, hnew(m), this%ibound(m), ihc, hnew, &
2163 if (irewet == 1)
then
2164 call this%wdmsg(2, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2170 do n = 1, this%dis%nodes
2173 if (this%ibound(n) == 0) cycle
2174 if (this%icelltype(n) == 0) cycle
2177 bbot = this%dis%bot(n)
2178 ttop = this%dis%top(n)
2179 if (bbot > ttop)
then
2180 write (errmsg, fmtnct) n
2182 write (errmsg, fmttopbot) ttop, bbot
2188 if (this%icelltype(n) /= 0)
then
2189 if (hnew(n) < ttop) ttop = hnew(n)
2194 if (thick <= dzero)
then
2195 call this%wdmsg(1, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2197 if (this%ibound(n) < 0)
then
2198 errmsg =
'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED'
2200 write (errmsg, fmttopbotthk) ttop, bbot, thick
2202 call this%dis%noder_to_string(n, nodestr)
2203 write (errmsg, fmtni) trim(adjustl(nodestr)), kiter,
kstp,
kper
2212 call this%wdmsg(0, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2215 do n = 1, this%dis%nodes
2216 if (this%ibound(n) == 30000) this%ibound(n) = 1
2230 subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet)
2233 integer(I4B),
intent(in) :: kiter
2234 integer(I4B),
intent(in) :: node
2235 real(DP),
intent(in) :: hm
2236 integer(I4B),
intent(in) :: ibdm
2237 integer(I4B),
intent(in) :: ihc
2238 real(DP),
intent(inout),
dimension(:) :: hnew
2239 integer(I4B),
intent(out) :: irewet
2241 integer(I4B) :: itflg
2242 real(DP) :: wd, awd, turnon, bbot
2247 if (this%irewet > 0)
then
2248 itflg = mod(kiter, this%iwetit)
2249 if (itflg == 0)
then
2250 if (this%ibound(node) == 0 .and. this%wetdry(node) /=
dzero)
then
2253 bbot = this%dis%bot(node)
2254 wd = this%wetdry(node)
2256 if (wd < 0) awd = -wd
2264 if (ibdm > 0 .and. hm >= turnon) irewet = 1
2266 if (wd >
dzero)
then
2269 if (ibdm > 0 .and. hm >= turnon) irewet = 1
2273 if (irewet == 1)
then
2276 if (this%ihdwet == 0)
then
2277 hnew(node) = bbot + this%wetfct * (hm - bbot)
2279 hnew(node) = bbot + this%wetfct * awd
2281 this%ibound(node) = 30000
2299 integer(I4B),
intent(in) :: icode
2300 integer(I4B),
intent(inout) :: ncnvrt
2301 character(len=30),
dimension(5),
intent(inout) :: nodcnvrt
2302 character(len=3),
dimension(5),
intent(inout) :: acnvrt
2303 integer(I4B),
intent(inout) :: ihdcnv
2304 integer(I4B),
intent(in) :: kiter
2305 integer(I4B),
intent(in) :: n
2309 character(len=*),
parameter :: fmtcnvtn = &
2310 "(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I0, &
2311 &' STEP=',I0,' PERIOD=',I0,' (NODE or LRC)')"
2312 character(len=*),
parameter :: fmtnode =
"(1X,3X,5(A4, A20))"
2317 call this%dis%noder_to_string(n, nodcnvrt(ncnvrt))
2318 if (icode == 1)
then
2319 acnvrt(ncnvrt) =
'DRY'
2321 acnvrt(ncnvrt) =
'WET'
2327 if (ncnvrt == 5 .or. (icode == 0 .and. ncnvrt > 0))
then
2328 if (ihdcnv == 0)
write (this%iout, fmtcnvtn) kiter,
kstp,
kper
2330 write (this%iout, fmtnode) &
2331 (acnvrt(l), trim(adjustl(nodcnvrt(l))), l=1, ncnvrt)
2349 function hy_eff(this, n, m, ihc, ipos, vg)
result(hy)
2354 integer(I4B),
intent(in) :: n
2355 integer(I4B),
intent(in) :: m
2356 integer(I4B),
intent(in) :: ihc
2357 integer(I4B),
intent(in),
optional :: ipos
2358 real(dp),
dimension(3),
intent(in),
optional :: vg
2360 integer(I4B) :: iipos
2361 real(dp) :: hy11, hy22, hy33
2362 real(dp) :: ang1, ang2, ang3
2363 real(dp) :: vg1, vg2, vg3
2367 if (
present(ipos)) iipos = ipos
2381 if (this%iangle2 > 0)
then
2382 if (
present(vg))
then
2387 call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2389 ang1 = this%angle1(n)
2390 ang2 = this%angle2(n)
2392 if (this%iangle3 > 0) ang3 = this%angle3(n)
2393 hy =
hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2401 if (this%ik22 > 0)
then
2402 if (
present(vg))
then
2407 call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2412 if (this%iangle1 > 0)
then
2413 ang1 = this%angle1(n)
2414 if (this%iangle2 > 0)
then
2415 ang2 = this%angle2(n)
2416 if (this%iangle3 > 0) ang3 = this%angle3(n)
2419 hy =
hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2436 real(DP),
intent(in),
dimension(:) :: flowja
2440 integer(I4B) :: ipos
2441 integer(I4B) :: isympos
2469 real(DP),
allocatable,
dimension(:) :: vi
2470 real(DP),
allocatable,
dimension(:) :: di
2471 real(DP),
allocatable,
dimension(:) :: viz
2472 real(DP),
allocatable,
dimension(:) :: diz
2473 real(DP),
allocatable,
dimension(:) :: nix
2474 real(DP),
allocatable,
dimension(:) :: niy
2475 real(DP),
allocatable,
dimension(:) :: wix
2476 real(DP),
allocatable,
dimension(:) :: wiy
2477 real(DP),
allocatable,
dimension(:) :: wiz
2478 real(DP),
allocatable,
dimension(:) :: bix
2479 real(DP),
allocatable,
dimension(:) :: biy
2480 logical :: nozee = .true.
2483 if (this%icalcspdis /= 0 .and. this%dis%con%ianglex == 0)
then
2484 call store_error(
'Error. ANGLDEGX not provided in '// &
2485 'discretization file. ANGLDEGX required for '// &
2486 'calculation of specific discharge.', terminate=.true.)
2491 do n = 1, this%dis%nodes
2494 ic = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
2497 do m = 1, this%nedges
2498 if (this%nodedge(m) == n)
then
2504 if (ic > nc) nc = ic
2521 do n = 1, this%dis%nodes
2533 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2534 m = this%dis%con%ja(ipos)
2535 isympos = this%dis%con%jas(ipos)
2536 ihc = this%dis%con%ihc(isympos)
2537 area = this%dis%con%hwva(isympos)
2543 call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2544 ihc, xc, yc, zc, dltot)
2545 cl1 = this%dis%con%cl1(isympos)
2546 cl2 = this%dis%con%cl2(isympos)
2548 cl1 = this%dis%con%cl2(isympos)
2549 cl2 = this%dis%con%cl1(isympos)
2551 ooclsum =
done / (cl1 + cl2)
2552 diz(iz) = dltot * cl1 * ooclsum
2560 dz =
thksatnm(this%ibound(n), this%ibound(m), &
2561 this%icelltype(n), this%icelltype(m), &
2562 this%inewton, ihc, &
2563 this%hnew(n), this%hnew(m), this%sat(n), this%sat(m), &
2564 this%dis%top(n), this%dis%top(m), this%dis%bot(n), &
2567 call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
2568 call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2569 ihc, xc, yc, zc, dltot)
2570 cl1 = this%dis%con%cl1(isympos)
2571 cl2 = this%dis%con%cl2(isympos)
2573 cl1 = this%dis%con%cl2(isympos)
2574 cl2 = this%dis%con%cl1(isympos)
2576 ooclsum =
done / (cl1 + cl2)
2579 di(ic) = dltot * cl1 * ooclsum
2580 if (area >
dzero)
then
2581 vi(ic) = flowja(ipos) / area
2590 do m = 1, this%nedges
2591 if (this%nodedge(m) == n)
then
2594 ihc = this%ihcedge(m)
2595 area = this%propsedge(2, m)
2598 viz(iz) = this%propsedge(1, m) / area
2599 diz(iz) = this%propsedge(5, m)
2602 nix(ic) = -this%propsedge(3, m)
2603 niy(ic) = -this%propsedge(4, m)
2604 di(ic) = this%propsedge(5, m)
2605 if (area >
dzero)
then
2606 vi(ic) = this%propsedge(1, m) / area
2624 dsumz = dsumz + diz(iz)
2626 denom = (ncz -
done)
2628 dsumz = dsumz +
dem10 * dsumz
2630 if (dsumz >
dzero) wiz(iz) =
done - diz(iz) / dsumz
2632 wiz(iz) = wiz(iz) / denom
2640 vz = vz + wiz(iz) * viz(iz)
2649 wix(ic) = di(ic) * abs(nix(ic))
2650 wiy(ic) = di(ic) * abs(niy(ic))
2651 dsumx = dsumx + wix(ic)
2652 dsumy = dsumy + wiy(ic)
2659 dsumx = dsumx +
dem10 * dsumx
2660 dsumy = dsumy +
dem10 * dsumy
2662 wix(ic) = (dsumx - wix(ic)) * abs(nix(ic))
2663 wiy(ic) = (dsumy - wiy(ic)) * abs(niy(ic))
2670 bix(ic) = wix(ic) * sign(
done, nix(ic))
2671 biy(ic) = wiy(ic) * sign(
done, niy(ic))
2672 dsumx = dsumx + wix(ic) * abs(nix(ic))
2673 dsumy = dsumy + wiy(ic) * abs(niy(ic))
2680 bix(ic) = bix(ic) * dsumx
2681 biy(ic) = biy(ic) * dsumy
2682 axy = axy + bix(ic) * niy(ic)
2683 ayx = ayx + biy(ic) * nix(ic)
2696 vx = vx + (bix(ic) - axy * biy(ic)) * vi(ic)
2697 vy = vy + (biy(ic) - ayx * bix(ic)) * vi(ic)
2699 denom =
done - axy * ayx
2700 if (denom /=
dzero)
then
2705 this%spdis(1, n) = vx
2706 this%spdis(2, n) = vy
2707 this%spdis(3, n) = vz
2731 integer(I4B),
intent(in) :: ibinun
2733 character(len=16) :: text
2734 character(len=16),
dimension(3) :: auxtxt
2736 integer(I4B) :: naux
2739 text =
' DATA-SPDIS'
2741 auxtxt(:) = [
' qx',
' qy',
' qz']
2742 call this%dis%record_srcdst_list_header(text, this%name_model, &
2743 this%packName, this%name_model, &
2744 this%packName, naux, auxtxt, ibinun, &
2745 this%dis%nodes, this%iout)
2748 do n = 1, this%dis%nodes
2749 call this%dis%record_mf6_list_entry(ibinun, n, n,
dzero, naux, &
2762 integer(I4B),
intent(in) :: ibinun
2764 character(len=16) :: text
2765 character(len=16),
dimension(1) :: auxtxt
2766 real(DP),
dimension(1) :: a
2768 integer(I4B) :: naux
2773 auxtxt(:) = [
' sat']
2774 call this%dis%record_srcdst_list_header(text, this%name_model, &
2775 this%packName, this%name_model, &
2776 this%packName, naux, auxtxt, ibinun, &
2777 this%dis%nodes, this%iout)
2780 do n = 1, this%dis%nodes
2782 call this%dis%record_mf6_list_entry(ibinun, n, n,
dzero, naux, a)
2797 integer(I4B),
intent(in) :: nedges
2799 this%nedges = this%nedges + nedges
2811 integer(I4B),
intent(in) :: nodedge
2812 integer(I4B),
intent(in) :: ihcedge
2813 real(DP),
intent(in) :: q
2814 real(DP),
intent(in) :: area
2815 real(DP),
intent(in) :: nx
2816 real(DP),
intent(in) :: ny
2817 real(DP),
intent(in) :: distance
2819 integer(I4B) :: lastedge
2821 this%lastedge = this%lastedge + 1
2822 lastedge = this%lastedge
2823 this%nodedge(lastedge) = nodedge
2824 this%ihcedge(lastedge) = ihcedge
2825 this%propsedge(1, lastedge) = q
2826 this%propsedge(2, lastedge) = area
2827 this%propsedge(3, lastedge) = nx
2828 this%propsedge(4, lastedge) = ny
2829 this%propsedge(5, lastedge) = distance
2833 if (this%lastedge == this%nedges) this%lastedge = 0
2848 real(dp) :: satthickness
2850 satthickness =
thksatnm(this%ibound(n), &
2852 this%icelltype(n), &
2853 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.
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 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 memorylist_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...