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