MODFLOW 6  version 6.7.0.dev0
USGS Modular Hydrologic Model
gwe-uze.f90
Go to the documentation of this file.
1 ! -- Unsaturated Zone Flow Energy Transport Module
2 ! -- todo: save the uze temperature into the uze aux variable?
3 ! -- todo: calculate the uzf DENSE aux variable using temperature?
4 ! -- todo: GWD and GWD-TO-MVR do not seem to be included; prob in UZF?
5 !
6 ! UZF flows (flowbudptr) index var UZE term Transport Type
7 !---------------------------------------------------------------------------------
8 
9 ! -- terms from UZF that will be handled by parent APT Package
10 ! FLOW-JA-FACE idxbudfjf FLOW-JA-FACE cv2cv
11 ! GWF (aux FLOW-AREA) idxbudgwf GWF uzf2gwf
12 ! STORAGE (aux VOLUME) idxbudsto none used for water volumes
13 ! FROM-MVR idxbudfmvr FROM-MVR q * cext = this%qfrommvr(:)
14 ! AUXILIARY none none none
15 ! none none STORAGE (aux MASS)
16 ! none none AUXILIARY none
17 
18 ! -- terms from UZF that need to be handled here
19 ! INFILTRATION idxbudinfl INFILTRATION q < 0: q * cwell, else q * cuser
20 ! REJ-INF idxbudrinf REJ-INF q * cuze
21 ! UZET idxbuduzet UZET q * cet
22 ! REJ-INF-TO-MVR idxbudritm REJ-INF-TO-MVR q * cinfil?
23 ! THERMAL-EQUIL idxbudtheq THERMAL-EQUIL residual
24 
25 ! -- terms from UZF that should be skipped
26 
28 
29  use kindmodule, only: dp, i4b
31  use simmodule, only: store_error
32  use bndmodule, only: bndtype, getbndfromlist
33  use tspfmimodule, only: tspfmitype
34  use uzfmodule, only: uzftype
35  use observemodule, only: observetype
40 
41  implicit none
42 
43  public uze_create
44 
45  character(len=*), parameter :: ftype = 'UZE'
46  character(len=*), parameter :: flowtype = 'UZF'
47  character(len=16) :: text = ' UZE'
48 
49  type, extends(tspapttype) :: gweuzetype
50 
51  type(gweinputdatatype), pointer :: gwecommon => null() !< pointer to shared gwe data used by multiple packages but set in mst
52 
53  integer(I4B), pointer :: idxbudinfl => null() !< index of uzf infiltration terms in flowbudptr
54  integer(I4B), pointer :: idxbudrinf => null() !< index of rejected infiltration terms in flowbudptr
55  integer(I4B), pointer :: idxbuduzet => null() !< index of unsat et terms in flowbudptr
56  integer(I4B), pointer :: idxbudritm => null() !< index of rej infil to mover rate to mover terms in flowbudptr
57  integer(I4B), pointer :: idxbudtheq => null() !< index of thermal equilibration terms in flowbudptr
58 
59  real(dp), dimension(:), pointer, contiguous :: tempinfl => null() !< infiltration temperature
60  real(dp), dimension(:), pointer, contiguous :: tempuzet => null() !< unsat et temperature
61 
62  contains
63 
64  procedure :: bnd_da => uze_da
65  procedure :: allocate_scalars
66  procedure :: apt_allocate_arrays => uze_allocate_arrays
67  procedure :: find_apt_package => find_uze_package
68  procedure :: apt_fc_expanded => uze_fc_expanded
69  procedure :: pak_solve => uze_solve
70  procedure :: pak_get_nbudterms => uze_get_nbudterms
71  procedure :: pak_setup_budobj => uze_setup_budobj
72  procedure :: pak_fill_budobj => uze_fill_budobj
73  procedure :: uze_infl_term
74  procedure :: uze_rinf_term
75  procedure :: uze_uzet_term
76  procedure :: uze_ritm_term
77  procedure :: uze_theq_term
78  procedure :: pak_df_obs => uze_df_obs
79  procedure :: pak_rp_obs => uze_rp_obs
80  procedure :: pak_bd_obs => uze_bd_obs
81  procedure :: pak_set_stressperiod => uze_set_stressperiod
82  procedure :: apt_ad_chk => uze_ad_chk
83  procedure :: bnd_ac => uze_ac
84  procedure :: bnd_mc => uze_mc
85  procedure :: get_mvr_depvar
86  procedure, private :: area_error
87 
88  end type gweuzetype
89 
90 contains
91 
92  !> @brief Create a new UZE package
93  !<
94  subroutine uze_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
95  fmi, eqnsclfac, gwecommon, dvt, dvu, dvua)
96  ! -- dummy
97  class(bndtype), pointer :: packobj
98  integer(I4B), intent(in) :: id
99  integer(I4B), intent(in) :: ibcnum
100  integer(I4B), intent(in) :: inunit
101  integer(I4B), intent(in) :: iout
102  character(len=*), intent(in) :: namemodel
103  character(len=*), intent(in) :: pakname
104  type(tspfmitype), pointer :: fmi
105  real(dp), intent(in), pointer :: eqnsclfac !< governing equation scale factor
106  type(gweinputdatatype), intent(in), target :: gwecommon !< shared data container for use by multiple GWE packages
107  character(len=*), intent(in) :: dvt !< For GWE, set to "TEMPERATURE" in TspAptType
108  character(len=*), intent(in) :: dvu !< For GWE, set to "energy" in TspAptType
109  character(len=*), intent(in) :: dvua !< For GWE, set to "E" in TspAptType
110  ! -- local
111  type(gweuzetype), pointer :: uzeobj
112  !
113  ! -- Allocate the object and assign values to object variables
114  allocate (uzeobj)
115  packobj => uzeobj
116  !
117  ! -- Create name and memory path
118  call packobj%set_names(ibcnum, namemodel, pakname, ftype)
119  packobj%text = text
120  !
121  ! -- Allocate scalars
122  call uzeobj%allocate_scalars()
123  !
124  ! -- Initialize package
125  call packobj%pack_initialize()
126  !
127  packobj%inunit = inunit
128  packobj%iout = iout
129  packobj%id = id
130  packobj%ibcnum = ibcnum
131  packobj%ncolbnd = 1
132  packobj%iscloc = 1
133 
134  ! -- Store pointer to flow model interface. When the GwfGwt exchange is
135  ! created, it sets fmi%bndlist so that the GWT model has access to all
136  ! the flow packages
137  uzeobj%fmi => fmi
138  !
139  ! -- Store pointer to governing equation scale factor
140  uzeobj%eqnsclfac => eqnsclfac
141  !
142  ! -- Store pointer to shared data module for accessing cpw, rhow
143  ! for the budget calculations, and for accessing the latent heat of
144  ! vaporization
145  uzeobj%gwecommon => gwecommon
146  !
147  ! -- Set labels that will be used in generalized APT class
148  uzeobj%depvartype = dvt
149  uzeobj%depvarunit = dvu
150  uzeobj%depvarunitabbrev = dvua
151  end subroutine uze_create
152 
153  !> @brief Find corresponding uze package
154  !<
155  subroutine find_uze_package(this)
156  ! -- modules
158  use simvariablesmodule, only: errmsg
159  ! -- dummy
160  class(gweuzetype) :: this
161  ! -- local
162  class(bndtype), pointer :: packobj
163  integer(I4B) :: ip, icount
164  integer(I4B) :: nbudterm
165  logical :: found
166  !
167  ! -- Initialize found to false, and error later if flow package cannot
168  ! be found
169  found = .false.
170  !
171  ! -- If user is specifying flows in a binary budget file, then set up
172  ! the budget file reader, otherwise set a pointer to the flow package
173  ! budobj
174  if (this%fmi%flows_from_file) then
175  call this%fmi%set_aptbudobj_pointer(this%flowpackagename, this%flowbudptr)
176  if (associated(this%flowbudptr)) then
177  found = .true.
178  end if
179  !
180  else
181  if (associated(this%fmi%gwfbndlist)) then
182  ! -- Look through gwfbndlist for a flow package with the same name as
183  ! this transport package name
184  do ip = 1, this%fmi%gwfbndlist%Count()
185  packobj => getbndfromlist(this%fmi%gwfbndlist, ip)
186  if (packobj%packName == this%flowpackagename) then
187  found = .true.
188  !
189  ! -- Store BndType pointer to packobj, and then
190  ! use the select type to point to the budobj in flow package
191  this%flowpackagebnd => packobj
192  select type (packobj)
193  type is (uzftype)
194  this%flowbudptr => packobj%budobj
195  end select
196  end if
197  if (found) exit
198  end do
199  end if
200  end if
201  !
202  ! -- Error if flow package not found
203  if (.not. found) then
204  write (errmsg, '(a)') 'Could not find flow package with name '&
205  &//trim(adjustl(this%flowpackagename))//'.'
206  call store_error(errmsg)
207  call this%parser%StoreErrorUnit()
208  end if
209  !
210  ! -- Allocate space for idxbudssm, which indicates whether this is a
211  ! special budget term or one that is a general source and sink
212  nbudterm = this%flowbudptr%nbudterm
213  call mem_allocate(this%idxbudssm, nbudterm, 'IDXBUDSSM', this%memoryPath)
214  !
215  ! -- Process budget terms and identify special budget terms
216  write (this%iout, '(/, a, a)') &
217  'PROCESSING '//ftype//' INFORMATION FOR ', this%packName
218  write (this%iout, '(a)') ' IDENTIFYING FLOW TERMS IN '//flowtype//' PACKAGE'
219  write (this%iout, '(a, i0)') &
220  ' NUMBER OF '//flowtype//' = ', this%flowbudptr%ncv
221  icount = 1
222  do ip = 1, this%flowbudptr%nbudterm
223  select case (trim(adjustl(this%flowbudptr%budterm(ip)%flowtype)))
224  case ('FLOW-JA-FACE')
225  this%idxbudfjf = ip
226  this%idxbudssm(ip) = 0
227  case ('GWF')
228  this%idxbudgwf = ip
229  this%idxbudssm(ip) = 0
230  case ('STORAGE')
231  this%idxbudsto = ip
232  this%idxbudssm(ip) = 0
233  case ('INFILTRATION')
234  this%idxbudinfl = ip
235  this%idxbudssm(ip) = 0
236  case ('REJ-INF')
237  this%idxbudrinf = ip
238  this%idxbudssm(ip) = 0
239  case ('UZET')
240  this%idxbuduzet = ip
241  this%idxbudssm(ip) = 0
242  case ('REJ-INF-TO-MVR')
243  this%idxbudritm = ip
244  this%idxbudssm(ip) = 0
245  case ('TO-MVR')
246  this%idxbudtmvr = ip
247  this%idxbudssm(ip) = 0
248  case ('FROM-MVR')
249  this%idxbudfmvr = ip
250  this%idxbudssm(ip) = 0
251  case ('AUXILIARY')
252  this%idxbudaux = ip
253  this%idxbudssm(ip) = 0
254  case default
255  !
256  ! -- Set idxbudssm equal to a column index for where the temperatures
257  ! are stored in the tempbud(nbudssm, ncv) array
258  this%idxbudssm(ip) = icount
259  icount = icount + 1
260  end select
261  !
262  write (this%iout, '(a, i0, " = ", a,/, a, i0)') &
263  ' TERM ', ip, trim(adjustl(this%flowbudptr%budterm(ip)%flowtype)), &
264  ' MAX NO. OF ENTRIES = ', this%flowbudptr%budterm(ip)%maxlist
265  end do
266  write (this%iout, '(a, //)') 'DONE PROCESSING '//ftype//' INFORMATION'
267  !
268  ! -- Thermal equilibration term
269  this%idxbudtheq = this%flowbudptr%nbudterm + 1
270  end subroutine find_uze_package
271 
272  !> @brief Add package connection to matrix.
273  !!
274  !! Overrides apt_ac to fold the UZE heat balance terms into the row
275  !! corresponding to the host cell and enforce thermal equilibrium between
276  !! UZE and the GWE cell.
277  !<
278  subroutine uze_ac(this, moffset, sparse)
280  use sparsemodule, only: sparsematrix
281  ! -- dummy
282  class(gweuzetype), intent(inout) :: this
283  integer(I4B), intent(in) :: moffset !< current models starting position in global matrix
284  type(sparsematrix), intent(inout) :: sparse
285  ! -- local
286  integer(I4B) :: i, ii
287  integer(I4B) :: n !< index of a uze object within the complete list of uze objects
288  integer(I4B) :: jj !<
289  integer(I4B) :: nglo !< index of uze object in global matrix
290  integer(I4B) :: jglo !< host cell's position in global matrix for a uze object
291  integer(I4B) :: idxn !< used for cross-check
292  integer(I4B) :: idxjj !< used for cross-check
293  integer(I4B) :: idxnglo !< used for cross-check
294  integer(I4B) :: idxjglo !< used for cross-check
295  !
296  ! -- Add package rows to sparse
297  if (this%imatrows /= 0) then
298  !
299  ! -- Diagonal on the row assoc. with the uze feature
300  do n = 1, this%ncv
301  nglo = moffset + this%dis%nodes + this%ioffset + n
302  call sparse%addconnection(nglo, nglo, 1)
303  end do
304  !
305  ! -- Add uze-to-gwe connections. For uze, this particular do loop
306  ! is the same as its counterpart in apt_ac.
307  ! nlist: number of gwe cells with a connection to at least one uze object
308  do i = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
309  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(i) !< uze object position within uze object list
310  jj = this%flowbudptr%budterm(this%idxbudgwf)%id2(i) !< position of gwe cell to which uze feature is connected
311  nglo = moffset + this%dis%nodes + this%ioffset + n !< uze feature position
312  jglo = moffset + jj !< gwe cell position
313  call sparse%addconnection(nglo, jglo, 1)
314  call sparse%addconnection(jglo, nglo, 1)
315  end do
316  !
317  ! -- For uze, add feature-to-feature connections (i.e.,
318  ! vertically stacked UZ objects) to row corresponding
319  ! to the host cell. Terms added to the row associated with host
320  ! cell are added to columns associated with other uze features.
321  ! This approach deviates from the approach taken in apt_ac.
322  if (this%idxbudfjf /= 0) then
323  do i = 1, this%flowbudptr%budterm(this%idxbudfjf)%maxlist
324  n = this%flowbudptr%budterm(this%idxbudfjf)%id1(i) !< position of currently considered uze feature
325  jj = this%flowbudptr%budterm(this%idxbudfjf)%id2(i) !< position of connected uze feature
326  nglo = moffset + this%dis%nodes + this%ioffset + n !< global position of currently considered uze feature
327  jglo = moffset + this%dis%nodes + this%ioffset + jj !< global position of connected uze feature
328  ! -- If connected uze feature is upstream, find cell that hosts currently
329  ! considered uze feature and add connection to that cell's row
330  do ii = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist !< uze object id among uze objects
331  idxn = this%flowbudptr%budterm(this%idxbudgwf)%id1(ii) !< uze object position within uze object list
332  idxjj = this%flowbudptr%budterm(this%idxbudgwf)%id2(ii) !< position of gwe cell to which uze feature is connected
333  idxnglo = moffset + this%dis%nodes + this%ioffset + idxn !< uze feature global position
334  idxjglo = moffset + idxjj !< gwe cell global position
335  if (nglo == idxnglo) exit
336  end do
337  call sparse%addconnection(idxjglo, jglo, 1)
338  end do
339  end if
340  end if
341  end subroutine uze_ac
342 
343  !> @brief Map package connection to matrix
344  !<
345  subroutine uze_mc(this, moffset, matrix_sln)
346  use sparsemodule, only: sparsematrix
347  ! -- dummy
348  class(gweuzetype), intent(inout) :: this
349  integer(I4B), intent(in) :: moffset
350  class(matrixbasetype), pointer :: matrix_sln
351  ! -- local
352  integer(I4B) :: n, j, iglo, jglo
353  integer(I4B) :: idxn, idxj, idxiglo, idxjglo
354  integer(I4B) :: ipos, idxpos
355  !
356  ! -- Allocate memory for index arrays
357  call this%apt_allocate_index_arrays()
358  !
359  ! -- Store index positions
360  if (this%imatrows /= 0) then
361  !
362  ! -- Find the position of each connection in the global ia, ja structure
363  ! and store them in idxglo. idxglo allows this model to insert or
364  ! retrieve values into or from the global A matrix
365  ! apt rows
366  !
367  ! -- Feature diagonal in global matrix
368  do n = 1, this%ncv
369  iglo = moffset + this%dis%nodes + this%ioffset + n
370  this%idxpakdiag(n) = matrix_sln%get_position_diag(iglo)
371  end do
372  !
373  ! -- Cell to feature connection in global matrix
374  do ipos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
375  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(ipos) !< feature number
376  j = this%flowbudptr%budterm(this%idxbudgwf)%id2(ipos) !< cell number
377  iglo = moffset + this%dis%nodes + this%ioffset + n !< feature row index
378  jglo = j + moffset !< cell row index
379  ! -- Note that this is where idxlocnode is set for uze; it is set
380  ! to the host cell local row index rather than the feature local
381  ! row index
382  this%idxlocnode(n) = j ! kluge note: do we want to introduce a new array instead of co-opting idxlocnode???
383  ! -- For connection ipos in list of feature-cell connections,
384  ! global positions of feature-row diagonal and off-diagonal
385  ! corresponding to the cell
386  this%idxdglo(ipos) = matrix_sln%get_position_diag(iglo)
387  this%idxoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
388  end do
389  !
390  ! -- Feature to cell connection in global matrix
391  do ipos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
392  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(ipos) !< feature number
393  j = this%flowbudptr%budterm(this%idxbudgwf)%id2(ipos) !< cell number
394  iglo = j + moffset !< cell row index
395  jglo = moffset + this%dis%nodes + this%ioffset + n !< feature row index
396  ! -- For connection ipos in list of feature-cell connections,
397  ! global positions of cell-row diagonal and off-diagonal
398  ! corresponding to the feature
399  this%idxsymdglo(ipos) = matrix_sln%get_position_diag(iglo)
400  this%idxsymoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
401  end do
402  !
403  ! -- Feature to feature connection in global matrix
404  if (this%idxbudfjf /= 0) then
405  do ipos = 1, this%flowbudptr%budterm(this%idxbudfjf)%nlist
406  n = this%flowbudptr%budterm(this%idxbudfjf)%id1(ipos) !< number of currently considered uze feature
407  j = this%flowbudptr%budterm(this%idxbudfjf)%id2(ipos) !< number of connected uze feature
408  iglo = moffset + this%dis%nodes + this%ioffset + n !< global position of currently considered uze feature
409  jglo = moffset + this%dis%nodes + this%ioffset + j !< global position of connected uze feature
410  ! -- If connected uze feature is upstream, find cell that hosts currently
411  ! considered uze feature and map connection to that cell's row
412  do idxpos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
413  idxn = this%flowbudptr%budterm(this%idxbudgwf)%id1(idxpos) !< feature number
414  idxj = this%flowbudptr%budterm(this%idxbudgwf)%id2(idxpos) !< cell number)
415  idxjglo = moffset + this%dis%nodes + this%ioffset + idxn !< feature row index
416  idxiglo = moffset + idxj !< cell row index
417  if (idxjglo == iglo) exit
418  end do
419  ! -- For connection ipos in list of feature-feature connections,
420  ! global positions of host-cell-row entries corresponding to
421  ! (in the same columns as) the feature-id1-row diagonal and the
422  ! feature-id1-row off-diagonal corresponding to feature id2
423  this%idxfjfdglo(ipos) = matrix_sln%get_position_diag(idxiglo)
424  this%idxfjfoffdglo(ipos) = matrix_sln%get_position(idxiglo, jglo)
425  end do
426  end if
427  end if
428  end subroutine uze_mc
429 
430  !> @brief Add matrix terms related to UZE
431  !!
432  !! This will be called from TspAptType%apt_fc_expanded()
433  !! in order to add matrix terms specifically for this package
434  !<
435  subroutine uze_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
436  ! -- dummy
437  class(gweuzetype) :: this
438  real(DP), dimension(:), intent(inout) :: rhs
439  integer(I4B), dimension(:), intent(in) :: ia
440  integer(I4B), dimension(:), intent(in) :: idxglo
441  class(matrixbasetype), pointer :: matrix_sln
442  ! -- local
443  integer(I4B) :: j, n, n1, n2
444  integer(I4B) :: iloc
445  integer(I4B) :: iposd, iposoffd
446  integer(I4B) :: ipossymoffd
447  real(DP) :: cold
448  real(DP) :: qbnd
449  real(DP) :: omega
450  real(DP) :: rrate
451  real(DP) :: rhsval
452  real(DP) :: hcofval
453  !
454  ! -- Add infiltration contribution
455  ! uze does not put feature balance coefficients in the row
456  ! associated with the feature. Instead, these coefficients are
457  ! moved into the row associated with cell hosting the uze feature
458  if (this%idxbudinfl /= 0) then
459  do j = 1, this%flowbudptr%budterm(this%idxbudinfl)%nlist
460  call this%uze_infl_term(j, n1, n2, rrate, rhsval, hcofval)
461  iloc = this%idxlocnode(n1) !< for uze idxlocnode stores the host cell local row index
462  ipossymoffd = this%idxsymoffdglo(j)
463  call matrix_sln%add_value_pos(ipossymoffd, hcofval)
464  rhs(iloc) = rhs(iloc) + rhsval
465  end do
466  end if
467  !
468  ! -- Add rejected infiltration contribution
469  if (this%idxbudrinf /= 0) then
470  do j = 1, this%flowbudptr%budterm(this%idxbudrinf)%nlist
471  call this%uze_rinf_term(j, n1, n2, rrate, rhsval, hcofval)
472  iloc = this%idxlocnode(n1) ! for uze idxlocnode stores the host cell local row index
473  ipossymoffd = this%idxsymoffdglo(j)
474  call matrix_sln%add_value_pos(ipossymoffd, hcofval)
475  rhs(iloc) = rhs(iloc) + rhsval
476  end do
477  end if
478  !
479  ! -- Add unsaturated et contribution
480  if (this%idxbuduzet /= 0) then
481  do j = 1, this%flowbudptr%budterm(this%idxbuduzet)%nlist
482  call this%uze_uzet_term(j, n1, n2, rrate, rhsval, hcofval)
483  iloc = this%idxlocnode(n1) ! for uze idxlocnode stores the host cell local row index
484  ipossymoffd = this%idxsymoffdglo(j)
485  call matrix_sln%add_value_pos(ipossymoffd, hcofval)
486  rhs(iloc) = rhs(iloc) + rhsval
487  end do
488  end if
489  !
490  ! -- Add rejected infiltration to mover contribution
491  if (this%idxbudritm /= 0) then
492  do j = 1, this%flowbudptr%budterm(this%idxbudritm)%nlist
493  call this%uze_ritm_term(j, n1, n2, rrate, rhsval, hcofval)
494  iloc = this%idxlocnode(n1) ! for uze idxlocnode stores the host cell local row index
495  ipossymoffd = this%idxsymoffdglo(j)
496  call matrix_sln%add_value_pos(ipossymoffd, hcofval)
497  rhs(iloc) = rhs(iloc) + rhsval
498  end do
499  end if
500  !
501  ! -- For UZE, content of apt_fc_expanded placed here as the approach is to
502  ! completely override apt_fc_expanded() with what follows
503  !
504  ! -- Mass (or energy) storage in features
505  do n = 1, this%ncv
506  cold = this%xoldpak(n)
507  iloc = this%idxlocnode(n) ! for uze idxlocnode stores the host cell local row index
508  ipossymoffd = this%idxsymoffdglo(n)
509  call this%apt_stor_term(n, n1, n2, rrate, rhsval, hcofval)
510  call matrix_sln%add_value_pos(ipossymoffd, hcofval)
511  rhs(iloc) = rhs(iloc) + rhsval
512  end do
513  !
514  ! -- Add to mover contribution
515  if (this%idxbudtmvr /= 0) then
516  do j = 1, this%flowbudptr%budterm(this%idxbudtmvr)%nlist
517  call this%apt_tmvr_term(j, n1, n2, rrate, rhsval, hcofval)
518  iloc = this%idxlocnode(n1) ! for uze, idxlocnode stores the host cell local row index
519  ipossymoffd = this%idxsymoffdglo(j)
520  call matrix_sln%add_value_pos(ipossymoffd, hcofval)
521  rhs(iloc) = rhs(iloc) + rhsval
522  end do
523  end if
524  !
525  ! -- Add from mover contribution
526  if (this%idxbudfmvr /= 0) then
527  do n = 1, this%ncv
528  rhsval = this%qmfrommvr(n) ! kluge note: presumably already in terms of energy
529  iloc = this%idxlocnode(n) ! for uze idxlocnode stores the host cell local row index
530  rhs(iloc) = rhs(iloc) - rhsval
531  end do
532  end if
533  !
534  ! -- Go through each apt-gwf connection
535  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
536  !
537  ! -- Set n to feature number and process if active feature
538  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
539  if (this%iboundpak(n) /= 0) then
540  !
541  ! -- This code altered from its counterpart appearing in apt; this equates
542  ! uze temperature to cell temperature using the feature's row
543  iposd = this%idxdglo(j)
544  iposoffd = this%idxoffdglo(j)
545  call matrix_sln%add_value_pos(iposd, done)
546  call matrix_sln%add_value_pos(iposoffd, -done)
547  end if
548  end do
549  !
550  ! -- Go through each apt-apt connection
551  if (this%idxbudfjf /= 0) then
552  do j = 1, this%flowbudptr%budterm(this%idxbudfjf)%nlist
553  n1 = this%flowbudptr%budterm(this%idxbudfjf)%id1(j)
554  n2 = this%flowbudptr%budterm(this%idxbudfjf)%id2(j)
555  qbnd = this%flowbudptr%budterm(this%idxbudfjf)%flow(j)
556  if (qbnd <= dzero) then
557  omega = done
558  else
559  omega = dzero
560  end if
561  iposd = this%idxfjfdglo(j) !< position of feature-id1 column in feature id1's host-cell row
562  iposoffd = this%idxfjfoffdglo(j) !< position of feature-id2 column in feature id1's host-cell row
563  call matrix_sln%add_value_pos(iposd, omega * qbnd * this%eqnsclfac)
564  call matrix_sln%add_value_pos(iposoffd, &
565  (done - omega) * qbnd * this%eqnsclfac)
566  end do
567  end if
568  end subroutine uze_fc_expanded
569 
570  !> @brief Explicit solve
571  !!
572  !! There should be no explicit solve for uze. However, if there were, then
573  !! this subroutine would add terms specific to the unsaturated zone to the
574  !! explicit unsaturated-zone solve
575  subroutine uze_solve(this)
576  ! -- dummy
577  class(gweuzetype) :: this
578  ! -- local
579  integer(I4B) :: j
580  integer(I4B) :: n1, n2
581  real(DP) :: rrate
582  !
583  ! -- Add infiltration contribution
584  if (this%idxbudinfl /= 0) then
585  do j = 1, this%flowbudptr%budterm(this%idxbudinfl)%nlist
586  call this%uze_infl_term(j, n1, n2, rrate)
587  this%dbuff(n1) = this%dbuff(n1) + rrate
588  end do
589  end if
590  !
591  ! -- Add rejected infiltration contribution
592  if (this%idxbudrinf /= 0) then
593  do j = 1, this%flowbudptr%budterm(this%idxbudrinf)%nlist
594  call this%uze_rinf_term(j, n1, n2, rrate)
595  this%dbuff(n1) = this%dbuff(n1) + rrate
596  end do
597  end if
598  !
599  ! -- Add unsaturated et contribution
600  if (this%idxbuduzet /= 0) then
601  do j = 1, this%flowbudptr%budterm(this%idxbuduzet)%nlist
602  call this%uze_uzet_term(j, n1, n2, rrate)
603  this%dbuff(n1) = this%dbuff(n1) + rrate
604  end do
605  end if
606  !
607  ! -- Add rejected infiltration to mover contribution
608  if (this%idxbudritm /= 0) then
609  do j = 1, this%flowbudptr%budterm(this%idxbudritm)%nlist
610  call this%uze_ritm_term(j, n1, n2, rrate)
611  this%dbuff(n1) = this%dbuff(n1) + rrate
612  end do
613  end if
614  end subroutine uze_solve
615 
616  !> @brief Return the number of UZE-specific budget terms
617  !!
618  !! Function to return the number of budget terms just for this package.
619  !! This overrides function in parent.
620  !<
621  function uze_get_nbudterms(this) result(nbudterms)
622  ! -- dummy
623  class(gweuzetype) :: this
624  ! -- Return
625  integer(I4B) :: nbudterms
626  !
627  ! -- Number of budget terms is 5
628  nbudterms = 0
629  if (this%idxbudinfl /= 0) nbudterms = nbudterms + 1
630  if (this%idxbudrinf /= 0) nbudterms = nbudterms + 1
631  if (this%idxbuduzet /= 0) nbudterms = nbudterms + 1
632  if (this%idxbudritm /= 0) nbudterms = nbudterms + 1
633  if (this%idxbudtheq /= 0) nbudterms = nbudterms + 1
634  end function uze_get_nbudterms
635 
636  !> @brief Override similarly named function in APT
637  !!
638  !! Set the temperature to be used by MVE as the user-specified
639  !! temperature applied to the infiltration
640  !<
641  function get_mvr_depvar(this)
642  ! -- dummy
643  class(gweuzetype) :: this
644  ! -- return
645  real(dp), dimension(:), contiguous, pointer :: get_mvr_depvar
646  !
647  get_mvr_depvar => this%tempinfl
648  end function get_mvr_depvar
649 
650  !> @brief Setup budget object
651  !!
652  !! Set up the budget object that stores all the unsaturated-zone
653  !! flows
654  !<
655  subroutine uze_setup_budobj(this, idx)
656  ! -- modules
657  use constantsmodule, only: lenbudtxt
658  ! -- dummy
659  class(gweuzetype) :: this
660  integer(I4B), intent(inout) :: idx
661  ! -- local
662  integer(I4B) :: maxlist, naux
663  character(len=LENBUDTXT) :: text
664  !
665  ! -- Infiltration
666  text = ' INFILTRATION'
667  idx = idx + 1
668  maxlist = this%flowbudptr%budterm(this%idxbudinfl)%maxlist
669  naux = 0
670  call this%budobj%budterm(idx)%initialize(text, &
671  this%name_model, &
672  this%packName, &
673  this%name_model, &
674  this%packName, &
675  maxlist, .false., .false., &
676  naux)
677  !
678  ! -- Rejected infiltration (Hortonian flow)
679  if (this%idxbudrinf /= 0) then
680  text = ' REJ-INF'
681  idx = idx + 1
682  maxlist = this%flowbudptr%budterm(this%idxbudrinf)%maxlist
683  naux = 0
684  call this%budobj%budterm(idx)%initialize(text, &
685  this%name_model, &
686  this%packName, &
687  this%name_model, &
688  this%packName, &
689  maxlist, .false., .false., &
690  naux)
691  end if
692  !
693  ! -- Evapotranspiration from the unsaturated zone
694  if (this%idxbuduzet /= 0) then
695  text = ' UZET'
696  idx = idx + 1
697  maxlist = this%flowbudptr%budterm(this%idxbuduzet)%maxlist
698  naux = 0
699  call this%budobj%budterm(idx)%initialize(text, &
700  this%name_model, &
701  this%packName, &
702  this%name_model, &
703  this%packName, &
704  maxlist, .false., .false., &
705  naux)
706  end if
707  !
708  ! -- Rejected infiltration that is subsequently transferred by MVR
709  if (this%idxbudritm /= 0) then
710  text = ' INF-REJ-TO-MVR'
711  idx = idx + 1
712  maxlist = this%flowbudptr%budterm(this%idxbudritm)%maxlist
713  naux = 0
714  call this%budobj%budterm(idx)%initialize(text, &
715  this%name_model, &
716  this%packName, &
717  this%name_model, &
718  this%packName, &
719  maxlist, .false., .false., &
720  naux)
721  end if
722  !
723  ! -- Energy transferred to solid phase by the thermal equilibrium assumption
724  text = ' THERMAL-EQUIL'
725  idx = idx + 1
726  ! -- use dimension of GWF term
727  maxlist = this%flowbudptr%budterm(this%idxbudgwf)%maxlist
728  naux = 0
729  call this%budobj%budterm(idx)%initialize(text, &
730  this%name_model, &
731  this%packName, &
732  this%name_model, &
733  this%packName, &
734  maxlist, .false., .false., &
735  naux)
736  end subroutine uze_setup_budobj
737 
738  !> @brief Fill UZE budget object
739  !!
740  !! Copy flow terms into this%budobj
741  !<
742  subroutine uze_fill_budobj(this, idx, x, flowja, ccratin, ccratout)
743  ! -- modules
744  ! -- dummy
745  class(gweuzetype) :: this
746  integer(I4B), intent(inout) :: idx
747  real(DP), dimension(:), intent(in) :: x
748  real(DP), dimension(:), contiguous, intent(inout) :: flowja
749  real(DP), intent(inout) :: ccratin
750  real(DP), intent(inout) :: ccratout
751  ! -- local
752  integer(I4B) :: j, n1, n2, indx
753  integer(I4B) :: nlist, nlen
754  integer(I4B) :: igwfnode
755  integer(I4B) :: idiag
756  real(DP) :: q
757  real(DP), dimension(:), allocatable :: budresid
758  !
759  allocate (budresid(this%ncv))
760  do n1 = 1, this%ncv
761  budresid(n1) = dzero
762  end do
763  !
764  indx = 0
765  !
766  ! -- FLOW JA FACE into budresid
767  nlen = 0
768  if (this%idxbudfjf /= 0) then
769  nlen = this%flowbudptr%budterm(this%idxbudfjf)%maxlist
770  end if
771  if (nlen > 0) then
772  indx = indx + 1
773  nlist = this%budobj%budterm(indx)%nlist
774  do j = 1, nlist
775  n1 = this%budobj%budterm(indx)%id1(j)
776  n2 = this%budobj%budterm(indx)%id2(j)
777  if (n1 < n2) then
778  q = this%budobj%budterm(indx)%flow(j)
779  budresid(n1) = budresid(n1) + q
780  budresid(n2) = budresid(n2) - q
781  end if
782  end do
783  end if
784  !
785  ! -- GWF (LEAKAGE) into budresid
786  indx = indx + 1
787  nlist = this%budobj%budterm(indx)%nlist
788  do j = 1, nlist
789  n1 = this%budobj%budterm(indx)%id1(j)
790  q = this%budobj%budterm(indx)%flow(j)
791  budresid(n1) = budresid(n1) + q
792  end do
793  !
794  ! -- Skip individual package terms
795  indx = this%idxlastpak
796  !
797  ! -- STORAGE into budresid
798  indx = indx + 1
799  do n1 = 1, this%ncv
800  q = this%budobj%budterm(indx)%flow(n1)
801  budresid(n1) = budresid(n1) + q
802  end do
803  !
804  ! -- TO MOVER into budresid
805  if (this%idxbudtmvr /= 0) then
806  indx = indx + 1
807  nlist = this%budobj%budterm(indx)%nlist
808  do j = 1, nlist
809  n1 = this%budobj%budterm(indx)%id1(j)
810  q = this%budobj%budterm(indx)%flow(j)
811  budresid(n1) = budresid(n1) + q
812  end do
813  end if
814  !
815  ! -- FROM MOVER into budresid
816  if (this%idxbudfmvr /= 0) then
817  indx = indx + 1
818  nlist = this%budobj%budterm(indx)%nlist
819  do j = 1, nlist
820  n1 = this%budobj%budterm(indx)%id1(j)
821  q = this%budobj%budterm(indx)%flow(j)
822  budresid(n1) = budresid(n1) + q
823  end do
824  end if
825  !
826  ! -- CONSTANT FLOW into budresid
827  indx = indx + 1
828  do n1 = 1, this%ncv
829  q = this%budobj%budterm(indx)%flow(n1)
830  budresid(n1) = budresid(n1) + q
831  end do
832  !
833  ! -- AUXILIARY VARIABLES into budresid
834  ! -- (No flows associated with these)
835  !
836  ! -- Individual package terms processed last
837  !
838  ! -- Infiltration
839  idx = idx + 1
840  nlist = this%flowbudptr%budterm(this%idxbudinfl)%nlist
841  call this%budobj%budterm(idx)%reset(nlist)
842  do j = 1, nlist
843  call this%uze_infl_term(j, n1, n2, q)
844  call this%budobj%budterm(idx)%update_term(n1, n2, q)
845  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
846  budresid(n1) = budresid(n1) + q
847  end do
848  !
849  ! -- Rej-Inf
850  if (this%idxbudrinf /= 0) then
851  idx = idx + 1
852  nlist = this%flowbudptr%budterm(this%idxbudrinf)%nlist
853  call this%budobj%budterm(idx)%reset(nlist)
854  do j = 1, nlist
855  call this%uze_rinf_term(j, n1, n2, q)
856  call this%budobj%budterm(idx)%update_term(n1, n2, q)
857  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
858  budresid(n1) = budresid(n1) + q
859  end do
860  end if
861  !
862  ! -- UZET
863  if (this%idxbuduzet /= 0) then
864  idx = idx + 1
865  nlist = this%flowbudptr%budterm(this%idxbuduzet)%nlist
866  call this%budobj%budterm(idx)%reset(nlist)
867  do j = 1, nlist
868  call this%uze_uzet_term(j, n1, n2, q)
869  call this%budobj%budterm(idx)%update_term(n1, n2, q)
870  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
871  budresid(n1) = budresid(n1) + q
872  end do
873  end if
874  !
875  ! -- Rej-Inf-To-MVR
876  if (this%idxbudritm /= 0) then
877  idx = idx + 1
878  nlist = this%flowbudptr%budterm(this%idxbudritm)%nlist
879  call this%budobj%budterm(idx)%reset(nlist)
880  do j = 1, nlist
881  call this%uze_ritm_term(j, n1, n2, q)
882  call this%budobj%budterm(idx)%update_term(n1, n2, q)
883  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
884  budresid(n1) = budresid(n1) + q
885  end do
886  end if
887  !
888  ! -- Thermal-Equil
889  ! -- processed last because it is calculated from the residual
890  idx = idx + 1
891  nlist = this%flowbudptr%budterm(this%idxbudgwf)%nlist
892  call this%budobj%budterm(idx)%reset(nlist)
893  do j = 1, nlist
894  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
895  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(j)
896  q = -budresid(n1)
897  call this%uze_theq_term(j, n1, igwfnode, q)
898  call this%budobj%budterm(idx)%update_term(n1, igwfnode, q)
899  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
900  if (this%iboundpak(n1) /= 0) then
901  ! -- Contribution to gwe cell budget
902  this%simvals(n1) = this%simvals(n1) - q
903  idiag = this%dis%con%ia(igwfnode)
904  flowja(idiag) = flowja(idiag) - q
905  end if
906  end do
907  !
908  deallocate (budresid)
909  end subroutine uze_fill_budobj
910 
911  !> @brief Allocate scalars
912  !!
913  !! Allocate scalars specific to UZE package
914  !<
915  subroutine allocate_scalars(this)
916  ! -- modules
918  ! -- dummy
919  class(gweuzetype) :: this
920  ! -- local
921  !
922  ! -- Allocate scalars in TspAptType
923  call this%TspAptType%allocate_scalars()
924  !
925  ! -- Allocate
926  call mem_allocate(this%idxbudinfl, 'IDXBUDINFL', this%memoryPath)
927  call mem_allocate(this%idxbudrinf, 'IDXBUDRINF', this%memoryPath)
928  call mem_allocate(this%idxbuduzet, 'IDXBUDUZET', this%memoryPath)
929  call mem_allocate(this%idxbudritm, 'IDXBUDRITM', this%memoryPath)
930  call mem_allocate(this%idxbudtheq, 'IDXBUDTHEQ', this%memoryPath)
931  !
932  ! -- Initialize
933  this%idxbudinfl = 0
934  this%idxbudrinf = 0
935  this%idxbuduzet = 0
936  this%idxbudritm = 0
937  this%idxbudtheq = 0
938  end subroutine allocate_scalars
939 
940  !> @brief Allocate arrays
941  !!
942  !! Allocate arrays used by the UZE package
943  !<
944  subroutine uze_allocate_arrays(this)
945  ! -- modules
947  ! -- dummy
948  class(gweuzetype), intent(inout) :: this
949  ! -- local
950  integer(I4B) :: n
951  !
952  ! -- Time series
953  call mem_allocate(this%tempinfl, this%ncv, 'TEMPINFL', this%memoryPath)
954  call mem_allocate(this%tempuzet, this%ncv, 'TEMPUZET', this%memoryPath)
955  !
956  ! -- Call standard TspAptType allocate arrays
957  call this%TspAptType%apt_allocate_arrays()
958  !
959  ! -- Initialize
960  do n = 1, this%ncv
961  this%tempinfl(n) = dzero
962  this%tempuzet(n) = dzero
963  end do
964  end subroutine uze_allocate_arrays
965 
966  !> @brief Deallocate memory
967  !<
968  subroutine uze_da(this)
969  ! -- modules
971  ! -- dummy
972  class(gweuzetype) :: this
973  !
974  ! -- Deallocate scalars
975  call mem_deallocate(this%idxbudinfl)
976  call mem_deallocate(this%idxbudrinf)
977  call mem_deallocate(this%idxbuduzet)
978  call mem_deallocate(this%idxbudritm)
979  call mem_deallocate(this%idxbudtheq)
980  !
981  ! -- Deallocate time series
982  call mem_deallocate(this%tempinfl)
983  call mem_deallocate(this%tempuzet)
984  !
985  ! -- Deallocate scalars in TspAptType
986  call this%TspAptType%bnd_da()
987  end subroutine uze_da
988 
989  !> @brief Infiltration term
990  !!
991  !! Accounts for energy added to the subsurface via infiltration, for example,
992  !! energy entering the model domain via rainfall or irrigation.
993  !<
994  subroutine uze_infl_term(this, ientry, n1, n2, rrate, &
995  rhsval, hcofval)
996  ! -- dummy
997  class(gweuzetype) :: this
998  integer(I4B), intent(in) :: ientry
999  integer(I4B), intent(inout) :: n1
1000  integer(I4B), intent(inout) :: n2
1001  real(DP), intent(inout), optional :: rrate
1002  real(DP), intent(inout), optional :: rhsval
1003  real(DP), intent(inout), optional :: hcofval
1004  ! -- local
1005  real(DP) :: qbnd
1006  real(DP) :: ctmp
1007  real(DP) :: h, r
1008  !
1009  n1 = this%flowbudptr%budterm(this%idxbudinfl)%id1(ientry)
1010  n2 = this%flowbudptr%budterm(this%idxbudinfl)%id2(ientry)
1011  !
1012  ! -- Note that qbnd is negative for negative infiltration
1013  qbnd = this%flowbudptr%budterm(this%idxbudinfl)%flow(ientry)
1014  if (qbnd < dzero) then
1015  ctmp = this%xnewpak(n1)
1016  h = qbnd
1017  r = dzero
1018  else
1019  ctmp = this%tempinfl(n1)
1020  h = dzero
1021  r = -qbnd * ctmp
1022  end if
1023  if (present(rrate)) rrate = qbnd * ctmp * this%eqnsclfac
1024  if (present(rhsval)) rhsval = r * this%eqnsclfac
1025  if (present(hcofval)) hcofval = h * this%eqnsclfac
1026  end subroutine uze_infl_term
1027 
1028  !> @brief Rejected infiltration term
1029  !!
1030  !! Accounts for energy that is added to the model from specifying an
1031  !! infiltration rate and temperature, but is subsequently removed from
1032  !! the model as that portion of the infiltration that is rejected (and
1033  !! NOT transferred to another advanced package via the MVR/MVT packages).
1034  !<
1035  subroutine uze_rinf_term(this, ientry, n1, n2, rrate, &
1036  rhsval, hcofval)
1037  ! -- dummy
1038  class(gweuzetype) :: this
1039  integer(I4B), intent(in) :: ientry
1040  integer(I4B), intent(inout) :: n1
1041  integer(I4B), intent(inout) :: n2
1042  real(DP), intent(inout), optional :: rrate
1043  real(DP), intent(inout), optional :: rhsval
1044  real(DP), intent(inout), optional :: hcofval
1045  ! -- local
1046  real(DP) :: qbnd
1047  real(DP) :: ctmp
1048  !
1049  n1 = this%flowbudptr%budterm(this%idxbudrinf)%id1(ientry)
1050  n2 = this%flowbudptr%budterm(this%idxbudrinf)%id2(ientry)
1051  qbnd = this%flowbudptr%budterm(this%idxbudrinf)%flow(ientry)
1052  ctmp = this%tempinfl(n1)
1053  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
1054  if (present(rhsval)) rhsval = dzero
1055  if (present(hcofval)) hcofval = qbnd * this%eqnsclfac
1056  end subroutine uze_rinf_term
1057 
1058  !> @brief Evapotranspiration from the unsaturated-zone term
1059  !!
1060  !! Accounts for thermal cooling in the unsaturated zone as a result of
1061  !! evapotranspiration from the unsaturated zone. Amount of water converted
1062  !! to vapor phase (UZET) determined by GWF model
1063  !<
1064  subroutine uze_uzet_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
1065  ! -- dummy
1066  class(gweuzetype) :: this
1067  integer(I4B), intent(in) :: ientry
1068  integer(I4B), intent(inout) :: n1
1069  integer(I4B), intent(inout) :: n2
1070  real(DP), intent(inout), optional :: rrate
1071  real(DP), intent(inout), optional :: rhsval
1072  real(DP), intent(inout), optional :: hcofval
1073  ! -- local
1074  real(DP) :: qbnd
1075  real(DP) :: ctmp
1076  real(DP) :: omega
1077  !
1078  n1 = this%flowbudptr%budterm(this%idxbuduzet)%id1(ientry)
1079  n2 = this%flowbudptr%budterm(this%idxbuduzet)%id2(ientry)
1080  ! -- Note that qbnd is negative for uzet
1081  qbnd = this%flowbudptr%budterm(this%idxbuduzet)%flow(ientry)
1082  ctmp = this%tempuzet(n1)
1083  if (this%xnewpak(n1) < ctmp) then
1084  omega = done
1085  else
1086  omega = dzero
1087  end if
1088  if (present(rrate)) &
1089  rrate = (omega * qbnd * this%xnewpak(n1) + &
1090  (done - omega) * qbnd * ctmp) * this%eqnsclfac
1091  if (present(rhsval)) rhsval = -(done - omega) * qbnd * ctmp * this%eqnsclfac
1092  if (present(hcofval)) hcofval = omega * qbnd * this%eqnsclfac
1093  end subroutine uze_uzet_term
1094 
1095  !> @brief Rejected infiltration to MVR/MVT term
1096  !!
1097  !! Accounts for energy that is added to the model from specifying an
1098  !! infiltration rate and temperature, but does not infiltrate into the
1099  !! subsurface. This subroutine is called when the rejected infiltration
1100  !! is transferred to another advanced package via the MVR/MVT packages.
1101  !<
1102  subroutine uze_ritm_term(this, ientry, n1, n2, rrate, &
1103  rhsval, hcofval)
1104  ! -- dummy
1105  class(gweuzetype) :: this
1106  integer(I4B), intent(in) :: ientry
1107  integer(I4B), intent(inout) :: n1
1108  integer(I4B), intent(inout) :: n2
1109  real(DP), intent(inout), optional :: rrate
1110  real(DP), intent(inout), optional :: rhsval
1111  real(DP), intent(inout), optional :: hcofval
1112  ! -- local
1113  real(DP) :: qbnd
1114  real(DP) :: ctmp
1115  !
1116  n1 = this%flowbudptr%budterm(this%idxbudritm)%id1(ientry)
1117  n2 = this%flowbudptr%budterm(this%idxbudritm)%id2(ientry)
1118  qbnd = this%flowbudptr%budterm(this%idxbudritm)%flow(ientry)
1119  ctmp = this%tempinfl(n1)
1120  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
1121  if (present(rhsval)) rhsval = dzero
1122  if (present(hcofval)) hcofval = qbnd * this%eqnsclfac
1123  end subroutine uze_ritm_term
1124 
1125  !> @brief Heat transferred through thermal equilibrium with the solid phase
1126  !!
1127  !! Accounts for the transfer of energy from the liquid phase to the solid
1128  !! phase as a result of the instantaneous thermal equilibrium assumption.
1129  !<
1130  subroutine uze_theq_term(this, ientry, n1, n2, rrate)
1131  ! -- modules
1132  use constantsmodule, only: lenbudtxt
1133  ! -- dummy
1134  class(gweuzetype) :: this
1135  integer(I4B), intent(in) :: ientry
1136  integer(I4B), intent(inout) :: n1
1137  integer(I4B), intent(inout) :: n2
1138  real(DP), intent(inout) :: rrate
1139  ! -- local
1140  real(DP) :: r
1141  integer(I4B) :: i
1142  character(len=LENBUDTXT) :: flowtype
1143  !
1144  r = dzero
1145  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(ientry)
1146  n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(ientry)
1147  if (this%iboundpak(n1) /= 0) then
1148  do i = 1, this%budobj%nbudterm
1149  flowtype = this%budobj%budterm(i)%flowtype
1150  select case (trim(adjustl(flowtype)))
1151  case ('THERMAL-EQUIL')
1152  ! -- Skip
1153  continue
1154  case default
1155  r = r - this%budobj%budterm(i)%flow(ientry)
1156  end select
1157  end do
1158  end if
1159  rrate = r
1160  end subroutine uze_theq_term
1161 
1162  !> @brief Define UZE Observation
1163  !!
1164  !! This subroutine:
1165  !! - Stores observation types supported by the parent APT package.
1166  !! - Overrides BndType%bnd_df_obs
1167  !<
1168  subroutine uze_df_obs(this)
1169  ! -- dummy
1170  class(gweuzetype) :: this
1171  ! -- local
1172  integer(I4B) :: indx
1173  !
1174  ! -- Store obs type and assign procedure pointer
1175  ! for temperature observation type.
1176  call this%obs%StoreObsType('temperature', .false., indx)
1177  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1178  !
1179  ! -- Store obs type and assign procedure pointer
1180  ! for flow between uze cells.
1181  call this%obs%StoreObsType('flow-ja-face', .true., indx)
1182  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid12
1183  !
1184  ! -- Store obs type and assign procedure pointer
1185  ! for from-mvr observation type.
1186  call this%obs%StoreObsType('from-mvr', .true., indx)
1187  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1188  !
1189  ! -- to-mvr not supported for uze
1190  !call this%obs%StoreObsType('to-mvr', .true., indx)
1191  !this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID
1192  !
1193  ! -- Store obs type and assign procedure pointer
1194  ! for storage observation type.
1195  call this%obs%StoreObsType('storage', .true., indx)
1196  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1197  !
1198  ! -- Store obs type and assign procedure pointer
1199  ! for constant observation type.
1200  call this%obs%StoreObsType('constant', .true., indx)
1201  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1202  !
1203  ! -- Store obs type and assign procedure pointer
1204  ! for observation type: uze
1205  call this%obs%StoreObsType('uze', .true., indx)
1206  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1207  !
1208  ! -- Store obs type and assign procedure pointer
1209  ! for observation type.
1210  call this%obs%StoreObsType('infiltration', .true., indx)
1211  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1212  !
1213  ! -- Store obs type and assign procedure pointer
1214  ! for observation type.
1215  call this%obs%StoreObsType('rej-inf', .true., indx)
1216  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1217  !
1218  ! -- Store obs type and assign procedure pointer
1219  ! for observation type.
1220  call this%obs%StoreObsType('uzet', .true., indx)
1221  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1222  !
1223  ! -- Store obs type and assign procedure pointer
1224  ! for observation type.
1225  call this%obs%StoreObsType('rej-inf-to-mvr', .true., indx)
1226  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1227  !
1228  ! -- Store obs type and assign procedure pointer
1229  ! for observation type.
1230  call this%obs%StoreObsType('thermal-equil', .true., indx)
1231  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
1232  end subroutine uze_df_obs
1233 
1234  !> @brief Process package specific obs
1235  !!
1236  !! Method to process specific observations for this package.
1237  !<
1238  subroutine uze_rp_obs(this, obsrv, found)
1239  ! -- dummy
1240  class(gweuzetype), intent(inout) :: this !< package class
1241  type(observetype), intent(inout) :: obsrv !< observation object
1242  logical, intent(inout) :: found !< indicate whether observation was found
1243  !
1244  found = .true.
1245  select case (obsrv%ObsTypeId)
1246  case ('INFILTRATION')
1247  call this%rp_obs_byfeature(obsrv)
1248  case ('REJ-INF')
1249  call this%rp_obs_byfeature(obsrv)
1250  case ('UZET')
1251  call this%rp_obs_byfeature(obsrv)
1252  case ('REJ-INF-TO-MVR')
1253  call this%rp_obs_byfeature(obsrv)
1254  case ('THERMAL-EQUIL')
1255  call this%rp_obs_byfeature(obsrv)
1256  case default
1257  found = .false.
1258  end select
1259  end subroutine uze_rp_obs
1260 
1261  !> @brief Calculate observation value and pass it back to APT
1262  !<
1263  subroutine uze_bd_obs(this, obstypeid, jj, v, found)
1264  ! -- dummy
1265  class(gweuzetype), intent(inout) :: this
1266  character(len=*), intent(in) :: obstypeid
1267  real(DP), intent(inout) :: v
1268  integer(I4B), intent(in) :: jj
1269  logical, intent(inout) :: found
1270  ! -- local
1271  integer(I4B) :: n1, n2
1272  !
1273  found = .true.
1274  select case (obstypeid)
1275  case ('INFILTRATION')
1276  if (this%iboundpak(jj) /= 0 .and. this%idxbudinfl > 0) then
1277  call this%uze_infl_term(jj, n1, n2, v)
1278  end if
1279  case ('REJ-INF')
1280  if (this%iboundpak(jj) /= 0 .and. this%idxbudrinf > 0) then
1281  call this%uze_rinf_term(jj, n1, n2, v)
1282  end if
1283  case ('UZET')
1284  if (this%iboundpak(jj) /= 0 .and. this%idxbuduzet > 0) then
1285  call this%uze_uzet_term(jj, n1, n2, v)
1286  end if
1287  case ('REJ-INF-TO-MVR')
1288  if (this%iboundpak(jj) /= 0 .and. this%idxbudritm > 0) then
1289  call this%uze_ritm_term(jj, n1, n2, v)
1290  end if
1291  case ('THERMAL-EQUIL')
1292  if (this%iboundpak(jj) /= 0 .and. this%idxbudtheq > 0) then
1293  call this%uze_theq_term(jj, n1, n2, v)
1294  end if
1295  case default
1296  found = .false.
1297  end select
1298  end subroutine uze_bd_obs
1299 
1300  !> @brief Check if UZF object area is not equal to the cell area
1301  !!
1302  !! GWE equilibrates the temperature of a UZE object with the host cell. UZE
1303  !! assumes that the area associated with a specific UZE object is equal to
1304  !! the area of the host cell. When this condition is not true, the code
1305  !! exits with an appropriate message.
1306  !<
1307  subroutine uze_ad_chk(this)
1308  ! -- modules
1309  use constantsmodule, only: izero
1310  use mathutilmodule, only: is_close
1311  use simmodule, only: count_errors
1312  ! -- dummy
1313  class(gweuzetype), intent(inout) :: this
1314  ! -- local
1315  integer(I4B) :: nuz
1316  integer(I4B) :: n
1317  integer(I4B) :: igwfnode
1318  real(DP) :: carea
1319  real(DP) :: uzarea
1320 
1321  nuz = this%flowbudptr%budterm(this%idxbudgwf)%maxlist
1322 
1323  ! cycle through uze objects, stop at first occurrence of more than one
1324  ! uze object in a cell
1325  do n = 1, nuz
1326  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(n)
1327  carea = this%dis%area(igwfnode)
1328  uzarea = this%flowbudptr%budterm(this%idxbudgwf)%auxvar(1, n)
1329  ! compare areas
1330  if (.not. is_close(carea, uzarea)) then
1331  call this%area_error(igwfnode)
1332  end if
1333  end do
1334  if (count_errors() > 0) then
1335  call this%parser%StoreErrorUnit()
1336  end if
1337  end subroutine uze_ad_chk
1338 
1339  !> @brief Print and store error msg indicating area of UZF object is not
1340  !! equal to that of the host cell
1341  !<
1342  subroutine area_error(this, iloc)
1343  ! -- modules
1344  use simvariablesmodule, only: errmsg
1345  ! -- dummy
1346  class(gweuzetype) :: this
1347  integer(I4B) :: iloc
1348  ! -- local
1349  character(len=30) :: nodestr
1350  !
1351  call this%dis%noder_to_string(iloc, nodestr)
1352  write (errmsg, '(3a)') &
1353  'In a GWE model, the area of every UZF object must be equal to that of &
1354  &the host cell. This condition is violated in cell ', &
1355  trim(adjustl(nodestr)), '. Check use of AUXMULTNAME option in UZF &
1356  &package.'
1357  call store_error(errmsg)
1358  end subroutine area_error
1359 
1360  !> @brief Sets the stress period attributes for keyword use.
1361  !<
1362  subroutine uze_set_stressperiod(this, itemno, keyword, found)
1363  ! -- modules
1365  ! -- dummy
1366  class(gweuzetype), intent(inout) :: this
1367  integer(I4B), intent(in) :: itemno
1368  character(len=*), intent(in) :: keyword
1369  logical, intent(inout) :: found
1370  ! -- local
1371  character(len=LINELENGTH) :: temp_text
1372  integer(I4B) :: ierr
1373  integer(I4B) :: jj
1374  real(DP), pointer :: bndElem => null()
1375  !
1376  ! INFILTRATION <infiltration>
1377  ! UZET <uzet>
1378  !
1379  found = .true.
1380  select case (keyword)
1381  case ('INFILTRATION')
1382  ierr = this%apt_check_valid(itemno)
1383  if (ierr /= 0) then
1384  goto 999
1385  end if
1386  call this%parser%GetString(temp_text)
1387  jj = 1
1388  bndelem => this%tempinfl(itemno)
1389  call read_value_or_time_series_adv(temp_text, itemno, jj, bndelem, &
1390  this%packName, 'BND', this%tsManager, &
1391  this%iprpak, 'INFILTRATION')
1392  case ('UZET')
1393  ierr = this%apt_check_valid(itemno)
1394  if (ierr /= 0) then
1395  goto 999
1396  end if
1397  call this%parser%GetString(temp_text)
1398  jj = 1
1399  bndelem => this%tempuzet(itemno)
1400  call read_value_or_time_series_adv(temp_text, itemno, jj, bndelem, &
1401  this%packName, 'BND', this%tsManager, &
1402  this%iprpak, 'UZET')
1403  case default
1404  !
1405  ! -- Keyword not recognized so return to caller with found = .false.
1406  found = .false.
1407  end select
1408  !
1409 999 continue
1410  end subroutine uze_set_stressperiod
1411 
1412 end module gweuzemodule
This module contains the base boundary package.
class(bndtype) function, pointer, public getbndfromlist(list, idx)
Get boundary from package list.
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
integer(i4b), parameter izero
integer constant zero
Definition: Constants.f90:51
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter lenbudtxt
maximum length of a budget component names
Definition: Constants.f90:37
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
subroutine uze_ad_chk(this)
Check if UZF object area is not equal to the cell area.
Definition: gwe-uze.f90:1308
subroutine uze_rp_obs(this, obsrv, found)
Process package specific obs.
Definition: gwe-uze.f90:1239
character(len= *), parameter ftype
Definition: gwe-uze.f90:45
subroutine uze_mc(this, moffset, matrix_sln)
Map package connection to matrix.
Definition: gwe-uze.f90:346
subroutine uze_ac(this, moffset, sparse)
Add package connection to matrix.
Definition: gwe-uze.f90:279
subroutine allocate_scalars(this)
Allocate scalars.
Definition: gwe-uze.f90:916
character(len=16) text
Definition: gwe-uze.f90:47
subroutine uze_setup_budobj(this, idx)
Setup budget object.
Definition: gwe-uze.f90:656
subroutine uze_ritm_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Rejected infiltration to MVR/MVT term.
Definition: gwe-uze.f90:1104
subroutine uze_df_obs(this)
Define UZE Observation.
Definition: gwe-uze.f90:1169
subroutine uze_set_stressperiod(this, itemno, keyword, found)
Sets the stress period attributes for keyword use.
Definition: gwe-uze.f90:1363
subroutine uze_bd_obs(this, obstypeid, jj, v, found)
Calculate observation value and pass it back to APT.
Definition: gwe-uze.f90:1264
subroutine uze_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
Add matrix terms related to UZE.
Definition: gwe-uze.f90:436
subroutine uze_theq_term(this, ientry, n1, n2, rrate)
Heat transferred through thermal equilibrium with the solid phase.
Definition: gwe-uze.f90:1131
subroutine, public uze_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, fmi, eqnsclfac, gwecommon, dvt, dvu, dvua)
Create a new UZE package.
Definition: gwe-uze.f90:96
subroutine uze_rinf_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Rejected infiltration term.
Definition: gwe-uze.f90:1037
integer(i4b) function uze_get_nbudterms(this)
Return the number of UZE-specific budget terms.
Definition: gwe-uze.f90:622
subroutine uze_solve(this)
Explicit solve.
Definition: gwe-uze.f90:576
character(len= *), parameter flowtype
Definition: gwe-uze.f90:46
subroutine uze_da(this)
Deallocate memory.
Definition: gwe-uze.f90:969
real(dp) function, dimension(:), pointer, contiguous get_mvr_depvar(this)
Override similarly named function in APT.
Definition: gwe-uze.f90:642
subroutine find_uze_package(this)
Find corresponding uze package.
Definition: gwe-uze.f90:156
subroutine uze_fill_budobj(this, idx, x, flowja, ccratin, ccratout)
Fill UZE budget object.
Definition: gwe-uze.f90:743
subroutine uze_allocate_arrays(this)
Allocate arrays.
Definition: gwe-uze.f90:945
subroutine area_error(this, iloc)
Print and store error msg indicating area of UZF object is not equal to that of the host cell.
Definition: gwe-uze.f90:1343
subroutine uze_infl_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Infiltration term.
Definition: gwe-uze.f90:996
subroutine uze_uzet_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Evapotranspiration from the unsaturated-zone term.
Definition: gwe-uze.f90:1065
This module defines variable data types.
Definition: kind.f90:8
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
Definition: MathUtil.f90:46
This module contains the derived types ObserveType and ObsDataType.
Definition: Observe.f90:15
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
subroutine, public read_value_or_time_series_adv(textInput, ii, jj, bndElem, pkgName, auxOrBnd, tsManager, iprpak, varName)
Call this subroutine from advanced packages to define timeseries link for a variable (varName).
subroutine, public apt_process_obsid(obsrv, dis, inunitobs, iout)
Process observation IDs for an advanced package.
Definition: tsp-apt.f90:2851
subroutine, public apt_process_obsid12(obsrv, dis, inunitobs, iout)
Process observation IDs for a package.
Definition: tsp-apt.f90:2894
@ brief BndType
Data for sharing among multiple packages. Originally read in from.