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