21 integer(I4B),
pointer :: maxhfb => null()
22 integer(I4B),
pointer :: nhfb => null()
23 integer(I4B),
dimension(:),
pointer,
contiguous :: noden => null()
24 integer(I4B),
dimension(:),
pointer,
contiguous :: nodem => null()
25 integer(I4B),
dimension(:),
pointer,
contiguous :: idxloc => null()
26 real(dp),
dimension(:),
pointer,
contiguous :: hydchr => null()
27 real(dp),
dimension(:),
pointer,
contiguous :: csatsav => null()
28 real(dp),
dimension(:),
pointer,
contiguous :: condsav => null()
31 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
32 integer(I4B),
dimension(:),
pointer,
contiguous :: icelltype => null()
33 integer(I4B),
dimension(:),
pointer,
contiguous :: ihc => null()
34 integer(I4B),
dimension(:),
pointer,
contiguous :: ia => null()
35 integer(I4B),
dimension(:),
pointer,
contiguous :: ja => null()
36 integer(I4B),
dimension(:),
pointer,
contiguous :: jas => null()
37 integer(I4B),
dimension(:),
pointer,
contiguous :: isym => null()
38 real(dp),
dimension(:),
pointer,
contiguous :: condsat => null()
39 real(dp),
dimension(:),
pointer,
contiguous :: top => null()
40 real(dp),
dimension(:),
pointer,
contiguous :: bot => null()
41 real(dp),
dimension(:),
pointer,
contiguous :: hwva => null()
42 real(dp),
dimension(:),
pointer,
contiguous :: hnew => null()
45 integer(I4B),
pointer :: ivsc => null()
69 subroutine hfb_cr(hfbobj, name_model, inunit, iout)
72 character(len=*),
intent(in) :: name_model
73 integer(I4B),
intent(in) :: inunit
74 integer(I4B),
intent(in) :: iout
80 call hfbobj%set_names(1, name_model,
'HFB',
'HFB')
83 call hfbobj%allocate_scalars()
86 hfbobj%inunit = inunit
90 call hfbobj%parser%Initialize(hfbobj%inunit, hfbobj%iout)
98 subroutine hfb_ar(this, ibound, xt3d, dis, invsc, vsc)
104 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
107 integer(I4B),
pointer :: invsc
110 character(len=*),
parameter :: fmtheader = &
111 "(1x, /1x, 'HFB -- HORIZONTAL FLOW BARRIER PACKAGE, VERSION 8, ', &
112 &'4/24/2015 INPUT READ FROM UNIT ', i4, //)"
115 write (this%iout, fmtheader) this%inunit
119 this%ibound => ibound
122 call mem_setptr(this%icelltype,
'ICELLTYPE', &
135 call this%read_options()
136 call this%read_dimensions()
137 call this%allocate_arrays()
145 write (this%iout,
'(/1x,a,a)')
'Viscosity active in ', &
146 trim(this%filtyp)//
' Package calculations: '//trim(adjustl(this%packName))
163 character(len=LINELENGTH) :: line, errmsg
167 character(len=*),
parameter :: fmtblkerr = &
168 &
"('Error. Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
169 character(len=*),
parameter :: fmtlsp = &
170 &
"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
174 if (this%ionper <
kper)
then
177 call this%parser%GetBlock(
'PERIOD', isfound, ierr, &
178 supportopenclose=.true., &
179 blockrequired=.false.)
183 call this%read_check_ionper()
189 this%ionper =
nper + 1
192 call this%parser%GetCurrentLine(line)
193 write (errmsg, fmtblkerr) adjustl(trim(line))
195 call this%parser%StoreErrorUnit()
200 if (this%ionper ==
kper)
then
201 call this%condsat_reset()
202 call this%read_data()
203 call this%condsat_modify()
205 write (this%iout, fmtlsp)
'HFB'
220 subroutine hfb_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
225 integer(I4B) :: kiter
227 integer(I4B),
intent(in),
dimension(:) :: idxglo
228 real(DP),
intent(inout),
dimension(:) :: rhs
229 real(DP),
intent(inout),
dimension(:) :: hnew
231 integer(I4B) :: nodes, nja
232 integer(I4B) :: ihfb, n, m
234 integer(I4B) :: idiag, isymcon
235 integer(I4B) :: ixt3d
236 real(DP) :: cond, condhfb, aterm
237 real(DP) :: fawidth, faheight
238 real(DP) :: topn, topm, botn, botm
239 real(DP) :: viscratio
243 nodes = this%dis%nodes
244 nja = this%dis%con%nja
245 if (
associated(this%xt3d%ixt3d))
then
246 ixt3d = this%xt3d%ixt3d
253 do ihfb = 1, this%nhfb
254 n = min(this%noden(ihfb), this%nodem(ihfb))
255 m = max(this%noden(ihfb), this%nodem(ihfb))
257 if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle
259 if (this%ivsc /= 0)
then
260 call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio)
263 if (this%hydchr(ihfb) >
dzero)
then
264 if (this%inewton == 0)
then
265 ipos = this%idxloc(ihfb)
270 if (this%icelltype(n) == 1)
then
271 if (hnew(n) < topn) topn = hnew(n)
273 if (this%icelltype(m) == 1)
then
274 if (hnew(m) < topm) topm = hnew(m)
276 if (this%ihc(this%jas(ipos)) == 2)
then
277 faheight = min(topn, topm) - max(botn, botm)
279 faheight =
dhalf * ((topn - botn) + (topm - botm))
281 fawidth = this%hwva(this%jas(ipos))
282 condhfb = this%hydchr(ihfb) * viscratio * &
285 condhfb = this%hydchr(ihfb) * viscratio
288 condhfb = this%hydchr(ihfb)
291 call this%xt3d%xt3d_fhfb(kiter, nodes, nja, matrix_sln, idxglo, &
292 rhs, hnew, n, m, condhfb)
298 if (this%inewton == 0)
then
299 do ihfb = 1, this%nhfb
300 ipos = this%idxloc(ihfb)
301 aterm = matrix_sln%get_value_pos(idxglo(ipos))
304 if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle
306 if (this%ivsc /= 0)
then
307 call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio)
310 if (this%icelltype(n) == 1 .or. this%icelltype(m) == 1 .or. &
318 if (this%icelltype(n) == 1)
then
319 if (hnew(n) < topn) topn = hnew(n)
321 if (this%icelltype(m) == 1)
then
322 if (hnew(m) < topm) topm = hnew(m)
324 if (this%ihc(this%jas(ipos)) == 2)
then
325 faheight = min(topn, topm) - max(botn, botm)
327 faheight =
dhalf * ((topn - botn) + (topm - botm))
329 if (this%hydchr(ihfb) >
dzero)
then
330 fawidth = this%hwva(this%jas(ipos))
331 condhfb = this%hydchr(ihfb) * viscratio * &
333 cond = aterm * condhfb / (aterm + condhfb)
335 cond = -aterm * this%hydchr(ihfb)
339 this%condsav(ihfb) = cond
343 call matrix_sln%add_value_pos(idxglo(idiag), aterm - cond)
344 call matrix_sln%set_value_pos(idxglo(ipos), cond)
347 isymcon = this%isym(ipos)
349 call matrix_sln%add_value_pos(idxglo(idiag), aterm - cond)
350 call matrix_sln%set_value_pos(idxglo(isymcon), cond)
372 real(DP),
intent(inout),
dimension(:) :: hnew
373 real(DP),
intent(inout),
dimension(:) :: flowja
375 integer(I4B) :: ihfb, n, m
379 integer(I4B) :: ixt3d
381 real(DP) :: fawidth, faheight
382 real(DP) :: topn, topm, botn, botm
383 real(DP) :: viscratio
388 if (
associated(this%xt3d%ixt3d))
then
389 ixt3d = this%xt3d%ixt3d
396 do ihfb = 1, this%nhfb
397 n = min(this%noden(ihfb), this%nodem(ihfb))
398 m = max(this%noden(ihfb), this%nodem(ihfb))
400 if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle
402 if (this%ivsc /= 0)
then
403 call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio)
407 if (this%hydchr(ihfb) >
dzero)
then
408 if (this%inewton == 0)
then
409 ipos = this%idxloc(ihfb)
414 if (this%icelltype(n) == 1)
then
415 if (hnew(n) < topn) topn = hnew(n)
417 if (this%icelltype(m) == 1)
then
418 if (hnew(m) < topm) topm = hnew(m)
420 if (this%ihc(this%jas(ipos)) == 2)
then
421 faheight = min(topn, topm) - max(botn, botm)
423 faheight =
dhalf * ((topn - botn) + (topm - botm))
425 fawidth = this%hwva(this%jas(ipos))
426 condhfb = this%hydchr(ihfb) * viscratio * &
429 condhfb = this%hydchr(ihfb)
432 condhfb = this%hydchr(ihfb)
435 call this%xt3d%xt3d_flowjahfb(n, m, hnew, flowja, condhfb)
441 if (this%inewton == 0)
then
442 do ihfb = 1, this%nhfb
445 if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle
446 if (this%icelltype(n) == 1 .or. this%icelltype(m) == 1 .or. &
448 ipos = this%dis%con%getjaindex(n, m)
451 cond = this%condsav(ihfb)
452 qnm = cond * (hnew(m) - hnew(n))
454 ipos = this%dis%con%getjaindex(m, n)
481 if (this%inunit > 0)
then
491 call this%NumericalPackageType%da()
495 this%inewton => null()
496 this%ibound => null()
497 this%icelltype => null()
503 this%condsat => null()
522 call this%NumericalPackageType%allocate_scalars()
525 call mem_allocate(this%maxhfb,
'MAXHFB', this%memoryPath)
550 call mem_allocate(this%noden, this%maxhfb,
'NODEN', this%memoryPath)
551 call mem_allocate(this%nodem, this%maxhfb,
'NODEM', this%memoryPath)
552 call mem_allocate(this%hydchr, this%maxhfb,
'HYDCHR', this%memoryPath)
553 call mem_allocate(this%idxloc, this%maxhfb,
'IDXLOC', this%memoryPath)
554 call mem_allocate(this%csatsav, this%maxhfb,
'CSATSAV', this%memoryPath)
555 call mem_allocate(this%condsav, this%maxhfb,
'CONDSAV', this%memoryPath)
558 do ihfb = 1, this%maxhfb
559 this%idxloc(ihfb) = 0
575 character(len=LINELENGTH) :: errmsg, keyword
577 logical :: isfound, endOfBlock
580 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, &
581 supportopenclose=.true., blockrequired=.false.)
585 write (this%iout,
'(1x,a)')
'PROCESSING HFB OPTIONS'
587 call this%parser%GetNextLine(endofblock)
589 call this%parser%GetStringCaps(keyword)
590 select case (keyword)
593 write (this%iout,
'(4x,a)') &
594 'THE LIST OF HFBS WILL BE PRINTED.'
596 write (errmsg,
'(a,a)')
'Unknown HFB option: ', &
599 call this%parser%StoreErrorUnit()
602 write (this%iout,
'(1x,a)')
'END OF HFB OPTIONS'
617 character(len=LINELENGTH) :: errmsg, keyword
619 logical :: isfound, endOfBlock
622 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
623 supportopenclose=.true.)
627 write (this%iout,
'(/1x,a)')
'PROCESSING HFB DIMENSIONS'
629 call this%parser%GetNextLine(endofblock)
631 call this%parser%GetStringCaps(keyword)
632 select case (keyword)
634 this%maxhfb = this%parser%GetInteger()
635 write (this%iout,
'(4x,a,i7)')
'MAXHFB = ', this%maxhfb
637 write (errmsg,
'(a,a)') &
638 'Unknown HFB dimension: ', trim(keyword)
640 call this%parser%StoreErrorUnit()
644 write (this%iout,
'(1x,a)')
'END OF HFB DIMENSIONS'
646 call store_error(
'Required DIMENSIONS block not found.')
647 call this%parser%StoreErrorUnit()
651 if (this%maxhfb <= 0)
then
652 write (errmsg,
'(a)') &
653 'MAXHFB must be specified with value greater than zero.'
655 call this%parser%StoreErrorUnit()
677 character(len=LINELENGTH) :: nodenstr, nodemstr, cellidm, cellidn
678 integer(I4B) :: ihfb, nerr
679 logical :: endOfBlock
681 character(len=*),
parameter :: fmthfb =
"(i10, 2a10, 1(1pg15.6))"
683 write (this%iout,
'(//,1x,a)')
'READING HFB DATA'
684 if (this%iprpak > 0)
then
685 write (this%iout,
'(3a10, 1a15)')
'HFB NUM',
'CELL1',
'CELL2', &
694 call this%parser%GetNextLine(endofblock)
699 if (ihfb > this%maxhfb)
then
701 call this%parser%StoreErrorUnit()
703 call this%parser%GetCellid(this%dis%ndim, cellidn)
704 this%noden(ihfb) = this%dis%noder_from_cellid(cellidn, &
705 this%parser%iuactive, &
707 call this%parser%GetCellid(this%dis%ndim, cellidm)
708 this%nodem(ihfb) = this%dis%noder_from_cellid(cellidm, &
709 this%parser%iuactive, &
711 this%hydchr(ihfb) = this%parser%GetDouble()
714 if (this%iprpak /= 0)
then
715 call this%dis%noder_to_string(this%noden(ihfb), nodenstr)
716 call this%dis%noder_to_string(this%nodem(ihfb), nodemstr)
717 write (this%iout, fmthfb) ihfb, trim(adjustl(nodenstr)), &
718 trim(adjustl(nodemstr)), this%hydchr(ihfb)
727 call store_error(
'Errors encountered in HFB input file.')
728 call this%parser%StoreErrorUnit()
731 write (this%iout,
'(3x,i0,a,i0)') this%nhfb, &
732 ' HFBs READ FOR STRESS PERIOD ',
kper
733 call this%check_data()
734 write (this%iout,
'(1x,a)')
'END READING HFB DATA'
751 integer(I4B) :: ihfb, n, m
753 character(len=LINELENGTH) :: nodenstr, nodemstr
754 character(len=LINELENGTH) :: errmsg
757 character(len=*),
parameter :: fmterr =
"(1x, 'HFB no. ',i0, &
758 &' is between two unconnected cells: ', a, ' and ', a)"
759 character(len=*),
parameter :: fmtverr =
"(1x, 'HFB no. ',i0, &
760 &' is between two cells not horizontally connected: ', a, ' and ', a)"
762 do ihfb = 1, this%nhfb
766 do ipos = this%ia(n) + 1, this%ia(n + 1) - 1
767 if (m == this%ja(ipos))
then
769 this%idxloc(ihfb) = ipos
775 if (.not. found)
then
776 call this%dis%noder_to_string(n, nodenstr)
777 call this%dis%noder_to_string(m, nodemstr)
778 write (errmsg, fmterr) ihfb, trim(adjustl(nodenstr)), &
779 trim(adjustl(nodemstr))
784 ipos = this%idxloc(ihfb)
785 if (this%ihc(this%jas(ipos)) == 0)
then
786 call this%dis%noder_to_string(n, nodenstr)
787 call this%dis%noder_to_string(m, nodemstr)
788 write (errmsg, fmtverr) ihfb, trim(adjustl(nodenstr)), &
789 trim(adjustl(nodemstr))
813 do ihfb = 1, this%nhfb
814 ipos = this%idxloc(ihfb)
815 this%condsat(this%jas(ipos)) = this%csatsav(ihfb)
834 integer(I4B) :: ihfb, n, m
836 real(DP) :: cond, condhfb
837 real(DP) :: fawidth, faheight
838 real(DP) :: topn, topm, botn, botm
840 do ihfb = 1, this%nhfb
841 ipos = this%idxloc(ihfb)
842 cond = this%condsat(this%jas(ipos))
843 this%csatsav(ihfb) = cond
847 if (this%inewton == 1 .or. &
848 (this%icelltype(n) == 0 .and. this%icelltype(m) == 0))
then
855 if (this%ihc(this%jas(ipos)) == 2)
then
856 faheight = min(topn, topm) - max(botn, botm)
858 faheight =
dhalf * ((topn - botn) + (topm - botm))
860 if (this%hydchr(ihfb) >
dzero)
then
861 fawidth = this%hwva(this%jas(ipos))
862 condhfb = this%hydchr(ihfb) * &
864 cond = cond * condhfb / (cond + condhfb)
866 cond = -cond * this%hydchr(ihfb)
868 this%condsat(this%jas(ipos)) = cond
This module contains block parser methods.
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dzero
real constant zero
real(dp), parameter done
real constant 1
subroutine hfb_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
Fill matrix terms.
subroutine check_data(this)
Check for hfb's between two unconnected cells and write a warning.
subroutine hfb_da(this)
Deallocate memory.
subroutine hfb_cq(this, hnew, flowja)
flowja will automatically include the effects of the hfb for confined and newton cases when xt3d is n...
subroutine read_dimensions(this)
Read the dimensions for this package.
subroutine read_options(this)
Read a hfb options block.
subroutine allocate_scalars(this)
Allocate package scalars.
subroutine condsat_modify(this)
Modify condsat.
subroutine condsat_reset(this)
Reset condsat to its value prior to being modified by hfb's.
subroutine, public hfb_cr(hfbobj, name_model, inunit, iout)
Create a new hfb object.
subroutine hfb_rp(this)
Check for new HFB stress period data.
subroutine allocate_arrays(this)
Allocate package arrays.
subroutine read_data(this)
Read HFB period block.
subroutine hfb_ar(this, ibound, xt3d, dis, invsc, vsc)
Allocate and read.
This module defines variable data types.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
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_unit(iunit, terminate)
Store the file unit number.
integer(i4b), pointer, public kper
current stress period number
integer(i4b), pointer, public nper
number of stress period