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