22 character(len=LENFTYPE) ::
ftype =
'EVT'
23 character(len=LENPACKAGENAME) ::
text =
' EVT'
24 character(len=LENPACKAGENAME) ::
texta =
' EVTA'
28 logical,
pointer,
private :: segsdefined
29 logical,
pointer,
private :: fixed_cell
30 logical,
pointer,
private :: read_as_arrays
31 logical,
pointer,
private :: surfratespecified
33 integer(I4B),
pointer,
private :: nseg => null()
35 real(dp),
dimension(:),
pointer,
contiguous :: surface => null()
36 real(dp),
dimension(:),
pointer,
contiguous :: rate => null()
37 real(dp),
dimension(:),
pointer,
contiguous :: depth => null()
38 real(dp),
dimension(:, :),
pointer,
contiguous :: pxdp => null()
39 real(dp),
dimension(:, :),
pointer,
contiguous :: petm => null()
40 real(dp),
dimension(:),
pointer,
contiguous :: petm0 => null()
41 integer(I4B),
dimension(:),
pointer,
contiguous :: nodesontop => null()
68 subroutine evt_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
71 class(
bndtype),
pointer :: packobj
72 integer(I4B),
intent(in) :: id
73 integer(I4B),
intent(in) :: ibcnum
74 integer(I4B),
intent(in) :: inunit
75 integer(I4B),
intent(in) :: iout
76 character(len=*),
intent(in) :: namemodel
77 character(len=*),
intent(in) :: pakname
78 character(len=*),
intent(in) :: mempath
80 type(
evttype),
pointer :: evtobj
87 call packobj%set_names(ibcnum, namemodel, pakname,
ftype, mempath)
91 call evtobj%evt_allocate_scalars()
94 call packobj%pack_initialize()
96 packobj%inunit = inunit
99 packobj%ibcnum = ibcnum
112 class(
evttype),
intent(inout) :: this
115 call this%BndExtType%allocate_scalars()
121 allocate (this%segsdefined)
122 allocate (this%fixed_cell)
123 allocate (this%read_as_arrays)
124 allocate (this%surfratespecified)
128 this%segsdefined = .true.
129 this%fixed_cell = .false.
130 this%read_as_arrays = .false.
131 this%surfratespecified = .false.
144 integer(I4B),
dimension(:),
pointer,
contiguous,
optional :: nodelist
145 real(DP),
dimension(:, :),
pointer,
contiguous,
optional :: auxvar
148 call this%BndExtType%allocate_arrays(nodelist, auxvar)
151 call mem_setptr(this%surface,
'SURFACE', this%input_mempath)
152 call mem_setptr(this%rate,
'RATE', this%input_mempath)
153 call mem_setptr(this%depth,
'DEPTH', this%input_mempath)
156 call mem_checkin(this%surface,
'SURFACE', this%memoryPath, &
157 'SURFACE', this%input_mempath)
158 call mem_checkin(this%rate,
'RATE', this%memoryPath, &
159 'RATE', this%input_mempath)
160 call mem_checkin(this%depth,
'DEPTH', this%memoryPath, &
161 'DEPTH', this%input_mempath)
164 if (.not. this%read_as_arrays)
then
165 if (this%nseg > 1)
then
168 call mem_setptr(this%pxdp,
'PXDP', this%input_mempath)
169 call mem_setptr(this%petm,
'PETM', this%input_mempath)
172 call mem_checkin(this%pxdp,
'PXDP', this%memoryPath, &
173 'PXDP', this%input_mempath)
174 call mem_checkin(this%petm,
'PETM', this%memoryPath, &
175 'PETM', this%input_mempath)
178 if (this%surfratespecified)
then
181 call mem_setptr(this%petm0,
'PETM0', this%input_mempath)
184 call mem_checkin(this%petm0,
'PETM0', this%memoryPath, &
185 'PETM0', this%input_mempath)
199 class(
evttype),
intent(inout) :: this
201 logical(LGP) :: found_fixed_cell = .false.
202 logical(LGP) :: found_readasarrays = .false.
203 logical(LGP) :: found_surfratespec = .false.
206 call this%BndExtType%source_options()
210 this%input_mempath, found_fixed_cell)
212 this%input_mempath, found_readasarrays)
214 this%input_mempath, found_surfratespec)
216 if (found_readasarrays)
then
217 if (this%dis%supports_layers())
then
220 errmsg =
'READASARRAYS option is not compatible with selected'// &
221 ' discretization type.'
227 if (found_readasarrays .and. found_surfratespec)
then
228 if (this%read_as_arrays)
then
229 errmsg =
'READASARRAYS option is not compatible with the'// &
230 ' SURF_RATE_SPECIFIED option.'
237 call this%evt_log_options(found_fixed_cell, found_readasarrays, &
253 class(
evttype),
intent(inout) :: this
254 logical(LGP),
intent(in) :: found_fixed_cell
255 logical(LGP),
intent(in) :: found_readasarrays
256 logical(LGP),
intent(in) :: found_surfratespec
258 character(len=*),
parameter :: fmtihact = &
259 &
"(4x, 'EVAPOTRANSPIRATION WILL BE APPLIED TO HIGHEST ACTIVE CELL.')"
260 character(len=*),
parameter :: fmtfixedcell = &
261 &
"(4x, 'EVAPOTRANSPIRATION WILL BE APPLIED TO SPECIFIED CELL.')"
262 character(len=*),
parameter :: fmtreadasarrays = &
263 &
"(4x, 'EVAPOTRANSPIRATION INPUT WILL BE READ AS ARRAYS.')"
264 character(len=*),
parameter :: fmtsrz = &
265 &
"(4x, 'ET RATE AT SURFACE WILL BE ZERO.')"
266 character(len=*),
parameter :: fmtsrs = &
267 &
"(4x, 'ET RATE AT SURFACE WILL BE AS SPECIFIED BY PETM0.')"
270 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text)) &
273 if (found_fixed_cell)
then
274 write (this%iout, fmtfixedcell)
277 if (found_readasarrays)
then
278 write (this%iout, fmtreadasarrays)
281 if (found_surfratespec)
then
282 write (this%iout, fmtsrs)
286 write (this%iout,
'(1x,a)') &
287 'END OF '//trim(adjustl(this%text))//
' OPTIONS'
299 class(
evttype),
intent(inout) :: this
301 logical(LGP) :: found_nseg = .false.
303 character(len=*),
parameter :: fmtnsegerr = &
304 &
"('Error: In EVT, NSEG must be > 0 but is specified as ',i0)"
309 if (this%read_as_arrays)
then
310 this%maxbound = this%dis%get_ncpl()
313 if (this%maxbound <= 0)
then
315 'MAXBOUND must be an integer greater than zero.'
323 call this%BndExtType%source_dimensions()
326 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%text)) &
330 call mem_set_value(this%nseg,
'NSEG', this%input_mempath, found_nseg)
334 write (this%iout,
'(4x,a,i0)')
'NSEG = ', this%nseg
336 if (this%nseg < 1)
then
337 write (
errmsg, fmtnsegerr) this%nseg
341 elseif (this%nseg > 1)
then
343 if (this%read_as_arrays)
then
344 errmsg =
'In the EVT package, NSEG cannot be greater than 1'// &
345 ' when READASARRAYS is used.'
354 write (this%iout,
'(1x,a)') &
355 'END OF '//trim(adjustl(this%text))//
' DIMENSIONS'
361 call this%define_listlabel()
373 class(
evttype),
intent(inout) :: this
375 if (this%read_as_arrays)
then
376 call this%default_nodelist()
391 class(
evttype),
intent(inout) :: this
393 if (this%iper /=
kper)
return
395 if (this%read_as_arrays)
then
399 this%dis, this%input_mempath)
404 call this%BndExtType%bnd_rp()
407 if (this%nseg > 1)
then
408 call this%check_pxdp()
412 if (this%iprpak /= 0)
then
413 call this%write_list()
419 if (.not. this%fixed_cell)
call this%set_nodesontop()
433 class(
evttype),
intent(inout) :: this
438 integer(I4B) :: ierrmono
441 character(len=15) :: nodestr
443 character(len=*),
parameter :: fmtpxdp0 = &
444 &
"('PXDP must be between 0 and 1. Found ', G0, ' for cell ', A)"
445 character(len=*),
parameter :: fmtpxdp = &
446 &
"('PXDP is not monotonically increasing for cell ', A)"
450 do n = 1, this%nbound
451 node = this%nodelist(n)
454 segloop:
do i = 1, this%nseg
457 if (i < this%nseg)
then
458 pxdp2 = this%pxdp(i, n)
459 if (pxdp2 <=
dzero .or. pxdp2 >=
done)
then
460 call this%dis%noder_to_string(node, nodestr)
461 write (
errmsg, fmtpxdp0) pxdp2, trim(nodestr)
469 if (pxdp2 - pxdp1 <
dzero)
then
470 if (ierrmono == 0)
then
472 call this%dis%noder_to_string(node, nodestr)
473 write (
errmsg, fmtpxdp) trim(nodestr)
484 call this%parser%StoreErrorUnit()
495 class(
evttype),
intent(inout) :: this
500 if (.not.
associated(this%nodesontop))
then
501 allocate (this%nodesontop(this%maxbound))
505 do n = 1, this%nbound
506 this%nodesontop(n) = this%nodelist(n)
519 integer(I4B) :: i, iseg, node
520 integer(I4B) :: idxdepth, idxrate
521 real(DP) :: c, d, h, s, x
523 real(DP) :: petm1, petm2, pxdp1, pxdp2, thcof, trhs
526 if (this%nbound == 0)
return
529 do i = 1, this%nbound
532 if (this%fixed_cell)
then
533 node = this%nodelist(i)
535 node = this%nodesontop(i)
546 if (.not. this%fixed_cell)
then
547 if (this%ibound(node) == 0) &
548 call this%dis%highest_active(node, this%ibound)
549 this%nodelist(i) = node
557 if (this%ibound(node) > 0 .and. this%ibound(node) /=
iwetlake)
then
559 if (this%iauxmultcol > 0)
then
560 c = this%rate(i) * this%dis%get_area(node) * &
561 this%auxvar(this%iauxmultcol, i)
563 c = this%rate(i) * this%dis%get_area(node)
567 if (this%surfratespecified)
then
568 petm0 = this%petm0(i)
573 if (this%surfratespecified)
then
575 this%rhs(i) = this%rhs(i) + petm0 * c
578 this%rhs(i) = this%rhs(i) + c
586 if (this%nseg > 1)
then
599 if (this%surfratespecified)
then
610 segloop:
do iseg = 1, this%nseg
613 if (iseg < this%nseg)
then
615 idxdepth = idxdepth + 1
616 idxrate = idxrate + 1
618 pxdp2 = this%pxdp(idxdepth, i)
619 petm2 = this%petm(idxrate, i)
624 if (d <= pxdp2 * x)
then
635 thcof = -(petm1 - petm2) * c / ((pxdp2 - pxdp1) * x)
636 trhs = thcof * (s - pxdp1 * x) + petm1 * c
643 this%rhs(i) = this%rhs(i) + trhs
644 this%hcof(i) = this%hcof(i) + thcof
657 subroutine evt_fc(this, rhs, ia, idxglo, matrix_sln)
660 real(DP),
dimension(:),
intent(inout) :: rhs
661 integer(I4B),
dimension(:),
intent(in) :: ia
662 integer(I4B),
dimension(:),
intent(in) :: idxglo
665 integer(I4B) :: i, n, ipos
668 do i = 1, this%nbound
672 if (this%ibound(n) ==
iwetlake)
then
677 rhs(n) = rhs(n) + this%rhs(i)
679 call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))
695 if (
associated(this%nodesontop))
deallocate (this%nodesontop)
700 if (.not. this%read_as_arrays)
then
701 if (this%nseg > 1)
then
706 if (this%surfratespecified)
then
713 deallocate (this%segsdefined)
714 deallocate (this%fixed_cell)
715 deallocate (this%read_as_arrays)
716 deallocate (this%surfratespecified)
719 call this%BndExtType%bnd_da()
730 class(
evttype),
intent(inout) :: this
732 integer(I4B) :: nsegm1, i
735 this%listlabel = trim(this%filtyp)//
' NO.'
736 if (this%dis%ndim == 3)
then
737 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
738 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'ROW'
739 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'COL'
740 elseif (this%dis%ndim == 2)
then
741 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
742 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'CELL2D'
744 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'NODE'
746 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'SURFACE'
747 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'MAX. RATE'
748 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'EXT. DEPTH'
751 nsegm1 = this%nseg - 1
754 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'PXDP'
757 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'PETM'
762 if (this%surfratespecified)
then
763 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'PETM0'
771 if (this%inamedbound == 1)
then
772 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'BOUNDARY NAME'
790 integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nodeu, noder, ipos
793 if (this%dis%ndim == 3)
then
794 nlay = this%dis%mshape(1)
795 nrow = this%dis%mshape(2)
796 ncol = this%dis%mshape(3)
797 elseif (this%dis%ndim == 2)
then
798 nlay = this%dis%mshape(1)
800 ncol = this%dis%mshape(2)
808 nodeu =
get_node(il, ir, ic, nlay, nrow, ncol)
809 noder = this%dis%get_nodenumber(nodeu, 0)
810 this%nodelist(ipos) = noder
816 this%nbound = ipos - 1
820 if (.not. this%fixed_cell)
call this%set_nodesontop()
852 call this%obs%StoreObsType(
'evt', .true., indx)
865 class(
evttype),
intent(inout) :: this
866 integer(I4B),
intent(in) :: col
867 integer(I4B),
intent(in) :: row
878 bndval = this%surface(row)
880 if (this%iauxmultcol > 0)
then
881 bndval = this%rate(row) * this%auxvar(this%iauxmultcol, row)
883 bndval = this%rate(row)
886 bndval = this%depth(row)
889 if (this%nseg > 1)
then
890 if (col < (3 + this%nseg))
then
892 bndval = this%pxdp(idx, row)
893 else if (col < (3 + (2 * this%nseg) - 1))
then
894 idx = col - (3 + this%nseg - 1)
895 bndval = this%petm(idx, row)
896 else if (col == (3 + (2 * this%nseg) - 1))
then
897 if (this%surfratespecified)
then
899 bndval = this%petm0(row)
902 else if (this%surfratespecified)
then
905 bndval = this%petm0(row)
912 write (
errmsg,
'(a,i0,a)') &
913 'Programming error. EVT bound value requested column '&
914 &
'outside range of ncolbnd (', this%ncolbnd,
').'
938 integer(I4B),
dimension(:),
contiguous, &
939 pointer,
intent(inout) :: nodelist
941 character(len=*),
intent(in) :: input_mempath
942 integer(I4B),
intent(inout) :: nbound
943 integer(I4B),
intent(in) :: maxbound
945 character(len=24) :: aname =
' LAYER OR NODE INDEX'
947 integer(I4B),
dimension(:),
contiguous,
pointer :: ievt => null()
948 integer(I4B),
pointer :: inievt => null()
951 call mem_setptr(inievt,
'INIEVT', input_mempath)
954 if (inievt == 1)
then
961 call dis%nlarray_to_nodelist(ievt, nodelist, &
962 maxbound, nbound, aname)
This module contains block parser methods.
This module contains the extended boundary package.
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 lenpackagename
maximum length of the package name
integer(i4b), parameter iwetlake
integer constant for a dry lake
integer(i4b), parameter lenftype
maximum length of a package type (DIS, WEL, OC, etc.)
real(dp), parameter dzero
real constant zero
integer(i4b), parameter maxcharlen
maximum length of char string
real(dp), parameter done
real constant 1
subroutine evt_log_options(this, found_fixed_cell, found_readasarrays, found_surfratespec)
Source options specific to EvtType.
subroutine evt_source_dimensions(this)
Source the dimensions for this package.
subroutine evt_rp(this)
Read and Prepare.
subroutine set_nodesontop(this)
Store nodelist in nodesontop.
subroutine evt_cf(this)
Formulate the HCOF and RHS terms.
character(len=lenpackagename) text
subroutine default_nodelist(this)
Assign default nodelist when READASARRAYS is specified.
subroutine evt_da(this)
Deallocate.
subroutine evt_df_obs(this)
Store observation type supported by EVT package.
subroutine evt_allocate_arrays(this, nodelist, auxvar)
Allocate package arrays.
real(dp) function evt_bound_value(this, col, row)
Return requested boundary value.
subroutine evt_read_initial_attr(this)
Part of allocate and read.
subroutine evt_source_options(this)
Source options specific to EvtType.
subroutine evt_allocate_scalars(this)
Allocate package scalar members.
subroutine, public evt_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, mempath)
Create a new Evapotranspiration Segments Package and point pakobj to the new package.
character(len=lenftype) ftype
subroutine nodelist_update(nodelist, nbound, maxbound, dis, input_mempath)
Update the nodelist based on IEVT input.
subroutine evt_fc(this, rhs, ia, idxglo, matrix_sln)
Copy rhs and hcof into solution rhs and amat.
logical function evt_obs_supported(this)
Return true because EVT package supports observations.
character(len=lenpackagename) texta
subroutine evt_define_listlabel(this)
Define the list heading that is written to iout when PRINT_INPUT option is used.
subroutine check_pxdp(this)
Subroutine to check pxdp.
integer(i4b) function, public get_node(ilay, irow, icol, nlay, nrow, ncol)
Get node number, given layer, row, and column indices for a structured grid. If any argument is inval...
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 derived type ObsType.
subroutine, public defaultobsidprocessor(obsrv, dis, inunitobs, iout)
@ brief Process IDstring provided for each observation
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
integer(i4b), pointer, public kper
current stress period number
This class is used to store a single deferred-length character string. It was designed to work in an ...