25 real(dp),
dimension(:),
pointer :: conc => null()
26 integer(I4B),
dimension(:),
pointer :: icbund => null()
31 integer(I4B),
pointer :: ioutdense => null()
32 integer(I4B),
pointer :: iform => null()
33 integer(I4B),
pointer :: ireadelev => null()
34 integer(I4B),
pointer :: ireadconcbuy => null()
35 integer(I4B),
pointer :: iconcset => null()
36 real(dp),
pointer :: denseref => null()
37 real(dp),
dimension(:),
pointer,
contiguous :: dense => null()
38 real(dp),
dimension(:),
pointer,
contiguous :: concbuy => null()
39 real(dp),
dimension(:),
pointer,
contiguous :: elev => null()
40 integer(I4B),
dimension(:),
pointer :: ibound => null()
42 integer(I4B),
pointer :: nrhospecies => null()
43 real(dp),
dimension(:),
pointer,
contiguous :: drhodc => null()
44 real(dp),
dimension(:),
pointer,
contiguous :: crhoref => null()
45 real(dp),
dimension(:),
pointer,
contiguous :: ctemp => null()
46 character(len=LENMODELNAME),
dimension(:),
allocatable :: cmodelname
47 character(len=LENAUXNAME),
dimension(:),
allocatable :: cauxspeciesname
81 function calcdens(denseref, drhodc, crhoref, conc)
result(dense)
83 real(dp),
intent(in) :: denseref
84 real(dp),
dimension(:),
intent(in) :: drhodc
85 real(dp),
dimension(:),
intent(in) :: crhoref
86 real(dp),
dimension(:),
intent(in) :: conc
90 integer(I4B) :: nrhospec
93 nrhospec =
size(drhodc)
96 dense = dense + drhodc(i) * (conc(i) - crhoref(i))
105 subroutine buy_cr(buyobj, name_model, inunit, iout)
108 character(len=*),
intent(in) :: name_model
109 integer(I4B),
intent(in) :: inunit
110 integer(I4B),
intent(in) :: iout
116 call buyobj%set_names(1, name_model,
'BUY',
'BUY')
119 call buyobj%allocate_scalars()
122 buyobj%inunit = inunit
126 call buyobj%parser%Initialize(buyobj%inunit, buyobj%iout)
141 character(len=*),
parameter :: fmtbuy = &
142 "(1x,/1x,'BUY -- Buoyancy Package, Version 1, 5/16/2018', &
143 &' input read from unit ', i0, //)"
146 write (this%iout, fmtbuy) this%inunit
151 if (.not.
present(buy_input))
then
154 call this%read_options()
157 call this%read_dimensions()
160 call this%set_options(buy_input)
161 this%nrhospecies = buy_input%nrhospecies
165 call this%allocate_arrays(dis%nodes)
167 if (.not.
present(buy_input))
then
170 call this%read_packagedata()
173 call this%set_packagedata(buy_input)
186 integer(I4B),
dimension(:),
pointer :: ibound
190 this%ibound => ibound
193 if (this%npf%ixt3d /= 0)
then
194 call store_error(
'Error in model '//trim(this%name_model)// &
195 '. The XT3D option cannot be used with the BUY Package.')
196 call this%parser%StoreErrorUnit()
200 call this%buy_calcelev()
218 class(
bndtype),
pointer :: packobj
219 real(DP),
intent(in),
dimension(:) :: hnew
222 select case (packobj%filtyp)
226 select type (packobj)
228 call packobj%lak_activate_density()
234 select type (packobj)
236 call packobj%sfr_activate_density()
242 select type (packobj)
244 call packobj%maw_activate_density()
264 character(len=LINELENGTH) :: errmsg
267 character(len=*),
parameter :: fmtc = &
268 "('Buoyancy Package does not have a concentration set &
269 &for species ',i0,'. One or more model names may be specified &
270 &incorrectly in the PACKAGEDATA block or a gwf-gwt exchange may need &
275 do i = 1, this%nrhospecies
276 if (.not.
associated(this%modelconc(i)%conc))
then
277 write (errmsg, fmtc) i
282 call this%parser%StoreErrorUnit()
297 call this%buy_calcdens()
308 integer(I4B) :: kiter
311 if (this%ireadelev == 0)
then
312 if (this%iform == 1 .or. this%iform == 2)
then
313 call this%buy_calcelev()
328 class(
bndtype),
pointer :: packobj
329 real(DP),
intent(in),
dimension(:) :: hnew
332 integer(I4B) :: n, locdense, locelev
333 integer(I4B),
dimension(:),
allocatable :: locconc
337 if (this%iform == 0)
return
342 allocate (locconc(this%nrhospecies))
346 do n = 1, packobj%naux
347 if (packobj%auxname(n) ==
'DENSITY')
then
349 else if (packobj%auxname(n) ==
'ELEVATION')
then
355 do i = 1, this%nrhospecies
357 do j = 1, packobj%naux
358 if (this%cauxspeciesname(i) == packobj%auxname(j))
then
363 if (locconc(i) == 0)
then
371 select case (packobj%filtyp)
375 call buy_cf_ghb(packobj, hnew, this%dense, this%elev, this%denseref, &
376 locelev, locdense, locconc, this%drhodc, this%crhoref, &
377 this%ctemp, this%iform)
381 call buy_cf_riv(packobj, hnew, this%dense, this%elev, this%denseref, &
382 locelev, locdense, locconc, this%drhodc, this%crhoref, &
383 this%ctemp, this%iform)
387 call buy_cf_drn(packobj, hnew, this%dense, this%denseref)
391 call buy_cf_lak(packobj, hnew, this%dense, this%elev, this%denseref, &
392 locdense, locconc, this%drhodc, this%crhoref, &
393 this%ctemp, this%iform)
397 call buy_cf_sfr(packobj, hnew, this%dense, this%elev, this%denseref, &
398 locdense, locconc, this%drhodc, this%crhoref, &
399 this%ctemp, this%iform)
403 call buy_cf_maw(packobj, hnew, this%dense, this%elev, this%denseref, &
404 locdense, locconc, this%drhodc, this%crhoref, &
405 this%ctemp, this%iform)
425 ctemp, auxvar)
result(densebnd)
427 integer(I4B),
intent(in) :: n
428 integer(I4B),
intent(in) :: locdense
429 integer(I4B),
dimension(:),
intent(in) :: locconc
430 real(dp),
intent(in) :: denseref
431 real(dp),
dimension(:),
intent(in) :: drhodc
432 real(dp),
dimension(:),
intent(in) :: crhoref
433 real(dp),
dimension(:),
intent(inout) :: ctemp
434 real(dp),
dimension(:, :),
intent(in) :: auxvar
441 if (locdense > 0)
then
443 densebnd = auxvar(locdense, n)
444 else if (locconc(1) > 0)
then
446 do i = 1,
size(locconc)
448 if (locconc(i) > 0)
then
449 ctemp(i) = auxvar(locconc(i), n)
452 densebnd =
calcdens(denseref, drhodc, crhoref, ctemp)
464 subroutine buy_cf_ghb(packobj, hnew, dense, elev, denseref, locelev, &
465 locdense, locconc, drhodc, crhoref, ctemp, &
470 class(
bndtype),
pointer :: packobj
472 real(DP),
intent(in),
dimension(:) :: hnew
473 real(DP),
intent(in),
dimension(:) :: dense
474 real(DP),
intent(in),
dimension(:) :: elev
475 real(DP),
intent(in) :: denseref
476 integer(I4B),
intent(in) :: locelev
477 integer(I4B),
intent(in) :: locdense
478 integer(I4B),
dimension(:),
intent(in) :: locconc
479 real(DP),
dimension(:),
intent(in) :: drhodc
480 real(DP),
dimension(:),
intent(in) :: crhoref
481 real(DP),
dimension(:),
intent(inout) :: ctemp
482 integer(I4B),
intent(in) :: iform
490 real(DP) :: hcofterm, rhsterm
493 select type (packobj)
495 do n = 1, packobj%nbound
496 node = packobj%nodelist(n)
497 if (packobj%ibound(node) <= 0) cycle
501 drhodc, crhoref, ctemp, packobj%auxvar)
505 if (locelev > 0) elevghb = packobj%auxvar(locelev, n)
508 hghb = packobj%bhead(n)
509 cond = packobj%cond(n)
513 elevghb, elev(node), hghb, hnew(node), &
514 cond, iform, rhsterm, hcofterm)
515 packobj%hcof(n) = packobj%hcof(n) + hcofterm
516 packobj%rhs(n) = packobj%rhs(n) - rhsterm
528 elevghb, elevnode, hghb, hnode, &
529 cond, iform, rhsterm, hcofterm)
531 real(DP),
intent(in) :: denseref
532 real(DP),
intent(in) :: denseghb
533 real(DP),
intent(in) :: densenode
534 real(DP),
intent(in) :: elevghb
535 real(DP),
intent(in) :: elevnode
536 real(DP),
intent(in) :: hghb
537 real(DP),
intent(in) :: hnode
538 real(DP),
intent(in) :: cond
539 integer(I4B),
intent(in) :: iform
540 real(DP),
intent(inout) :: rhsterm
541 real(DP),
intent(inout) :: hcofterm
544 real(DP) :: avgdense, avgelev
547 avgdense =
dhalf * denseghb +
dhalf * densenode
549 t1 = avgdense / denseref -
done
550 t2 = (denseghb - densenode) / denseref
553 hcofterm = -cond * t1
557 hcofterm = hcofterm +
dhalf * cond * t2
561 rhsterm = cond * t1 * hghb
562 rhsterm = rhsterm - cond * t2 * avgelev
563 rhsterm = rhsterm +
dhalf * cond * t2 * hghb
567 rhsterm = rhsterm +
dhalf * cond * t2 * hnode
576 subroutine buy_cf_riv(packobj, hnew, dense, elev, denseref, locelev, &
577 locdense, locconc, drhodc, crhoref, ctemp, &
582 class(
bndtype),
pointer :: packobj
584 real(DP),
intent(in),
dimension(:) :: hnew
585 real(DP),
intent(in),
dimension(:) :: dense
586 real(DP),
intent(in),
dimension(:) :: elev
587 real(DP),
intent(in) :: denseref
588 integer(I4B),
intent(in) :: locelev
589 integer(I4B),
intent(in) :: locdense
590 integer(I4B),
dimension(:),
intent(in) :: locconc
591 real(DP),
dimension(:),
intent(in) :: drhodc
592 real(DP),
dimension(:),
intent(in) :: crhoref
593 real(DP),
dimension(:),
intent(inout) :: ctemp
594 integer(I4B),
intent(in) :: iform
607 select type (packobj)
609 do n = 1, packobj%nbound
610 node = packobj%nodelist(n)
611 if (packobj%ibound(node) <= 0) cycle
615 drhodc, crhoref, ctemp, packobj%auxvar)
619 if (locelev > 0) elevriv = packobj%auxvar(locelev, n)
622 hriv = packobj%stage(n)
623 cond = packobj%cond(n)
624 rbot = packobj%rbot(n)
627 if (hnew(node) > rbot)
then
631 elevriv, elev(node), hriv, hnew(node), &
632 cond, iform, rhsterm, hcofterm)
635 rhsterm = cond * (denseriv / denseref -
done) * (hriv - rbot)
639 packobj%hcof(n) = packobj%hcof(n) + hcofterm
640 packobj%rhs(n) = packobj%rhs(n) - rhsterm
654 class(
bndtype),
pointer :: packobj
656 real(DP),
intent(in),
dimension(:) :: hnew
657 real(DP),
intent(in),
dimension(:) :: dense
658 real(DP),
intent(in) :: denseref
669 select type (packobj)
671 do n = 1, packobj%nbound
672 node = packobj%nodelist(n)
673 if (packobj%ibound(node) <= 0) cycle
675 hbnd = packobj%elev(n)
676 cond = packobj%cond(n)
677 if (hnew(node) > hbnd)
then
678 hcofterm = -cond * (rho / denseref -
done)
679 rhsterm = hcofterm * hbnd
680 packobj%hcof(n) = packobj%hcof(n) + hcofterm
681 packobj%rhs(n) = packobj%rhs(n) + rhsterm
694 subroutine buy_cf_lak(packobj, hnew, dense, elev, denseref, locdense, &
695 locconc, drhodc, crhoref, ctemp, iform)
699 class(
bndtype),
pointer :: packobj
701 real(DP),
intent(in),
dimension(:) :: hnew
702 real(DP),
intent(in),
dimension(:) :: dense
703 real(DP),
intent(in),
dimension(:) :: elev
704 real(DP),
intent(in) :: denseref
705 integer(I4B),
intent(in) :: locdense
706 integer(I4B),
dimension(:),
intent(in) :: locconc
707 real(DP),
dimension(:),
intent(in) :: drhodc
708 real(DP),
dimension(:),
intent(in) :: crhoref
709 real(DP),
dimension(:),
intent(inout) :: ctemp
710 integer(I4B),
intent(in) :: iform
718 select type (packobj)
720 do n = 1, packobj%nbound
723 node = packobj%nodelist(n)
724 if (packobj%ibound(node) <= 0) cycle
728 drhodc, crhoref, ctemp, packobj%auxvar)
731 packobj%denseterms(1, n) = denselak / denseref
734 packobj%denseterms(2, n) = dense(node) / denseref
737 packobj%denseterms(3, n) = elev(node)
750 subroutine buy_cf_sfr(packobj, hnew, dense, elev, denseref, locdense, &
751 locconc, drhodc, crhoref, ctemp, iform)
755 class(
bndtype),
pointer :: packobj
757 real(DP),
intent(in),
dimension(:) :: hnew
758 real(DP),
intent(in),
dimension(:) :: dense
759 real(DP),
intent(in),
dimension(:) :: elev
760 real(DP),
intent(in) :: denseref
761 integer(I4B),
intent(in) :: locdense
762 integer(I4B),
dimension(:),
intent(in) :: locconc
763 real(DP),
dimension(:),
intent(in) :: drhodc
764 real(DP),
dimension(:),
intent(in) :: crhoref
765 real(DP),
dimension(:),
intent(inout) :: ctemp
766 integer(I4B),
intent(in) :: iform
774 select type (packobj)
776 do n = 1, packobj%nbound
779 node = packobj%nodelist(n)
780 if (packobj%ibound(node) <= 0) cycle
784 drhodc, crhoref, ctemp, packobj%auxvar)
787 packobj%denseterms(1, n) = densesfr / denseref
790 packobj%denseterms(2, n) = dense(node) / denseref
793 packobj%denseterms(3, n) = elev(node)
806 subroutine buy_cf_maw(packobj, hnew, dense, elev, denseref, locdense, &
807 locconc, drhodc, crhoref, ctemp, iform)
811 class(
bndtype),
pointer :: packobj
813 real(DP),
intent(in),
dimension(:) :: hnew
814 real(DP),
intent(in),
dimension(:) :: dense
815 real(DP),
intent(in),
dimension(:) :: elev
816 real(DP),
intent(in) :: denseref
817 integer(I4B),
intent(in) :: locdense
818 integer(I4B),
dimension(:),
intent(in) :: locconc
819 real(DP),
dimension(:),
intent(in) :: drhodc
820 real(DP),
dimension(:),
intent(in) :: crhoref
821 real(DP),
dimension(:),
intent(inout) :: ctemp
822 integer(I4B),
intent(in) :: iform
830 select type (packobj)
832 do n = 1, packobj%nbound
835 node = packobj%nodelist(n)
836 if (packobj%ibound(node) <= 0) cycle
840 drhodc, crhoref, ctemp, packobj%auxvar)
843 packobj%denseterms(1, n) = densemaw / denseref
846 packobj%denseterms(2, n) = dense(node) / denseref
849 packobj%denseterms(3, n) = elev(node)
860 subroutine buy_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
863 integer(I4B) :: kiter
865 integer(I4B),
intent(in),
dimension(:) :: idxglo
866 real(DP),
dimension(:),
intent(inout) :: rhs
867 real(DP),
intent(inout),
dimension(:) :: hnew
869 integer(I4B) :: n, m, ipos, idiag
870 real(DP) :: rhsterm, amatnn, amatnm
877 do n = 1, this%dis%nodes
878 if (this%ibound(n) == 0) cycle
879 idiag = this%dis%con%ia(n)
880 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
881 m = this%dis%con%ja(ipos)
882 if (this%ibound(m) == 0) cycle
883 if (this%iform == 0)
then
884 call this%calcbuy(n, m, ipos, hnew(n), hnew(m), rhsterm)
886 call this%calchhterms(n, m, ipos, hnew(n), hnew(m), rhsterm, &
891 rhs(n) = rhs(n) - rhsterm
892 call matrix_sln%add_value_pos(idxglo(idiag), -amatnn)
893 call matrix_sln%add_value_pos(idxglo(ipos), amatnm)
906 integer(I4B),
intent(in) :: idvfl
908 character(len=1) :: cdatafmp =
' ', editdesc =
' '
909 integer(I4B) :: ibinun
910 integer(I4B) :: iprint
911 integer(I4B) :: nvaluesp
912 integer(I4B) :: nwidthp
916 if (this%ioutdense /= 0)
then
921 if (idvfl == 0) ibinun = 0
924 if (ibinun /= 0)
then
929 if (this%ioutdense /= 0)
then
930 ibinun = this%ioutdense
931 call this%dis%record_array(this%dense, this%iout, iprint, ibinun, &
932 ' DENSITY', cdatafmp, nvaluesp, &
933 nwidthp, editdesc, dinact)
946 real(DP),
intent(in),
dimension(:) :: hnew
947 real(DP),
intent(inout),
dimension(:) :: flowja
948 integer(I4B) :: n, m, ipos
950 real(DP) :: rhsterm, amatnn, amatnm
953 do n = 1, this%dis%nodes
954 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
955 m = this%dis%con%ja(ipos)
957 if (this%iform == 0)
then
959 call this%calcbuy(n, m, ipos, hnew(n), hnew(m), deltaq)
962 call this%calchhterms(n, m, ipos, hnew(n), hnew(m), rhsterm, &
964 deltaq = amatnm * hnew(m) - amatnn * hnew(n) + rhsterm
966 flowja(ipos) = flowja(ipos) + deltaq
967 flowja(this%dis%con%isym(ipos)) = flowja(this%dis%con%isym(ipos)) - &
983 if (this%inunit > 0)
then
990 deallocate (this%cmodelname)
991 deallocate (this%cauxspeciesname)
992 deallocate (this%modelconc)
1006 call this%NumericalPackageType%da()
1018 character(len=LINELENGTH) :: errmsg, keyword
1019 integer(I4B) :: ierr
1020 logical :: isfound, endOfBlock
1024 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
1025 supportopenclose=.true.)
1029 write (this%iout,
'(/1x,a)')
'Processing BUY DIMENSIONS block'
1031 call this%parser%GetNextLine(endofblock)
1032 if (endofblock)
exit
1033 call this%parser%GetStringCaps(keyword)
1034 select case (keyword)
1035 case (
'NRHOSPECIES')
1036 this%nrhospecies = this%parser%GetInteger()
1037 write (this%iout,
'(4x,a,i0)')
'NRHOSPECIES = ', this%nrhospecies
1039 write (errmsg,
'(a,a)') &
1040 'Unknown BUY dimension: ', trim(keyword)
1042 call this%parser%StoreErrorUnit()
1045 write (this%iout,
'(1x,a)')
'End of BUY DIMENSIONS block'
1047 call store_error(
'Required BUY DIMENSIONS block not found.')
1048 call this%parser%StoreErrorUnit()
1052 if (this%nrhospecies < 1)
then
1053 call store_error(
'NRHOSPECIES must be greater than zero.')
1054 call this%parser%StoreErrorUnit()
1067 character(len=LINELENGTH) :: errmsg
1068 character(len=LINELENGTH) :: line
1069 integer(I4B) :: ierr
1070 integer(I4B) :: irhospec
1071 logical :: isfound, endOfBlock
1072 logical :: blockrequired
1073 integer(I4B),
dimension(:),
allocatable :: itemp
1074 character(len=10) :: c10
1075 character(len=16) :: c16
1077 character(len=*),
parameter :: fmterr = &
1078 "('Invalid value for IRHOSPEC (',i0,') detected in BUY Package. &
1079 &IRHOSPEC must be > 0 and <= NRHOSPECIES, and duplicate values &
1080 &are not allowed.')"
1083 allocate (itemp(this%nrhospecies))
1087 blockrequired = .true.
1088 call this%parser%GetBlock(
'PACKAGEDATA', isfound, ierr, &
1089 blockrequired=blockrequired, &
1090 supportopenclose=.true.)
1094 write (this%iout,
'(1x,a)')
'Processing BUY PACKAGEDATA block'
1096 call this%parser%GetNextLine(endofblock)
1097 if (endofblock)
exit
1098 irhospec = this%parser%GetInteger()
1099 if (irhospec < 1 .or. irhospec > this%nrhospecies)
then
1100 write (errmsg, fmterr) irhospec
1103 if (itemp(irhospec) /= 0)
then
1104 write (errmsg, fmterr) irhospec
1108 this%drhodc(irhospec) = this%parser%GetDouble()
1109 this%crhoref(irhospec) = this%parser%GetDouble()
1110 call this%parser%GetStringCaps(this%cmodelname(irhospec))
1111 call this%parser%GetStringCaps(this%cauxspeciesname(irhospec))
1113 write (this%iout,
'(1x,a)')
'End of BUY PACKAGEDATA block'
1115 call store_error(
'Required BUY PACKAGEDATA block not found.')
1116 call this%parser%StoreErrorUnit()
1121 call this%parser%StoreErrorUnit()
1125 write (this%iout,
'(/,a)')
'Summary of species information in BUY Package'
1126 write (this%iout,
'(1a11, 4a17)') &
1127 'SPECIES',
'DRHODC',
'CRHOREF',
'MODEL', &
1129 do irhospec = 1, this%nrhospecies
1130 write (c10,
'(i0)') irhospec
1131 line =
' '//adjustr(c10)
1132 write (c16,
'(g15.6)') this%drhodc(irhospec)
1133 line = trim(line)//
' '//adjustr(c16)
1134 write (c16,
'(g15.6)') this%crhoref(irhospec)
1135 line = trim(line)//
' '//adjustr(c16)
1136 write (c16,
'(a)') this%cmodelname(irhospec)
1137 line = trim(line)//
' '//adjustr(c16)
1138 write (c16,
'(a)') this%cauxspeciesname(irhospec)
1139 line = trim(line)//
' '//adjustr(c16)
1140 write (this%iout,
'(a)') trim(line)
1157 integer(I4B) :: ispec
1159 do ispec = 1, this%nrhospecies
1160 this%drhodc(ispec) = input_data%drhodc(ispec)
1161 this%crhoref(ispec) = input_data%crhoref(ispec)
1162 this%cmodelname(ispec) = input_data%cmodelname(ispec)
1163 this%cauxspeciesname(ispec) = input_data%cauxspeciesname(ispec)
1177 integer(I4B),
intent(in) :: n
1178 integer(I4B),
intent(in) :: m
1179 integer(I4B),
intent(in) :: icon
1180 real(DP),
intent(in) :: hn
1181 real(DP),
intent(in) :: hm
1182 real(DP),
intent(inout) :: buy
1185 real(DP) :: densen, densem, cl1, cl2, avgdense, wt, elevn, elevm, &
1191 densen = this%dense(n)
1192 densem = this%dense(m)
1194 cl1 = this%dis%con%cl1(this%dis%con%jas(icon))
1195 cl2 = this%dis%con%cl2(this%dis%con%jas(icon))
1197 cl1 = this%dis%con%cl2(this%dis%con%jas(icon))
1198 cl2 = this%dis%con%cl1(this%dis%con%jas(icon))
1200 wt = cl1 / (cl1 + cl2)
1201 avgdense = wt * densen + (
done - wt) * densem
1204 if (this%ireadelev == 0)
then
1205 tp = this%dis%top(n)
1206 bt = this%dis%bot(n)
1207 elevn = bt +
dhalf * this%npf%sat(n) * (tp - bt)
1208 tp = this%dis%top(m)
1209 bt = this%dis%bot(m)
1210 elevm = bt +
dhalf * this%npf%sat(m) * (tp - bt)
1212 elevn = this%elev(n)
1213 elevm = this%elev(m)
1216 ihc = this%dis%con%ihc(this%dis%con%jas(icon))
1217 hyn = this%npf%hy_eff(n, m, ihc, ipos=icon)
1218 hym = this%npf%hy_eff(m, n, ihc, ipos=icon)
1221 if (this%dis%con%ihc(this%dis%con%jas(icon)) == 0)
then
1222 cond =
vcond(this%ibound(n), this%ibound(m), &
1223 this%npf%icelltype(n), this%npf%icelltype(m), &
1225 this%npf%ivarcv, this%npf%idewatcv, &
1226 this%npf%condsat(this%dis%con%jas(icon)), hn, hm, &
1228 this%npf%sat(n), this%npf%sat(m), &
1229 this%dis%top(n), this%dis%top(m), &
1230 this%dis%bot(n), this%dis%bot(m), &
1231 this%dis%con%hwva(this%dis%con%jas(icon)))
1233 cond =
hcond(this%ibound(n), this%ibound(m), &
1234 this%npf%icelltype(n), this%npf%icelltype(m), &
1236 this%dis%con%ihc(this%dis%con%jas(icon)), &
1237 this%npf%icellavg, &
1238 this%npf%condsat(this%dis%con%jas(icon)), &
1239 hn, hm, this%npf%sat(n), this%npf%sat(m), &
1241 this%dis%top(n), this%dis%top(m), &
1242 this%dis%bot(n), this%dis%bot(m), &
1243 this%dis%con%cl1(this%dis%con%jas(icon)), &
1244 this%dis%con%cl2(this%dis%con%jas(icon)), &
1245 this%dis%con%hwva(this%dis%con%jas(icon)))
1249 buy = cond * (avgdense - this%denseref) / this%denseref * (elevm - elevn)
1257 subroutine calchhterms(this, n, m, icon, hn, hm, rhsterm, amatnn, amatnm)
1262 integer(I4B),
intent(in) :: n
1263 integer(I4B),
intent(in) :: m
1264 integer(I4B),
intent(in) :: icon
1265 real(DP),
intent(in) :: hn
1266 real(DP),
intent(in) :: hm
1267 real(DP),
intent(inout) :: rhsterm
1268 real(DP),
intent(inout) :: amatnn
1269 real(DP),
intent(inout) :: amatnm
1272 real(DP) :: densen, densem, cl1, cl2, avgdense, wt, elevn, elevm, cond
1273 real(DP) :: rhonormn, rhonormm
1281 densen = this%dense(n)
1282 densem = this%dense(m)
1284 cl1 = this%dis%con%cl1(this%dis%con%jas(icon))
1285 cl2 = this%dis%con%cl2(this%dis%con%jas(icon))
1287 cl1 = this%dis%con%cl2(this%dis%con%jas(icon))
1288 cl2 = this%dis%con%cl1(this%dis%con%jas(icon))
1290 wt = cl1 / (cl1 + cl2)
1291 avgdense = wt * densen + (1.0 - wt) * densem
1294 elevn = this%elev(n)
1295 elevm = this%elev(m)
1296 elevnm = (
done - wt) * elevn + wt * elevm
1298 ihc = this%dis%con%ihc(this%dis%con%jas(icon))
1299 hyn = this%npf%hy_eff(n, m, ihc, ipos=icon)
1300 hym = this%npf%hy_eff(m, n, ihc, ipos=icon)
1304 cond =
vcond(this%ibound(n), this%ibound(m), &
1305 this%npf%icelltype(n), this%npf%icelltype(m), &
1307 this%npf%ivarcv, this%npf%idewatcv, &
1308 this%npf%condsat(this%dis%con%jas(icon)), hn, hm, &
1310 this%npf%sat(n), this%npf%sat(m), &
1311 this%dis%top(n), this%dis%top(m), &
1312 this%dis%bot(n), this%dis%bot(m), &
1313 this%dis%con%hwva(this%dis%con%jas(icon)))
1315 cond =
hcond(this%ibound(n), this%ibound(m), &
1316 this%npf%icelltype(n), this%npf%icelltype(m), &
1318 this%dis%con%ihc(this%dis%con%jas(icon)), &
1319 this%npf%icellavg, &
1320 this%npf%condsat(this%dis%con%jas(icon)), &
1321 hn, hm, this%npf%sat(n), this%npf%sat(m), &
1323 this%dis%top(n), this%dis%top(m), &
1324 this%dis%bot(n), this%dis%bot(m), &
1325 this%dis%con%cl1(this%dis%con%jas(icon)), &
1326 this%dis%con%cl2(this%dis%con%jas(icon)), &
1327 this%dis%con%hwva(this%dis%con%jas(icon)))
1331 rhonormn = densen / this%denseref
1332 rhonormm = densem / this%denseref
1333 rhoterm = wt * rhonormn + (
done - wt) * rhonormm
1334 amatnn = cond * (rhoterm -
done)
1336 rhsterm = -cond * (rhonormm - rhonormn) * elevnm
1337 if (this%iform == 1)
then
1339 hphi = (
done - wt) * hn + wt * hm
1340 rhsterm = rhsterm + cond * hphi * (rhonormm - rhonormn)
1343 amatnn = amatnn - cond * (
done - wt) * (rhonormm - rhonormn)
1344 amatnm = amatnm + cond * wt * (rhonormm - rhonormn)
1362 do n = 1, this%dis%nodes
1363 do i = 1, this%nrhospecies
1364 if (this%modelconc(i)%icbund(n) == 0)
then
1367 this%ctemp(i) = this%modelconc(i)%conc(n)
1370 this%dense(n) =
calcdens(this%denseref, this%drhodc, this%crhoref, &
1385 real(DP) :: tp, bt, frac
1388 do n = 1, this%dis%nodes
1389 tp = this%dis%top(n)
1390 bt = this%dis%bot(n)
1391 frac = this%npf%sat(n)
1392 this%elev(n) = bt +
dhalf * frac * (tp - bt)
1409 call this%NumericalPackageType%allocate_scalars()
1412 call mem_allocate(this%ioutdense,
'IOUTDENSE', this%memoryPath)
1413 call mem_allocate(this%iform,
'IFORM', this%memoryPath)
1414 call mem_allocate(this%ireadelev,
'IREADELEV', this%memoryPath)
1415 call mem_allocate(this%ireadconcbuy,
'IREADCONCBUY', this%memoryPath)
1416 call mem_allocate(this%iconcset,
'ICONCSET', this%memoryPath)
1417 call mem_allocate(this%denseref,
'DENSEREF', this%memoryPath)
1419 call mem_allocate(this%nrhospecies,
'NRHOSPECIES', this%memoryPath)
1425 this%ireadconcbuy = 0
1426 this%denseref = 1000.d0
1428 this%nrhospecies = 0
1443 integer(I4B),
intent(in) :: nodes
1448 call mem_allocate(this%dense, nodes,
'DENSE', this%memoryPath)
1449 call mem_allocate(this%concbuy, 0,
'CONCBUY', this%memoryPath)
1450 call mem_allocate(this%elev, nodes,
'ELEV', this%memoryPath)
1451 call mem_allocate(this%drhodc, this%nrhospecies,
'DRHODC', this%memoryPath)
1452 call mem_allocate(this%crhoref, this%nrhospecies,
'CRHOREF', this%memoryPath)
1453 call mem_allocate(this%ctemp, this%nrhospecies,
'CTEMP', this%memoryPath)
1454 allocate (this%cmodelname(this%nrhospecies))
1455 allocate (this%cauxspeciesname(this%nrhospecies))
1456 allocate (this%modelconc(this%nrhospecies))
1460 this%dense(i) = this%denseref
1461 this%elev(i) =
dzero
1465 do i = 1, this%nrhospecies
1466 this%drhodc(i) =
dzero
1467 this%crhoref(i) =
dzero
1468 this%ctemp(i) =
dzero
1469 this%cmodelname(i) =
''
1470 this%cauxspeciesname(i) =
''
1486 character(len=LINELENGTH) :: errmsg, keyword
1487 character(len=MAXCHARLEN) :: fname
1488 integer(I4B) :: ierr
1489 logical :: isfound, endOfBlock
1491 character(len=*),
parameter :: fmtfileout = &
1492 "(4x, 'BUY ', 1x, a, 1x, ' will be saved to file: ', &
1493 &a, /4x, 'opened on unit: ', I7)"
1496 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, &
1497 supportopenclose=.true., blockrequired=.false.)
1501 write (this%iout,
'(1x,a)')
'Processing BUY OPTIONS block'
1503 call this%parser%GetNextLine(endofblock)
1504 if (endofblock)
exit
1505 call this%parser%GetStringCaps(keyword)
1506 select case (keyword)
1507 case (
'HHFORMULATION_RHS')
1510 write (this%iout,
'(4x,a)') &
1511 'Hydraulic head formulation set to right-hand side'
1513 this%denseref = this%parser%GetDouble()
1514 write (this%iout,
'(4x,a,1pg15.6)') &
1515 'Reference density has been set to: ', &
1517 case (
'DEV_EFH_FORMULATION')
1518 call this%parser%DevOpt()
1521 write (this%iout,
'(4x,a)') &
1522 'Formulation set to equivalent freshwater head'
1524 call this%parser%GetStringCaps(keyword)
1525 if (keyword ==
'FILEOUT')
then
1526 call this%parser%GetString(fname)
1528 call openfile(this%ioutdense, this%iout, fname,
'DATA(BINARY)', &
1530 write (this%iout, fmtfileout) &
1531 'DENSITY', fname, this%ioutdense
1533 errmsg =
'Optional density keyword must be '// &
1534 'followed by FILEOUT'
1538 write (errmsg,
'(a,a)')
'Unknown BUY option: ', &
1541 call this%parser%StoreErrorUnit()
1544 write (this%iout,
'(1x,a)')
'End of BUY OPTIONS block'
1558 this%iform = input_data%iform
1559 this%denseref = input_data%denseref
1563 if (this%iform == 0 .or. this%iform == 1)
then
1580 character(len=LENMODELNAME),
intent(in) :: modelname
1581 real(DP),
dimension(:),
pointer :: conc
1582 integer(I4B),
dimension(:),
pointer :: icbund
1589 do i = 1, this%nrhospecies
1590 if (this%cmodelname(i) == modelname)
then
1591 this%modelconc(i)%conc => conc
1592 this%modelconc(i)%icbund => icbund
This module contains the base boundary package.
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
integer(i4b), parameter lenmodelname
maximum length of the model name
real(dp), parameter dhnoflo
real no flow constant
real(dp), parameter dhalf
real constant 1/2
integer(i4b), parameter lenauxname
maximum length of a aux variable
real(dp), parameter dzero
real constant zero
integer(i4b), parameter maxcharlen
maximum length of char string
real(dp), parameter done
real constant 1
real(dp) function calcdens(denseref, drhodc, crhoref, conc)
Generic function to calculate fluid density from concentration.
subroutine buy_cf(this, kiter)
Fill coefficients.
subroutine buy_cf_bnd(this, packobj, hnew)
Fill coefficients.
subroutine read_options(this)
Read package options.
subroutine buy_cf_lak(packobj, hnew, dense, elev, denseref, locdense, locconc, drhodc, crhoref, ctemp, iform)
Pass density information into lak package; density terms are calculated in the lake package as part o...
subroutine set_options(this, input_data)
Sets options as opposed to reading them from a file.
subroutine calchhterms(this, n, m, icon, hn, hm, rhsterm, amatnn, amatnm)
Calculate hydraulic head term for this connection.
subroutine buy_rp(this)
Check for new buy period data.
subroutine, public buy_cr(buyobj, name_model, inunit, iout)
Create a new BUY object.
subroutine buy_cf_maw(packobj, hnew, dense, elev, denseref, locdense, locconc, drhodc, crhoref, ctemp, iform)
Pass density information into maw package; density terms are calculated in the maw package as part of...
subroutine buy_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
Fill coefficients.
subroutine calc_ghb_hcof_rhs_terms(denseref, denseghb, densenode, elevghb, elevnode, hghb, hnode, cond, iform, rhsterm, hcofterm)
Calculate density hcof and rhs terms for ghb conditions.
subroutine allocate_scalars(this)
Allocate scalars used by the package.
subroutine buy_cq(this, hnew, flowja)
Add buy term to flowja.
subroutine buy_calcelev(this)
Calculate cell elevations to use in density flow equations.
subroutine set_concentration_pointer(this, modelname, conc, icbund)
Pass in a gwt model name, concentration array and ibound, and store a pointer to these in the BUY pac...
subroutine buy_cf_riv(packobj, hnew, dense, elev, denseref, locelev, locdense, locconc, drhodc, crhoref, ctemp, iform)
Fill riv coefficients.
subroutine set_packagedata(this, input_data)
Sets package data instead of reading from file.
subroutine buy_df(this, dis, buy_input)
Read options and package data, or set from argument.
subroutine read_packagedata(this)
Read PACKAGEDATA block.
subroutine buy_cf_sfr(packobj, hnew, dense, elev, denseref, locdense, locconc, drhodc, crhoref, ctemp, iform)
Pass density information into sfr package; density terms are calculated in the sfr package as part of...
subroutine buy_ot_dv(this, idvfl)
Save density array to binary file.
real(dp) function get_bnd_density(n, locdense, locconc, denseref, drhodc, crhoref, ctemp, auxvar)
Return the density of the boundary package using one of several different options in the following or...
subroutine buy_cf_ghb(packobj, hnew, dense, elev, denseref, locelev, locdense, locconc, drhodc, crhoref, ctemp, iform)
Fill ghb coefficients.
subroutine read_dimensions(this)
Read the dimensions for this package.
subroutine buy_ar(this, npf, ibound)
Allocate and Read.
subroutine buy_da(this)
Deallocate.
subroutine buy_cf_drn(packobj, hnew, dense, denseref)
Fill drn coefficients.
subroutine buy_ad(this)
Advance.
subroutine allocate_arrays(this, nodes)
Allocate arrays used by the package.
subroutine buy_ar_bnd(this, packobj, hnew)
Buoyancy ar_bnd routine to activate density in packages.
subroutine calcbuy(this, n, m, icon, hn, hm, buy)
Calculate buyancy term for this connection.
subroutine buy_calcdens(this)
calculate fluid density from concentration
This module contains stateless conductance functions.
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.
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.
This module defines variable data types.
This module contains the base numerical package type.
This module contains the SFR package methods.
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.
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number