MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
tsp-apt.f90
Go to the documentation of this file.
1 ! -- Advanced Package Transport Module
2 ! -- This module contains most of the routines for simulating transport
3 ! -- through the advanced packages.
4 ! -- Future work:
5 ! * support decay, sorption
6 ! * dispersion in SFT and UZT?
7 !
8 ! AFP flows (flowbudptr) index var ATP term Transport Type
9 !---------------------------------------------------------------------------------
10 
11 ! -- specialized terms in the flow budget
12 ! FLOW-JA-FACE idxbudfjf FLOW-JA-FACE cv2cv
13 ! GWF (aux FLOW-AREA) idxbudgwf GWF cv2gwf
14 ! STORAGE (aux VOLUME) idxbudsto none used for cv volumes
15 ! FROM-MVR idxbudfmvr FROM-MVR q * cext = this%qfrommvr(:) ! rhow*cpw is applied to various terms for heat transport
16 ! TO-MVR idxbudtmvr TO-MVR q * cfeat
17 
18 ! -- generalized source/sink terms (except ET?)
19 ! RAINFALL idxbudrain RAINFALL q * crain
20 ! EVAPORATION idxbudevap EVAPORATION cfeat<cevap: q*cfeat, else: q*cevap ! latent heat may be applied for evaporative cooling
21 ! RUNOFF idxbudroff RUNOFF q * croff
22 ! EXT-INFLOW idxbudiflw EXT-INFLOW q * ciflw
23 ! WITHDRAWAL idxbudwdrl WITHDRAWAL q * cfeat
24 ! EXT-OUTFLOW idxbudoutf EXT-OUTFLOW q * cfeat
25 
26 ! -- terms from a flow file that should be skipped
27 ! CONSTANT none none none
28 ! AUXILIARY none none none
29 
30 ! -- terms that are written to the transport budget file
31 ! none none STORAGE (aux MASS) dM/dt
32 ! none none AUXILIARY none
33 ! none none CONSTANT accumulate
34 !
35 !
37 
38  use kindmodule, only: dp, i4b, lgp
45  use simvariablesmodule, only: errmsg
46  use bndmodule, only: bndtype
47  use tspfmimodule, only: tspfmitype
50  use tablemodule, only: tabletype, table_cr
51  use observemodule, only: observetype
53  use basedismodule, only: disbasetype
55 
56  implicit none
57 
58  public :: tspapttype
59  public :: apt_process_obsid
60  public :: apt_process_obsid12
61 
62  character(len=LENFTYPE) :: ftype = 'APT'
63  character(len=LENVARNAME) :: text = ' APT'
64 
65  type, extends(bndtype) :: tspapttype
66 
67  character(len=LENPACKAGENAME) :: flowpackagename = '' !< name of corresponding flow package
68  character(len=8), &
69  dimension(:), pointer, contiguous :: status => null() !< active, inactive, constant
70  character(len=LENAUXNAME) :: cauxfpconc = '' !< name of aux column in flow package auxvar array for concentration (or temperature)
71  integer(I4B), pointer :: iauxfpconc => null() !< column in flow package bound array to insert concs
72  integer(I4B), pointer :: imatrows => null() !< if active, add new rows to matrix
73  integer(I4B), pointer :: iprconc => null() !< print conc to listing file
74  integer(I4B), pointer :: iconcout => null() !< unit number for conc output file
75  integer(I4B), pointer :: ibudgetout => null() !< unit number for budget output file
76  integer(I4B), pointer :: ibudcsv => null() !< unit number for csv budget output file
77  integer(I4B), pointer :: ncv => null() !< number of control volumes
78  integer(I4B), pointer :: igwfaptpak => null() !< package number of corresponding this package
79  integer(I4B), pointer :: idxprepak => null() !< budget-object index that precedes package-specific budget objects
80  integer(I4B), pointer :: idxlastpak => null() !< budget-object index of last package-specific budget object
81  real(dp), dimension(:), pointer, contiguous :: strt => null() !< starting feature concentration (or temperature)
82  real(dp), dimension(:), pointer, contiguous :: ktf => null() !< thermal conductivity between the apt and groundwater cell
83  real(dp), dimension(:), pointer, contiguous :: rfeatthk => null() !< thickness of streambed/lakebed/filter-pack material through which thermal conduction occurs
84  integer(I4B), dimension(:), pointer, contiguous :: idxlocnode => null() !< map position in global rhs and x array of pack entry
85  integer(I4B), dimension(:), pointer, contiguous :: idxpakdiag => null() !< map diag position of feature in global amat
86  integer(I4B), dimension(:), pointer, contiguous :: idxdglo => null() !< map position in global array of package diagonal row entries
87  integer(I4B), dimension(:), pointer, contiguous :: idxoffdglo => null() !< map position in global array of package off diagonal row entries
88  integer(I4B), dimension(:), pointer, contiguous :: idxsymdglo => null() !< map position in global array of package diagonal entries to model rows
89  integer(I4B), dimension(:), pointer, contiguous :: idxsymoffdglo => null() !< map position in global array of package off diagonal entries to model rows
90  integer(I4B), dimension(:), pointer, contiguous :: idxfjfdglo => null() !< map diagonal feature to feature in global amat
91  integer(I4B), dimension(:), pointer, contiguous :: idxfjfoffdglo => null() !< map off diagonal feature to feature in global amat
92  integer(I4B), dimension(:), pointer, contiguous :: iboundpak => null() !< package ibound
93  real(dp), dimension(:), pointer, contiguous :: xnewpak => null() !< feature concentration (or temperature) for current time step
94  real(dp), dimension(:), pointer, contiguous :: xoldpak => null() !< feature concentration (or temperature) from previous time step
95  real(dp), dimension(:), pointer, contiguous :: dbuff => null() !< temporary storage array
96  character(len=LENBOUNDNAME), &
97  dimension(:), pointer, contiguous :: featname => null()
98  real(dp), dimension(:), pointer, contiguous :: concfeat => null() !< concentration (or temperature) of the feature
99  real(dp), dimension(:, :), pointer, contiguous :: lauxvar => null() !< auxiliary variable
100  type(tspfmitype), pointer :: fmi => null() !< pointer to fmi object
101  real(dp), dimension(:), pointer, contiguous :: qsto => null() !< mass (or energy) flux due to storage change
102  real(dp), dimension(:), pointer, contiguous :: ccterm => null() !< mass (or energy) flux required to maintain constant concentration (or temperature)
103  integer(I4B), pointer :: idxbudfjf => null() !< index of flow ja face in flowbudptr
104  integer(I4B), pointer :: idxbudgwf => null() !< index of gwf terms in flowbudptr
105  integer(I4B), pointer :: idxbudsto => null() !< index of storage terms in flowbudptr
106  integer(I4B), pointer :: idxbudtmvr => null() !< index of to mover terms in flowbudptr
107  integer(I4B), pointer :: idxbudfmvr => null() !< index of from mover terms in flowbudptr
108  integer(I4B), pointer :: idxbudaux => null() !< index of auxiliary terms in flowbudptr
109  integer(I4B), dimension(:), pointer, contiguous :: idxbudssm => null() !< flag that flowbudptr%buditem is a general solute source/sink
110  integer(I4B), pointer :: nconcbudssm => null() !< number of concbudssm terms (columns)
111  real(dp), dimension(:, :), pointer, contiguous :: concbudssm => null() !< user specified concentrations (or temperatures) for flow terms
112  real(dp), dimension(:), pointer, contiguous :: qmfrommvr => null() !< a mass or energy flow coming from the mover that needs to be added
113  real(dp), pointer :: eqnsclfac => null() !< governing equation scale factor; =1. for solute; =rhow*cpw for energy
114  character(len=LENVARNAME) :: depvartype = '' !< stores string identifying dependent variable type, depending on model type
115  character(len=LENVARNAME) :: depvarunit = '' !< "mass" or "energy"
116  character(len=LENVARNAME) :: depvarunitabbrev = '' !< "M" or "E"
117  !
118  ! -- pointer to flow package boundary
119  type(bndtype), pointer :: flowpackagebnd => null()
120  !
121  ! -- budget objects
122  type(budgetobjecttype), pointer :: budobj => null() !< apt solute budget object
123  type(budgetobjecttype), pointer :: flowbudptr => null() !< GWF flow budget object
124  !
125  ! -- table objects
126  type(tabletype), pointer :: dvtab => null()
127 
128  contains
129 
130  procedure :: set_pointers => apt_set_pointers
131  procedure :: bnd_ac => apt_ac
132  procedure :: bnd_mc => apt_mc
133  procedure :: bnd_ar => apt_ar
134  procedure :: bnd_rp => apt_rp
135  procedure :: bnd_ad => apt_ad
136  procedure :: bnd_reset => apt_reset
137  procedure :: bnd_fc => apt_fc
138  procedure, public :: apt_fc_expanded ! Made public for uze
139  procedure, public :: apt_ad_chk
140  procedure :: pak_fc_expanded
141  procedure, private :: apt_fc_nonexpanded
142  procedure, public :: apt_cfupdate ! Made public for uze
143  procedure :: apt_check_valid
144  procedure :: apt_set_stressperiod
145  procedure :: pak_set_stressperiod
147  procedure :: bnd_cq => apt_cq
148  procedure :: bnd_ot_package_flows => apt_ot_package_flows
149  procedure :: bnd_ot_dv => apt_ot_dv
150  procedure :: bnd_ot_bdsummary => apt_ot_bdsummary
151  procedure :: bnd_da => apt_da
152  procedure :: allocate_scalars
154  procedure :: apt_allocate_arrays
155  procedure :: find_apt_package
156  procedure :: apt_solve
157  procedure :: pak_solve
158  procedure :: bnd_options => apt_options
159  procedure :: read_dimensions => apt_read_dimensions
160  procedure :: apt_read_cvs
161  procedure :: read_initial_attr => apt_read_initial_attr
162  procedure :: define_listlabel
163  ! -- methods for observations
164  procedure :: bnd_obs_supported => apt_obs_supported
165  procedure :: bnd_df_obs => apt_df_obs
166  procedure :: pak_df_obs
167  procedure :: pak_rp_obs
168  procedure :: bnd_rp_obs => apt_rp_obs
169  procedure :: rp_obs_byfeature
170  procedure :: rp_obs_budterm
171  procedure :: rp_obs_flowjaface
172  procedure :: bnd_bd_obs => apt_bd_obs
173  procedure :: pak_bd_obs
174  procedure :: get_volumes
175  procedure :: pak_get_nbudterms
176  procedure :: apt_setup_budobj
177  procedure :: pak_setup_budobj
178  procedure :: apt_fill_budobj
179  procedure :: pak_fill_budobj
180  procedure, public :: apt_stor_term
181  procedure, public :: apt_tmvr_term
182  procedure, public :: apt_fmvr_term ! Made public for uze
183  procedure, public :: apt_fjf_term ! Made public for uze
184  procedure, private :: apt_copy2flowp
185  procedure, private :: apt_setup_tableobj
186  procedure, public :: get_mvr_depvar
187 
188  end type tspapttype
189 
190 contains
191 
192  !> @brief Add package connection to matrix
193  !<
194  subroutine apt_ac(this, moffset, sparse)
196  use sparsemodule, only: sparsematrix
197  ! -- dummy
198  class(tspapttype), intent(inout) :: this
199  integer(I4B), intent(in) :: moffset
200  type(sparsematrix), intent(inout) :: sparse
201  ! -- local
202  integer(I4B) :: i, n
203  integer(I4B) :: jj, jglo
204  integer(I4B) :: nglo
205  ! -- format
206  !
207  ! -- Add package rows to sparse
208  if (this%imatrows /= 0) then
209  !
210  ! -- diagonal
211  do n = 1, this%ncv
212  nglo = moffset + this%dis%nodes + this%ioffset + n
213  call sparse%addconnection(nglo, nglo, 1)
214  end do
215  !
216  ! -- apt-gwf connections
217  do i = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
218  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(i)
219  jj = this%flowbudptr%budterm(this%idxbudgwf)%id2(i)
220  nglo = moffset + this%dis%nodes + this%ioffset + n
221  jglo = jj + moffset
222  call sparse%addconnection(nglo, jglo, 1)
223  call sparse%addconnection(jglo, nglo, 1)
224  end do
225  !
226  ! -- apt-apt connections
227  if (this%idxbudfjf /= 0) then
228  do i = 1, this%flowbudptr%budterm(this%idxbudfjf)%maxlist
229  n = this%flowbudptr%budterm(this%idxbudfjf)%id1(i)
230  jj = this%flowbudptr%budterm(this%idxbudfjf)%id2(i)
231  nglo = moffset + this%dis%nodes + this%ioffset + n
232  jglo = moffset + this%dis%nodes + this%ioffset + jj
233  call sparse%addconnection(nglo, jglo, 1)
234  end do
235  end if
236  end if
237  end subroutine apt_ac
238 
239  !> @brief Advanced package transport map package connections to matrix
240  !<
241  subroutine apt_mc(this, moffset, matrix_sln)
242  use sparsemodule, only: sparsematrix
243  ! -- dummy
244  class(tspapttype), intent(inout) :: this
245  integer(I4B), intent(in) :: moffset
246  class(matrixbasetype), pointer :: matrix_sln
247  ! -- local
248  integer(I4B) :: n, j, iglo, jglo
249  integer(I4B) :: ipos
250  ! -- format
251  !
252  ! -- allocate memory for index arrays
253  call this%apt_allocate_index_arrays()
254  !
255  ! -- store index positions
256  if (this%imatrows /= 0) then
257  !
258  ! -- Find the position of each connection in the global ia, ja structure
259  ! and store them in idxglo. idxglo allows this model to insert or
260  ! retrieve values into or from the global A matrix
261  ! -- apt rows
262  do n = 1, this%ncv
263  this%idxlocnode(n) = this%dis%nodes + this%ioffset + n
264  iglo = moffset + this%dis%nodes + this%ioffset + n
265  this%idxpakdiag(n) = matrix_sln%get_position_diag(iglo)
266  end do
267  do ipos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
268  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(ipos)
269  j = this%flowbudptr%budterm(this%idxbudgwf)%id2(ipos)
270  iglo = moffset + this%dis%nodes + this%ioffset + n
271  jglo = j + moffset
272  this%idxdglo(ipos) = matrix_sln%get_position_diag(iglo)
273  this%idxoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
274  end do
275  !
276  ! -- apt contributions to gwf portion of global matrix
277  do ipos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
278  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(ipos)
279  j = this%flowbudptr%budterm(this%idxbudgwf)%id2(ipos)
280  iglo = j + moffset
281  jglo = moffset + this%dis%nodes + this%ioffset + n
282  this%idxsymdglo(ipos) = matrix_sln%get_position_diag(iglo)
283  this%idxsymoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
284  end do
285  !
286  ! -- apt-apt contributions to gwf portion of global matrix
287  if (this%idxbudfjf /= 0) then
288  do ipos = 1, this%flowbudptr%budterm(this%idxbudfjf)%nlist
289  n = this%flowbudptr%budterm(this%idxbudfjf)%id1(ipos)
290  j = this%flowbudptr%budterm(this%idxbudfjf)%id2(ipos)
291  iglo = moffset + this%dis%nodes + this%ioffset + n
292  jglo = moffset + this%dis%nodes + this%ioffset + j
293  this%idxfjfdglo(ipos) = matrix_sln%get_position_diag(iglo)
294  this%idxfjfoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
295  end do
296  end if
297  end if
298  end subroutine apt_mc
299 
300  !> @brief Advanced package transport allocate and read (ar) routine
301  !<
302  subroutine apt_ar(this)
303  ! -- modules
304  ! -- dummy
305  class(tspapttype), intent(inout) :: this
306  ! -- local
307  integer(I4B) :: j
308  logical :: found
309  ! -- formats
310  character(len=*), parameter :: fmtapt = &
311  "(1x,/1x,'APT -- ADVANCED PACKAGE TRANSPORT, VERSION 1, 3/5/2020', &
312  &' INPUT READ FROM UNIT ', i0, //)"
313  !
314  ! -- Get obs setup
315  call this%obs%obs_ar()
316  !
317  ! --print a message identifying the apt package.
318  write (this%iout, fmtapt) this%inunit
319  !
320  ! -- Allocate arrays
321  call this%apt_allocate_arrays()
322  !
323  ! -- read optional initial package parameters
324  call this%read_initial_attr()
325  !
326  ! -- Find the package index in the GWF model or GWF budget file
327  ! for the corresponding apt flow package
328  call this%fmi%get_package_index(this%flowpackagename, this%igwfaptpak)
329  !
330  ! -- Tell fmi that this package is being handled by APT, otherwise
331  ! SSM would handle the flows into GWT from this pack. Then point the
332  ! fmi data for an advanced package to xnewpak and qmfrommvr
333  this%fmi%iatp(this%igwfaptpak) = 1
334  this%fmi%datp(this%igwfaptpak)%concpack => this%get_mvr_depvar()
335  this%fmi%datp(this%igwfaptpak)%qmfrommvr => this%qmfrommvr
336  !
337  ! -- If there is an associated flow package and the user wishes to put
338  ! simulated concentrations (or temperatures) into a aux variable
339  ! column, then find the column number.
340  if (associated(this%flowpackagebnd)) then
341  if (this%cauxfpconc /= '') then
342  found = .false.
343  do j = 1, this%flowpackagebnd%naux
344  if (this%flowpackagebnd%auxname(j) == this%cauxfpconc) then
345  this%iauxfpconc = j
346  found = .true.
347  exit
348  end if
349  end do
350  if (this%iauxfpconc == 0) then
351  errmsg = 'Could not find auxiliary variable '// &
352  trim(adjustl(this%cauxfpconc))//' in flow package '// &
353  trim(adjustl(this%flowpackagename))
354  call store_error(errmsg)
355  call this%parser%StoreErrorUnit()
356  else
357  ! -- tell package not to update this auxiliary variable
358  this%flowpackagebnd%noupdateauxvar(this%iauxfpconc) = 1
359  call this%apt_copy2flowp()
360  end if
361  end if
362  end if
363  end subroutine apt_ar
364 
365  !> @brief Advanced package transport read and prepare (rp) routine
366  !!
367  !! This subroutine calls the attached packages' read and prepare routines.
368  !<
369  subroutine apt_rp(this)
370  use tdismodule, only: kper, nper
371  ! -- dummy
372  class(tspapttype), intent(inout) :: this
373  ! -- local
374  integer(I4B) :: ierr
375  integer(I4B) :: n
376  logical :: isfound, endOfBlock
377  character(len=LINELENGTH) :: title
378  character(len=LINELENGTH) :: line
379  integer(I4B) :: itemno
380  integer(I4B) :: igwfnode
381  ! -- formats
382  character(len=*), parameter :: fmtblkerr = &
383  &"('Error. Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
384  character(len=*), parameter :: fmtlsp = &
385  &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
386  !
387  ! -- set nbound to maxbound
388  this%nbound = this%maxbound
389  !
390  ! -- Set ionper to the stress period number for which a new block of data
391  ! will be read.
392  if (this%inunit == 0) return
393  !
394  ! -- get stress period data
395  if (this%ionper < kper) then
396  !
397  ! -- get period block
398  call this%parser%GetBlock('PERIOD', isfound, ierr, &
399  supportopenclose=.true., &
400  blockrequired=.false.)
401  if (isfound) then
402  !
403  ! -- read ionper and check for increasing period numbers
404  call this%read_check_ionper()
405  else
406  !
407  ! -- PERIOD block not found
408  if (ierr < 0) then
409  ! -- End of file found; data applies for remainder of simulation.
410  this%ionper = nper + 1
411  else
412  ! -- Found invalid block
413  call this%parser%GetCurrentLine(line)
414  write (errmsg, fmtblkerr) adjustl(trim(line))
415  call store_error(errmsg)
416  call this%parser%StoreErrorUnit()
417  end if
418  end if
419  end if
420  !
421  ! -- Read data if ionper == kper
422  if (this%ionper == kper) then
423  !
424  ! -- setup table for period data
425  if (this%iprpak /= 0) then
426  !
427  ! -- reset the input table object
428  title = trim(adjustl(this%text))//' PACKAGE ('// &
429  trim(adjustl(this%packName))//') DATA FOR PERIOD'
430  write (title, '(a,1x,i6)') trim(adjustl(title)), kper
431  call table_cr(this%inputtab, this%packName, title)
432  call this%inputtab%table_df(1, 4, this%iout, finalize=.false.)
433  text = 'NUMBER'
434  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
435  text = 'KEYWORD'
436  call this%inputtab%initialize_column(text, 20, alignment=tableft)
437  do n = 1, 2
438  write (text, '(a,1x,i6)') 'VALUE', n
439  call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
440  end do
441  end if
442  !
443  ! -- read data
444  stressperiod: do
445  call this%parser%GetNextLine(endofblock)
446  if (endofblock) exit
447  !
448  ! -- get feature number
449  itemno = this%parser%GetInteger()
450  !
451  ! -- read data from the rest of the line
452  call this%apt_set_stressperiod(itemno)
453  !
454  ! -- write line to table
455  if (this%iprpak /= 0) then
456  call this%parser%GetCurrentLine(line)
457  call this%inputtab%line_to_columns(line)
458  end if
459  end do stressperiod
460 
461  if (this%iprpak /= 0) then
462  call this%inputtab%finalize_table()
463  end if
464  !
465  ! -- using stress period data from the previous stress period
466  else
467  write (this%iout, fmtlsp) trim(this%filtyp)
468  end if
469  !
470  ! -- write summary of stress period error messages
471  ierr = count_errors()
472  if (ierr > 0) then
473  call this%parser%StoreErrorUnit()
474  end if
475  !
476  ! -- fill arrays
477  do n = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
478  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(n)
479  this%nodelist(n) = igwfnode
480  end do
481  end subroutine apt_rp
482 
483  subroutine apt_ad_chk(this)
484  ! -- dummy
485  class(tspapttype), intent(inout) :: this
486  ! function available for override by packages
487  end subroutine apt_ad_chk
488 
489  !> @brief Advanced package transport set stress period routine.
490  !!
491  !! Set a stress period attribute for an advanced transport package feature
492  !! (itemno) using keywords.
493  !<
494  subroutine apt_set_stressperiod(this, itemno)
495  ! -- module
497  ! -- dummy
498  class(tspapttype), intent(inout) :: this
499  integer(I4B), intent(in) :: itemno
500  ! -- local
501  character(len=LINELENGTH) :: text
502  character(len=LINELENGTH) :: caux
503  character(len=LINELENGTH) :: keyword
504  integer(I4B) :: ierr
505  integer(I4B) :: ii
506  integer(I4B) :: jj
507  real(DP), pointer :: bndElem => null()
508  logical :: found
509  ! -- formats
510  !
511  ! -- Support these general options in LKT, SFT, MWT, UZT
512  ! STATUS <status>
513  ! CONCENTRATION <concentration> or TEMPERATURE <temperature>
514  ! WITHDRAWAL <withdrawal>
515  ! AUXILIARY <auxname> <auxval>
516  !
517  ! -- read line
518  call this%parser%GetStringCaps(keyword)
519  select case (keyword)
520  case ('STATUS')
521  ierr = this%apt_check_valid(itemno)
522  if (ierr /= 0) then
523  goto 999
524  end if
525  call this%parser%GetStringCaps(text)
526  this%status(itemno) = text(1:8)
527  if (text == 'CONSTANT') then
528  this%iboundpak(itemno) = -1
529  else if (text == 'INACTIVE') then
530  this%iboundpak(itemno) = 0
531  else if (text == 'ACTIVE') then
532  this%iboundpak(itemno) = 1
533  else
534  write (errmsg, '(a,a)') &
535  'Unknown '//trim(this%text)//' status keyword: ', text//'.'
536  call store_error(errmsg)
537  end if
538  case ('CONCENTRATION', 'TEMPERATURE')
539  ierr = this%apt_check_valid(itemno)
540  if (ierr /= 0) then
541  goto 999
542  end if
543  call this%parser%GetString(text)
544  jj = 1 ! For feature concentration
545  bndelem => this%concfeat(itemno)
546  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
547  this%packName, 'BND', this%tsManager, &
548  this%iprpak, this%depvartype)
549  case ('AUXILIARY')
550  ierr = this%apt_check_valid(itemno)
551  if (ierr /= 0) then
552  goto 999
553  end if
554  call this%parser%GetStringCaps(caux)
555  do jj = 1, this%naux
556  if (trim(adjustl(caux)) /= trim(adjustl(this%auxname(jj)))) cycle
557  call this%parser%GetString(text)
558  ii = itemno
559  bndelem => this%lauxvar(jj, ii)
560  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
561  this%packName, 'AUX', &
562  this%tsManager, this%iprpak, &
563  this%auxname(jj))
564  exit
565  end do
566  case default
567  !
568  ! -- call the specific package to look for stress period data
569  call this%pak_set_stressperiod(itemno, keyword, found)
570  !
571  ! -- terminate with error if data not valid
572  if (.not. found) then
573  write (errmsg, '(2a)') &
574  'Unknown '//trim(adjustl(this%text))//' data keyword: ', &
575  trim(keyword)//'.'
576  call store_error(errmsg)
577  end if
578  end select
579  !
580  ! -- terminate if any errors were detected
581 999 if (count_errors() > 0) then
582  call this%parser%StoreErrorUnit()
583  end if
584  end subroutine apt_set_stressperiod
585 
586  !> @brief Advanced package transport set stress period routine.
587  !!
588  !! Set a stress period attribute for an individual package. This routine
589  !! must be overridden.
590  !<
591  subroutine pak_set_stressperiod(this, itemno, keyword, found)
592  ! -- dummy
593  class(tspapttype), intent(inout) :: this
594  integer(I4B), intent(in) :: itemno
595  character(len=*), intent(in) :: keyword
596  logical, intent(inout) :: found
597  ! -- local
598 
599  ! -- formats
600  !
601  ! -- this routine should never be called
602  found = .false.
603  call store_error('Program error: pak_set_stressperiod not implemented.', &
604  terminate=.true.)
605  end subroutine pak_set_stressperiod
606 
607  !> @brief Advanced package transport routine
608  !!
609  !! Determine if a valid feature number has been specified.
610  !<
611  function apt_check_valid(this, itemno) result(ierr)
612  ! -- return
613  integer(I4B) :: ierr
614  ! -- dummy
615  class(tspapttype), intent(inout) :: this
616  integer(I4B), intent(in) :: itemno
617  ! -- formats
618  ierr = 0
619  if (itemno < 1 .or. itemno > this%ncv) then
620  write (errmsg, '(a,1x,i6,1x,a,1x,i6)') &
621  'Featureno ', itemno, 'must be > 0 and <= ', this%ncv
622  call store_error(errmsg)
623  ierr = 1
624  end if
625  end function apt_check_valid
626 
627  !> @brief Advanced package transport utility function
628  !!
629  !! Set the concentration (or temperature) to be used by either MVT or MVE
630  !<
631  function get_mvr_depvar(this)
632  ! -- dummy
633  class(tspapttype) :: this
634  ! -- return
635  real(dp), dimension(:), contiguous, pointer :: get_mvr_depvar
636  !
637  get_mvr_depvar => this%xnewpak
638  end function get_mvr_depvar
639 
640  !> @brief Advanced package transport routine
641  !!
642  !! Add package connections to matrix
643  !<
644  subroutine apt_ad(this)
645  ! -- modules
647  ! -- dummy
648  class(tspapttype) :: this
649  ! -- local
650  integer(I4B) :: n
651  integer(I4B) :: j, iaux
652  !
653  ! -- Advance the time series
654  call this%TsManager%ad()
655  !
656  ! -- update auxiliary variables by copying from the derived-type time
657  ! series variable into the bndpackage auxvar variable so that this
658  ! information is properly written to the GWF budget file
659  if (this%naux > 0) then
660  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
661  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
662  do iaux = 1, this%naux
663  this%auxvar(iaux, j) = this%lauxvar(iaux, n)
664  end do
665  end do
666  end if
667  !
668  ! -- copy xnew into xold and set xnewpak to specified concentration (or
669  ! temperature) for constant concentration/temperature features
670  if (ifailedstepretry == 0) then
671  do n = 1, this%ncv
672  this%xoldpak(n) = this%xnewpak(n)
673  if (this%iboundpak(n) < 0) then
674  this%xnewpak(n) = this%concfeat(n)
675  end if
676  end do
677  else
678  do n = 1, this%ncv
679  this%xnewpak(n) = this%xoldpak(n)
680  if (this%iboundpak(n) < 0) then
681  this%xnewpak(n) = this%concfeat(n)
682  end if
683  end do
684  end if
685  !
686  ! -- pakmvrobj ad
687  !if (this%imover == 1) then
688  ! call this%pakmvrobj%ad()
689  !end if
690  !
691  ! -- For each observation, push simulated value and corresponding
692  ! simulation time from "current" to "preceding" and reset
693  ! "current" value.
694  call this%obs%obs_ad()
695  !
696  ! -- run package-specific checks
697  call this%apt_ad_chk()
698  end subroutine apt_ad
699 
700  !> @brief Override bnd reset for custom mover logic
701  subroutine apt_reset(this)
702  class(tspapttype) :: this !< GwtAptType object
703  ! local
704  integer(I4B) :: i
705  !
706  do i = 1, size(this%qmfrommvr)
707  this%qmfrommvr(i) = dzero
708  end do
709  end subroutine apt_reset
710 
711  subroutine apt_fc(this, rhs, ia, idxglo, matrix_sln)
712  ! -- modules
713  ! -- dummy
714  class(tspapttype) :: this
715  real(DP), dimension(:), intent(inout) :: rhs
716  integer(I4B), dimension(:), intent(in) :: ia
717  integer(I4B), dimension(:), intent(in) :: idxglo
718  class(matrixbasetype), pointer :: matrix_sln
719  ! -- local
720  !
721  ! -- Call fc depending on whether or not a matrix is expanded or not
722  if (this%imatrows == 0) then
723  call this%apt_fc_nonexpanded(rhs, ia, idxglo, matrix_sln)
724  else
725  call this%apt_fc_expanded(rhs, ia, idxglo, matrix_sln)
726  end if
727  end subroutine apt_fc
728 
729  !> @brief Advanced package transport fill coefficient (fc) method
730  !!
731  !! Routine to formulate the nonexpanded matrix case in which feature
732  !! concentrations (or temperatures) are solved explicitly
733  !<
734  subroutine apt_fc_nonexpanded(this, rhs, ia, idxglo, matrix_sln)
735  ! -- modules
736  ! -- dummy
737  class(tspapttype) :: this
738  real(DP), dimension(:), intent(inout) :: rhs
739  integer(I4B), dimension(:), intent(in) :: ia
740  integer(I4B), dimension(:), intent(in) :: idxglo
741  class(matrixbasetype), pointer :: matrix_sln
742  ! -- local
743  integer(I4B) :: j, igwfnode, idiag
744  !
745  ! -- solve for concentration (or temperatures) in the features
746  call this%apt_solve()
747  !
748  ! -- add hcof and rhs terms (from apt_solve) to the gwf matrix
749  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
750  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(j)
751  if (this%ibound(igwfnode) < 1) cycle
752  idiag = idxglo(ia(igwfnode))
753  call matrix_sln%add_value_pos(idiag, this%hcof(j))
754  rhs(igwfnode) = rhs(igwfnode) + this%rhs(j)
755  end do
756  end subroutine apt_fc_nonexpanded
757 
758  !> @brief Advanced package transport fill coefficient (fc) method
759  !!
760  !! Routine to formulate the expanded matrix case in which new rows are added
761  !! to the system of equations for each advanced package transport feature
762  !<
763  subroutine apt_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
764  ! -- modules
765  ! -- dummy
766  class(tspapttype) :: this
767  real(DP), dimension(:), intent(inout) :: rhs
768  integer(I4B), dimension(:), intent(in) :: ia
769  integer(I4B), dimension(:), intent(in) :: idxglo
770  class(matrixbasetype), pointer :: matrix_sln
771  ! -- local
772  integer(I4B) :: j, n, n1, n2
773  integer(I4B) :: iloc
774  integer(I4B) :: iposd, iposoffd
775  integer(I4B) :: ipossymd, ipossymoffd
776  real(DP) :: cold
777  real(DP) :: qbnd, qbnd_scaled
778  real(DP) :: omega
779  real(DP) :: rrate
780  real(DP) :: rhsval
781  real(DP) :: hcofval
782  !
783  ! -- call the specific method for the advanced transport package, such as
784  ! what would be overridden by
785  ! GwtLktType, GwtSftType, GwtMwtType, GwtUztType
786  ! This routine will add terms for rainfall, runoff, or other terms
787  ! specific to the package
788  call this%pak_fc_expanded(rhs, ia, idxglo, matrix_sln)
789  !
790  ! -- mass (or energy) storage in features
791  do n = 1, this%ncv
792  cold = this%xoldpak(n)
793  iloc = this%idxlocnode(n)
794  iposd = this%idxpakdiag(n)
795  call this%apt_stor_term(n, n1, n2, rrate, rhsval, hcofval)
796  call matrix_sln%add_value_pos(iposd, hcofval)
797  rhs(iloc) = rhs(iloc) + rhsval
798  end do
799  !
800  ! -- add to mover contribution
801  if (this%idxbudtmvr /= 0) then
802  do j = 1, this%flowbudptr%budterm(this%idxbudtmvr)%nlist
803  call this%apt_tmvr_term(j, n1, n2, rrate, rhsval, hcofval)
804  iloc = this%idxlocnode(n1)
805  iposd = this%idxpakdiag(n1)
806  call matrix_sln%add_value_pos(iposd, hcofval)
807  rhs(iloc) = rhs(iloc) + rhsval
808  end do
809  end if
810  !
811  ! -- add from mover contribution
812  if (this%idxbudfmvr /= 0) then
813  do n = 1, this%ncv
814  rhsval = this%qmfrommvr(n) ! this will already be in terms of energy for heat transport
815  iloc = this%idxlocnode(n)
816  rhs(iloc) = rhs(iloc) - rhsval
817  end do
818  end if
819  !
820  ! -- go through each apt-gwf connection
821  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
822  !
823  ! -- set n to feature number and process if active feature
824  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
825  if (this%iboundpak(n) /= 0) then
826  !
827  ! -- set acoef and rhs to negative so they are relative to apt and not gwt
828  qbnd = this%flowbudptr%budterm(this%idxbudgwf)%flow(j)
829  omega = dzero
830  if (qbnd < dzero) omega = done
831  qbnd_scaled = qbnd * this%eqnsclfac
832  !
833  ! -- add to apt row
834  iposd = this%idxdglo(j)
835  iposoffd = this%idxoffdglo(j)
836  call matrix_sln%add_value_pos(iposd, omega * qbnd_scaled)
837  call matrix_sln%add_value_pos(iposoffd, (done - omega) * qbnd_scaled)
838  !
839  ! -- add to gwf row for apt connection
840  ipossymd = this%idxsymdglo(j)
841  ipossymoffd = this%idxsymoffdglo(j)
842  call matrix_sln%add_value_pos(ipossymd, -(done - omega) * qbnd_scaled)
843  call matrix_sln%add_value_pos(ipossymoffd, -omega * qbnd_scaled)
844  end if
845  end do
846  !
847  ! -- go through each apt-apt connection
848  if (this%idxbudfjf /= 0) then
849  do j = 1, this%flowbudptr%budterm(this%idxbudfjf)%nlist
850  n1 = this%flowbudptr%budterm(this%idxbudfjf)%id1(j)
851  n2 = this%flowbudptr%budterm(this%idxbudfjf)%id2(j)
852  qbnd = this%flowbudptr%budterm(this%idxbudfjf)%flow(j)
853  if (qbnd <= dzero) then
854  omega = done
855  else
856  omega = dzero
857  end if
858  qbnd_scaled = qbnd * this%eqnsclfac
859  iposd = this%idxfjfdglo(j)
860  iposoffd = this%idxfjfoffdglo(j)
861  call matrix_sln%add_value_pos(iposd, omega * qbnd_scaled)
862  call matrix_sln%add_value_pos(iposoffd, (done - omega) * qbnd_scaled)
863  end do
864  end if
865  end subroutine apt_fc_expanded
866 
867  !> @brief Advanced package transport fill coefficient (fc) method
868  !!
869  !! Routine to allow a subclass advanced transport package to inject
870  !! terms into the matrix assembly. This method must be overridden.
871  !<
872  subroutine pak_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
873  ! -- modules
874  ! -- dummy
875  class(tspapttype) :: this
876  real(DP), dimension(:), intent(inout) :: rhs
877  integer(I4B), dimension(:), intent(in) :: ia
878  integer(I4B), dimension(:), intent(in) :: idxglo
879  class(matrixbasetype), pointer :: matrix_sln
880  ! -- local
881  !
882  ! -- this routine should never be called
883  call store_error('Program error: pak_fc_expanded not implemented.', &
884  terminate=.true.)
885  end subroutine pak_fc_expanded
886 
887  !> @brief Advanced package transport routine
888  !!
889  !! Calculate advanced package transport hcof and rhs so transport budget is
890  !! calculated.
891  !<
892  subroutine apt_cfupdate(this)
893  ! -- modules
894  ! -- dummy
895  class(tspapttype) :: this
896  ! -- local
897  integer(I4B) :: j, n
898  real(DP) :: qbnd
899  real(DP) :: omega
900  !
901  ! -- Calculate hcof and rhs terms so GWF exchanges are calculated correctly
902  ! -- go through each apt-gwf connection and calculate
903  ! rhs and hcof terms for gwt/gwe matrix rows
904  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
905  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
906  this%hcof(j) = dzero
907  this%rhs(j) = dzero
908  if (this%iboundpak(n) /= 0) then
909  qbnd = this%flowbudptr%budterm(this%idxbudgwf)%flow(j)
910  omega = dzero
911  if (qbnd < dzero) omega = done
912  this%hcof(j) = -(done - omega) * qbnd * this%eqnsclfac
913  this%rhs(j) = omega * qbnd * this%xnewpak(n) * this%eqnsclfac
914  end if
915  end do
916  end subroutine apt_cfupdate
917 
918  !> @brief Advanced package transport calculate flows (cq) routine
919  !!
920  !! Calculate flows for the advanced package transport feature
921  !<
922  subroutine apt_cq(this, x, flowja, iadv)
923  ! -- modules
924  ! -- dummy
925  class(tspapttype), intent(inout) :: this
926  real(DP), dimension(:), intent(in) :: x
927  real(DP), dimension(:), contiguous, intent(inout) :: flowja
928  integer(I4B), optional, intent(in) :: iadv
929  ! -- local
930  integer(I4B) :: n, n1, n2
931  real(DP) :: rrate
932  !
933  ! -- Solve the feature concentrations (or temperatures) again or update
934  ! the feature hcof and rhs terms
935  if (this%imatrows == 0) then
936  call this%apt_solve()
937  else
938  call this%apt_cfupdate()
939  end if
940  !
941  ! -- call base functionality in bnd_cq
942  call this%BndType%bnd_cq(x, flowja)
943  !
944  ! -- calculate storage term
945  do n = 1, this%ncv
946  rrate = dzero
947  if (this%iboundpak(n) > 0) then
948  call this%apt_stor_term(n, n1, n2, rrate)
949  end if
950  this%qsto(n) = rrate
951  end do
952  !
953  ! -- Copy concentrations (or temperatures) into the flow package auxiliary variable
954  call this%apt_copy2flowp()
955  !
956  ! -- fill the budget object
957  call this%apt_fill_budobj(x, flowja)
958  end subroutine apt_cq
959 
960  !> @brief Save advanced package flows routine
961  !<
962  subroutine apt_ot_package_flows(this, icbcfl, ibudfl)
963  use tdismodule, only: kstp, kper, delt, pertim, totim
964  class(tspapttype) :: this
965  integer(I4B), intent(in) :: icbcfl
966  integer(I4B), intent(in) :: ibudfl
967  integer(I4B) :: ibinun
968  !
969  ! -- write the flows from the budobj
970  ibinun = 0
971  if (this%ibudgetout /= 0) then
972  ibinun = this%ibudgetout
973  end if
974  if (icbcfl == 0) ibinun = 0
975  if (ibinun > 0) then
976  call this%budobj%save_flows(this%dis, ibinun, kstp, kper, delt, &
977  pertim, totim, this%iout)
978  end if
979  !
980  ! -- Print lake flows table
981  if (ibudfl /= 0 .and. this%iprflow /= 0) then
982  call this%budobj%write_flowtable(this%dis, kstp, kper)
983  end if
984  end subroutine apt_ot_package_flows
985 
986  subroutine apt_ot_dv(this, idvsave, idvprint)
987  ! -- modules
988  use constantsmodule, only: lenbudtxt
989  use tdismodule, only: kstp, kper, pertim, totim
991  use inputoutputmodule, only: ulasav
992  ! -- dummy
993  class(tspapttype) :: this
994  integer(I4B), intent(in) :: idvsave
995  integer(I4B), intent(in) :: idvprint
996  ! -- local
997  integer(I4B) :: ibinun
998  integer(I4B) :: n
999  real(DP) :: c
1000  character(len=LENBUDTXT) :: text
1001  !
1002  ! -- set unit number for binary dependent variable output
1003  ibinun = 0
1004  if (this%iconcout /= 0) then
1005  ibinun = this%iconcout
1006  end if
1007  if (idvsave == 0) ibinun = 0
1008  !
1009  ! -- write binary output
1010  if (ibinun > 0) then
1011  do n = 1, this%ncv
1012  c = this%xnewpak(n)
1013  if (this%iboundpak(n) == 0) then
1014  c = dhnoflo
1015  end if
1016  this%dbuff(n) = c
1017  end do
1018  write (text, '(a)') str_pad_left(this%depvartype, lenvarname)
1019  call ulasav(this%dbuff, text, kstp, kper, pertim, totim, &
1020  this%ncv, 1, 1, ibinun)
1021  end if
1022  !
1023  ! -- write apt conc table
1024  if (idvprint /= 0 .and. this%iprconc /= 0) then
1025  !
1026  ! -- set table kstp and kper
1027  call this%dvtab%set_kstpkper(kstp, kper)
1028  !
1029  ! -- fill concentration data
1030  do n = 1, this%ncv
1031  if (this%inamedbound == 1) then
1032  call this%dvtab%add_term(this%featname(n))
1033  end if
1034  call this%dvtab%add_term(n)
1035  call this%dvtab%add_term(this%xnewpak(n))
1036  end do
1037  end if
1038  end subroutine apt_ot_dv
1039 
1040  !> @brief Print advanced package transport dependent variables
1041  !<
1042  subroutine apt_ot_bdsummary(this, kstp, kper, iout, ibudfl)
1043  ! -- module
1044  use tdismodule, only: totim, delt
1045  ! -- dummy
1046  class(tspapttype) :: this !< TspAptType object
1047  integer(I4B), intent(in) :: kstp !< time step number
1048  integer(I4B), intent(in) :: kper !< period number
1049  integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
1050  integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
1051  !
1052  call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
1053  end subroutine apt_ot_bdsummary
1054 
1055  !> @ brief Allocate scalars
1056  !!
1057  !! Allocate scalar variables for an advanced package
1058  !<
1059  subroutine allocate_scalars(this)
1060  ! -- modules
1062  ! -- dummy
1063  class(tspapttype) :: this
1064  ! -- local
1065  !
1066  ! -- allocate scalars in NumericalPackageType
1067  call this%BndType%allocate_scalars()
1068  !
1069  ! -- Allocate
1070  call mem_allocate(this%iauxfpconc, 'IAUXFPCONC', this%memoryPath)
1071  call mem_allocate(this%imatrows, 'IMATROWS', this%memoryPath)
1072  call mem_allocate(this%iprconc, 'IPRCONC', this%memoryPath)
1073  call mem_allocate(this%iconcout, 'ICONCOUT', this%memoryPath)
1074  call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
1075  call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
1076  call mem_allocate(this%igwfaptpak, 'IGWFAPTPAK', this%memoryPath)
1077  call mem_allocate(this%ncv, 'NCV', this%memoryPath)
1078  call mem_allocate(this%idxbudfjf, 'IDXBUDFJF', this%memoryPath)
1079  call mem_allocate(this%idxbudgwf, 'IDXBUDGWF', this%memoryPath)
1080  call mem_allocate(this%idxbudsto, 'IDXBUDSTO', this%memoryPath)
1081  call mem_allocate(this%idxbudtmvr, 'IDXBUDTMVR', this%memoryPath)
1082  call mem_allocate(this%idxbudfmvr, 'IDXBUDFMVR', this%memoryPath)
1083  call mem_allocate(this%idxbudaux, 'IDXBUDAUX', this%memoryPath)
1084  call mem_allocate(this%nconcbudssm, 'NCONCBUDSSM', this%memoryPath)
1085  call mem_allocate(this%idxprepak, 'IDXPREPAK', this%memoryPath)
1086  call mem_allocate(this%idxlastpak, 'IDXLASTPAK', this%memoryPath)
1087  !
1088  ! -- Initialize
1089  this%iauxfpconc = 0
1090  this%imatrows = 1
1091  this%iprconc = 0
1092  this%iconcout = 0
1093  this%ibudgetout = 0
1094  this%ibudcsv = 0
1095  this%igwfaptpak = 0
1096  this%ncv = 0
1097  this%idxbudfjf = 0
1098  this%idxbudgwf = 0
1099  this%idxbudsto = 0
1100  this%idxbudtmvr = 0
1101  this%idxbudfmvr = 0
1102  this%idxbudaux = 0
1103  this%nconcbudssm = 0
1104  this%idxprepak = 0
1105  this%idxlastpak = 0
1106  !
1107  ! -- set this package as causing asymmetric matrix terms
1108  this%iasym = 1
1109  end subroutine allocate_scalars
1110 
1111  !> @ brief Allocate index arrays
1112  !!
1113  !! Allocate arrays that map to locations in the numerical solution
1114  !<
1115  subroutine apt_allocate_index_arrays(this)
1116  ! -- modules
1118  ! -- dummy
1119  class(tspapttype), intent(inout) :: this
1120  ! -- local
1121  integer(I4B) :: n
1122  !
1123  if (this%imatrows /= 0) then
1124  !
1125  ! -- count number of flow-ja-face connections
1126  n = 0
1127  if (this%idxbudfjf /= 0) then
1128  n = this%flowbudptr%budterm(this%idxbudfjf)%maxlist
1129  end if
1130  !
1131  ! -- allocate pointers to global matrix
1132  call mem_allocate(this%idxlocnode, this%ncv, 'IDXLOCNODE', &
1133  this%memoryPath)
1134  call mem_allocate(this%idxpakdiag, this%ncv, 'IDXPAKDIAG', &
1135  this%memoryPath)
1136  call mem_allocate(this%idxdglo, this%maxbound, 'IDXGLO', &
1137  this%memoryPath)
1138  call mem_allocate(this%idxoffdglo, this%maxbound, 'IDXOFFDGLO', &
1139  this%memoryPath)
1140  call mem_allocate(this%idxsymdglo, this%maxbound, 'IDXSYMDGLO', &
1141  this%memoryPath)
1142  call mem_allocate(this%idxsymoffdglo, this%maxbound, 'IDXSYMOFFDGLO', &
1143  this%memoryPath)
1144  call mem_allocate(this%idxfjfdglo, n, 'IDXFJFDGLO', &
1145  this%memoryPath)
1146  call mem_allocate(this%idxfjfoffdglo, n, 'IDXFJFOFFDGLO', &
1147  this%memoryPath)
1148  else
1149  call mem_allocate(this%idxlocnode, 0, 'IDXLOCNODE', &
1150  this%memoryPath)
1151  call mem_allocate(this%idxpakdiag, 0, 'IDXPAKDIAG', &
1152  this%memoryPath)
1153  call mem_allocate(this%idxdglo, 0, 'IDXGLO', &
1154  this%memoryPath)
1155  call mem_allocate(this%idxoffdglo, 0, 'IDXOFFDGLO', &
1156  this%memoryPath)
1157  call mem_allocate(this%idxsymdglo, 0, 'IDXSYMDGLO', &
1158  this%memoryPath)
1159  call mem_allocate(this%idxsymoffdglo, 0, 'IDXSYMOFFDGLO', &
1160  this%memoryPath)
1161  call mem_allocate(this%idxfjfdglo, 0, 'IDXFJFDGLO', &
1162  this%memoryPath)
1163  call mem_allocate(this%idxfjfoffdglo, 0, 'IDXFJFOFFDGLO', &
1164  this%memoryPath)
1165  end if
1166  end subroutine apt_allocate_index_arrays
1167 
1168  !> @ brief Allocate arrays
1169  !!
1170  !! Allocate advanced package transport arrays
1171  !<
1172  subroutine apt_allocate_arrays(this)
1173  ! -- modules
1175  ! -- dummy
1176  class(tspapttype), intent(inout) :: this
1177  ! -- local
1178  integer(I4B) :: n
1179  !
1180  ! -- call standard BndType allocate scalars
1181  call this%BndType%allocate_arrays()
1182  !
1183  ! -- Allocate
1184  !
1185  ! -- allocate and initialize dbuff
1186  if (this%iconcout > 0) then
1187  call mem_allocate(this%dbuff, this%ncv, 'DBUFF', this%memoryPath)
1188  do n = 1, this%ncv
1189  this%dbuff(n) = dzero
1190  end do
1191  else
1192  call mem_allocate(this%dbuff, 0, 'DBUFF', this%memoryPath)
1193  end if
1194  !
1195  ! -- allocate character array for status
1196  allocate (this%status(this%ncv))
1197  !
1198  ! -- time series
1199  call mem_allocate(this%concfeat, this%ncv, 'CONCFEAT', this%memoryPath)
1200  !
1201  ! -- budget terms
1202  call mem_allocate(this%qsto, this%ncv, 'QSTO', this%memoryPath)
1203  call mem_allocate(this%ccterm, this%ncv, 'CCTERM', this%memoryPath)
1204  !
1205  ! -- concentration for budget terms
1206  call mem_allocate(this%concbudssm, this%nconcbudssm, this%ncv, &
1207  'CONCBUDSSM', this%memoryPath)
1208  !
1209  ! -- mass (or energy) added from the mover transport package
1210  call mem_allocate(this%qmfrommvr, this%ncv, 'QMFROMMVR', this%memoryPath)
1211  !
1212  ! -- initialize arrays
1213  do n = 1, this%ncv
1214  this%status(n) = 'ACTIVE'
1215  this%qsto(n) = dzero
1216  this%ccterm(n) = dzero
1217  this%qmfrommvr(n) = dzero
1218  this%concbudssm(:, n) = dzero
1219  this%concfeat(n) = dzero
1220  end do
1221  end subroutine apt_allocate_arrays
1222 
1223  !> @ brief Deallocate memory
1224  !!
1225  !! Deallocate memory associated with this package
1226  !<
1227  subroutine apt_da(this)
1228  ! -- modules
1230  ! -- dummy
1231  class(tspapttype) :: this
1232  ! -- local
1233  !
1234  ! -- deallocate arrays
1235  call mem_deallocate(this%dbuff)
1236  call mem_deallocate(this%qsto)
1237  call mem_deallocate(this%ccterm)
1238  call mem_deallocate(this%strt)
1239  call mem_deallocate(this%ktf)
1240  call mem_deallocate(this%rfeatthk)
1241  call mem_deallocate(this%lauxvar)
1242  call mem_deallocate(this%xoldpak)
1243  if (this%imatrows == 0) then
1244  call mem_deallocate(this%iboundpak)
1245  call mem_deallocate(this%xnewpak)
1246  end if
1247  call mem_deallocate(this%concbudssm)
1248  call mem_deallocate(this%concfeat)
1249  call mem_deallocate(this%qmfrommvr)
1250  deallocate (this%status)
1251  deallocate (this%featname)
1252  !
1253  ! -- budobj
1254  call this%budobj%budgetobject_da()
1255  deallocate (this%budobj)
1256  nullify (this%budobj)
1257  !
1258  ! -- conc table
1259  if (this%iprconc > 0) then
1260  call this%dvtab%table_da()
1261  deallocate (this%dvtab)
1262  nullify (this%dvtab)
1263  end if
1264  !
1265  ! -- index pointers
1266  call mem_deallocate(this%idxlocnode)
1267  call mem_deallocate(this%idxpakdiag)
1268  call mem_deallocate(this%idxdglo)
1269  call mem_deallocate(this%idxoffdglo)
1270  call mem_deallocate(this%idxsymdglo)
1271  call mem_deallocate(this%idxsymoffdglo)
1272  call mem_deallocate(this%idxfjfdglo)
1273  call mem_deallocate(this%idxfjfoffdglo)
1274  !
1275  ! -- deallocate scalars
1276  call mem_deallocate(this%iauxfpconc)
1277  call mem_deallocate(this%imatrows)
1278  call mem_deallocate(this%iprconc)
1279  call mem_deallocate(this%iconcout)
1280  call mem_deallocate(this%ibudgetout)
1281  call mem_deallocate(this%ibudcsv)
1282  call mem_deallocate(this%igwfaptpak)
1283  call mem_deallocate(this%ncv)
1284  call mem_deallocate(this%idxbudfjf)
1285  call mem_deallocate(this%idxbudgwf)
1286  call mem_deallocate(this%idxbudsto)
1287  call mem_deallocate(this%idxbudtmvr)
1288  call mem_deallocate(this%idxbudfmvr)
1289  call mem_deallocate(this%idxbudaux)
1290  call mem_deallocate(this%idxbudssm)
1291  call mem_deallocate(this%nconcbudssm)
1292  call mem_deallocate(this%idxprepak)
1293  call mem_deallocate(this%idxlastpak)
1294  !
1295  ! -- deallocate scalars in NumericalPackageType
1296  call this%BndType%bnd_da()
1297  end subroutine apt_da
1298 
1299  !> @brief Find corresponding advanced package transport package
1300  !<
1301  subroutine find_apt_package(this)
1302  ! -- modules
1304  ! -- dummy
1305  class(tspapttype) :: this
1306  ! -- local
1307  !
1308  ! -- this routine should never be called
1309  call store_error('Program error: pak_solve not implemented.', &
1310  terminate=.true.)
1311  end subroutine find_apt_package
1312 
1313  !> @brief Set options specific to the TspAptType
1314  !!
1315  !! This routine overrides BndType%bnd_options
1316  !<
1317  subroutine apt_options(this, option, found)
1318  use constantsmodule, only: maxcharlen, dzero
1319  use openspecmodule, only: access, form
1321  ! -- dummy
1322  class(tspapttype), intent(inout) :: this
1323  character(len=*), intent(inout) :: option
1324  logical, intent(inout) :: found
1325  ! -- local
1326  character(len=MAXCHARLEN) :: fname, keyword
1327  ! -- formats
1328  character(len=*), parameter :: fmtaptbin = &
1329  "(4x, a, 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, &
1330  &/4x, 'OPENED ON UNIT: ', I0)"
1331  !
1332  found = .true.
1333  select case (option)
1334  case ('FLOW_PACKAGE_NAME')
1335  call this%parser%GetStringCaps(this%flowpackagename)
1336  write (this%iout, '(4x,a)') &
1337  'THIS '//trim(adjustl(this%text))//' PACKAGE CORRESPONDS TO A GWF &
1338  &PACKAGE WITH THE NAME '//trim(adjustl(this%flowpackagename))
1339  case ('FLOW_PACKAGE_AUXILIARY_NAME')
1340  call this%parser%GetStringCaps(this%cauxfpconc)
1341  write (this%iout, '(4x,a)') &
1342  'SIMULATED CONCENTRATIONS WILL BE COPIED INTO THE FLOW PACKAGE &
1343  &AUXILIARY VARIABLE WITH THE NAME '//trim(adjustl(this%cauxfpconc))
1344  case ('DEV_NONEXPANDING_MATRIX')
1345  ! -- use an iterative solution where concentration is not solved
1346  ! as part of the matrix. It is instead solved separately with a
1347  ! general mixing equation and then added to the RHS of the GWT
1348  ! equations
1349  call this%parser%DevOpt()
1350  this%imatrows = 0
1351  write (this%iout, '(4x,a)') &
1352  trim(adjustl(this%text))// &
1353  ' WILL NOT ADD ADDITIONAL ROWS TO THE A MATRIX.'
1354  case ('PRINT_CONCENTRATION', 'PRINT_TEMPERATURE')
1355  this%iprconc = 1
1356  write (this%iout, '(4x,a,1x,a,1x,a)') trim(adjustl(this%text))// &
1357  trim(adjustl(this%depvartype))//'S WILL BE PRINTED TO LISTING &
1358  &FILE.'
1359  case ('CONCENTRATION', 'TEMPERATURE')
1360  call this%parser%GetStringCaps(keyword)
1361  if (keyword == 'FILEOUT') then
1362  call this%parser%GetString(fname)
1363  this%iconcout = getunit()
1364  call openfile(this%iconcout, this%iout, fname, 'DATA(BINARY)', &
1365  form, access, 'REPLACE')
1366  write (this%iout, fmtaptbin) &
1367  trim(adjustl(this%text)), trim(adjustl(this%depvartype)), &
1368  trim(fname), this%iconcout
1369  else
1370  write (errmsg, "('Optional', 1x, a, 1X, 'keyword must &
1371  &be followed by FILEOUT')") this%depvartype
1372  call store_error(errmsg)
1373  end if
1374  case ('BUDGET')
1375  call this%parser%GetStringCaps(keyword)
1376  if (keyword == 'FILEOUT') then
1377  call this%parser%GetString(fname)
1378  call assign_iounit(this%ibudgetout, this%inunit, "BUDGET fileout")
1379  call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
1380  form, access, 'REPLACE')
1381  write (this%iout, fmtaptbin) trim(adjustl(this%text)), 'BUDGET', &
1382  trim(fname), this%ibudgetout
1383  else
1384  call store_error('Optional BUDGET keyword must be followed by FILEOUT')
1385  end if
1386  case ('BUDGETCSV')
1387  call this%parser%GetStringCaps(keyword)
1388  if (keyword == 'FILEOUT') then
1389  call this%parser%GetString(fname)
1390  call assign_iounit(this%ibudcsv, this%inunit, "BUDGETCSV fileout")
1391  call openfile(this%ibudcsv, this%iout, fname, 'CSV', &
1392  filstat_opt='REPLACE')
1393  write (this%iout, fmtaptbin) trim(adjustl(this%text)), 'BUDGET CSV', &
1394  trim(fname), this%ibudcsv
1395  else
1396  call store_error('Optional BUDGETCSV keyword must be followed by &
1397  &FILEOUT')
1398  end if
1399  case default
1400  !
1401  ! -- No options found
1402  found = .false.
1403  end select
1404  end subroutine apt_options
1405 
1406  !> @brief Determine dimensions for this advanced package
1407  !<
1408  subroutine apt_read_dimensions(this)
1409  ! -- dummy
1410  class(tspapttype), intent(inout) :: this
1411  ! -- local
1412  integer(I4B) :: ierr
1413  ! -- format
1414  !
1415  ! -- Set a pointer to the GWF LAK Package budobj
1416  if (this%flowpackagename == '') then
1417  this%flowpackagename = this%packName
1418  write (this%iout, '(4x,a)') &
1419  'THE FLOW PACKAGE NAME FOR '//trim(adjustl(this%text))//' WAS NOT &
1420  &SPECIFIED. SETTING FLOW PACKAGE NAME TO '// &
1421  &trim(adjustl(this%flowpackagename))
1422 
1423  end if
1424  call this%find_apt_package()
1425  !
1426  ! -- Set dimensions from the GWF advanced package
1427  this%ncv = this%flowbudptr%ncv
1428  this%maxbound = this%flowbudptr%budterm(this%idxbudgwf)%maxlist
1429  this%nbound = this%maxbound
1430  write (this%iout, '(a, a)') 'SETTING DIMENSIONS FOR PACKAGE ', this%packName
1431  write (this%iout, '(2x,a,i0)') 'NUMBER OF CONTROL VOLUMES = ', this%ncv
1432  write (this%iout, '(2x,a,i0)') 'MAXBOUND = ', this%maxbound
1433  write (this%iout, '(2x,a,i0)') 'NBOUND = ', this%nbound
1434  if (this%imatrows /= 0) then
1435  this%npakeq = this%ncv
1436  write (this%iout, '(2x,a)') trim(adjustl(this%text))// &
1437  ' SOLVED AS PART OF GWT MATRIX EQUATIONS'
1438  else
1439  write (this%iout, '(2x,a)') trim(adjustl(this%text))// &
1440  ' SOLVED SEPARATELY FROM GWT MATRIX EQUATIONS '
1441  end if
1442  write (this%iout, '(a, //)') 'DONE SETTING DIMENSIONS FOR '// &
1443  trim(adjustl(this%text))
1444  !
1445  ! -- Check for errors
1446  if (this%ncv < 0) then
1447  write (errmsg, '(a)') &
1448  'Number of control volumes could not be determined correctly.'
1449  call store_error(errmsg)
1450  end if
1451  !
1452  ! -- stop if errors were encountered in the DIMENSIONS block
1453  ierr = count_errors()
1454  if (ierr > 0) then
1455  call store_error_unit(this%inunit)
1456  end if
1457  !
1458  ! -- read packagedata block
1459  call this%apt_read_cvs()
1460  !
1461  ! -- Call define_listlabel to construct the list label that is written
1462  ! when PRINT_INPUT option is used.
1463  call this%define_listlabel()
1464  !
1465  ! -- setup the budget object
1466  call this%apt_setup_budobj()
1467  !
1468  ! -- setup the conc table object
1469  call this%apt_setup_tableobj()
1470  end subroutine apt_read_dimensions
1471 
1472  !> @brief Read feature information for this advanced package
1473  !<
1474  subroutine apt_read_cvs(this)
1475  ! -- modules
1478  ! -- dummy
1479  class(tspapttype), intent(inout) :: this
1480  ! -- local
1481  character(len=LINELENGTH) :: text
1482  character(len=LENBOUNDNAME) :: bndName, bndNameTemp
1483  character(len=9) :: cno
1484  character(len=50), dimension(:), allocatable :: caux
1485  integer(I4B) :: ierr
1486  logical :: isfound, endOfBlock
1487  integer(I4B) :: n
1488  integer(I4B) :: ii, jj
1489  integer(I4B) :: iaux
1490  integer(I4B) :: itmp
1491  integer(I4B) :: nlak
1492  integer(I4B) :: nconn
1493  integer(I4B), dimension(:), pointer, contiguous :: nboundchk
1494  real(DP), pointer :: bndElem => null()
1495  !
1496  ! -- initialize itmp
1497  itmp = 0
1498  !
1499  ! -- allocate apt data
1500  call mem_allocate(this%strt, this%ncv, 'STRT', this%memoryPath)
1501  call mem_allocate(this%ktf, this%ncv, 'KTF', this%memoryPath)
1502  call mem_allocate(this%rfeatthk, this%ncv, 'RFEATTHK', this%memoryPath)
1503  call mem_allocate(this%lauxvar, this%naux, this%ncv, 'LAUXVAR', &
1504  this%memoryPath)
1505  !
1506  ! -- lake boundary and concentrations
1507  if (this%imatrows == 0) then
1508  call mem_allocate(this%iboundpak, this%ncv, 'IBOUND', this%memoryPath)
1509  call mem_allocate(this%xnewpak, this%ncv, 'XNEWPAK', this%memoryPath)
1510  end if
1511  call mem_allocate(this%xoldpak, this%ncv, 'XOLDPAK', this%memoryPath)
1512  !
1513  ! -- allocate character storage not managed by the memory manager
1514  allocate (this%featname(this%ncv)) ! ditch after boundnames allocated??
1515  !allocate(this%status(this%ncv))
1516  !
1517  do n = 1, this%ncv
1518  this%strt(n) = dep20
1519  this%ktf(n) = dzero
1520  this%rfeatthk(n) = dzero
1521  this%lauxvar(:, n) = dzero
1522  this%xoldpak(n) = dep20
1523  if (this%imatrows == 0) then
1524  this%iboundpak(n) = 1
1525  this%xnewpak(n) = dep20
1526  end if
1527  end do
1528  !
1529  ! -- allocate local storage for aux variables
1530  if (this%naux > 0) then
1531  allocate (caux(this%naux))
1532  end if
1533  !
1534  ! -- allocate and initialize temporary variables
1535  allocate (nboundchk(this%ncv))
1536  do n = 1, this%ncv
1537  nboundchk(n) = 0
1538  end do
1539  !
1540  ! -- get packagedata block
1541  call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
1542  supportopenclose=.true.)
1543  !
1544  ! -- parse locations block if detected
1545  if (isfound) then
1546  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
1547  ' PACKAGEDATA'
1548  nlak = 0
1549  nconn = 0
1550  do
1551  call this%parser%GetNextLine(endofblock)
1552  if (endofblock) exit
1553  n = this%parser%GetInteger()
1554 
1555  if (n < 1 .or. n > this%ncv) then
1556  write (errmsg, '(a,1x,i6)') &
1557  'Itemno must be > 0 and <= ', this%ncv
1558  call store_error(errmsg)
1559  cycle
1560  end if
1561  !
1562  ! -- increment nboundchk
1563  nboundchk(n) = nboundchk(n) + 1
1564  !
1565  ! -- strt
1566  this%strt(n) = this%parser%GetDouble()
1567  !
1568  ! -- If GWE model, read additional thermal conductivity terms
1569  if (this%depvartype == 'TEMPERATURE') then
1570  ! -- Skip for UZE
1571  if (trim(adjustl(this%text)) /= 'UZE') then
1572  this%ktf(n) = this%parser%GetDouble()
1573  this%rfeatthk(n) = this%parser%GetDouble()
1574  if (this%rfeatthk(n) <= dzero) then
1575  write (errmsg, '(4x,a)') &
1576  '****ERROR. Specified thickness used for thermal &
1577  &conduction MUST BE > 0 else divide by zero error occurs'
1578  call store_error(errmsg)
1579  cycle
1580  end if
1581  end if
1582  end if
1583  !
1584  ! -- get aux data
1585  do iaux = 1, this%naux
1586  call this%parser%GetString(caux(iaux))
1587  end do
1588 
1589  ! -- set default bndName
1590  write (cno, '(i9.9)') n
1591  bndname = 'Feature'//cno
1592 
1593  ! -- featname
1594  if (this%inamedbound /= 0) then
1595  call this%parser%GetStringCaps(bndnametemp)
1596  if (bndnametemp /= '') then
1597  bndname = bndnametemp
1598  end if
1599  end if
1600  this%featname(n) = bndname
1601 
1602  ! -- fill time series aware data
1603  ! -- fill aux data
1604  do jj = 1, this%naux
1605  text = caux(jj)
1606  ii = n
1607  bndelem => this%lauxvar(jj, ii)
1608  call read_value_or_time_series_adv(text, ii, jj, bndelem, &
1609  this%packName, 'AUX', &
1610  this%tsManager, this%iprpak, &
1611  this%auxname(jj))
1612  end do
1613  !
1614  nlak = nlak + 1
1615  end do
1616  !
1617  ! -- check for duplicate or missing lakes
1618  do n = 1, this%ncv
1619  if (nboundchk(n) == 0) then
1620  write (errmsg, '(a,1x,i0)') 'No data specified for feature', n
1621  call store_error(errmsg)
1622  else if (nboundchk(n) > 1) then
1623  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1624  'Data for feature', n, 'specified', nboundchk(n), 'times'
1625  call store_error(errmsg)
1626  end if
1627  end do
1628  !
1629  write (this%iout, '(1x,a)') &
1630  'END OF '//trim(adjustl(this%text))//' PACKAGEDATA'
1631  else
1632  call store_error('Required packagedata block not found.')
1633  end if
1634  !
1635  ! -- terminate if any errors were detected
1636  if (count_errors() > 0) then
1637  call this%parser%StoreErrorUnit()
1638  end if
1639  !
1640  ! -- deallocate local storage for aux variables
1641  if (this%naux > 0) then
1642  deallocate (caux)
1643  end if
1644  !
1645  ! -- deallocate local storage for nboundchk
1646  deallocate (nboundchk)
1647  end subroutine apt_read_cvs
1648 
1649  !> @brief Read the initial parameters for an advanced package
1650  !<
1651  subroutine apt_read_initial_attr(this)
1652  use constantsmodule, only: linelength
1653  use budgetmodule, only: budget_cr
1654  ! -- dummy
1655  class(tspapttype), intent(inout) :: this
1656  ! -- local
1657  !character(len=LINELENGTH) :: text
1658  integer(I4B) :: j, n
1659 
1660  !
1661  ! -- initialize xnewpak and set feature concentration (or temperature)
1662  ! -- todo: this should be a time series?
1663  do n = 1, this%ncv
1664  this%xnewpak(n) = this%strt(n)
1665  !
1666  ! -- todo: read aux
1667  !
1668  ! -- todo: read boundname
1669  end do
1670  !
1671  ! -- initialize status (iboundpak) of lakes to active
1672  do n = 1, this%ncv
1673  if (this%status(n) == 'CONSTANT') then
1674  this%iboundpak(n) = -1
1675  else if (this%status(n) == 'INACTIVE') then
1676  this%iboundpak(n) = 0
1677  else if (this%status(n) == 'ACTIVE ') then
1678  this%iboundpak(n) = 1
1679  end if
1680  end do
1681  !
1682  ! -- set boundname for each connection
1683  if (this%inamedbound /= 0) then
1684  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
1685  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
1686  this%boundname(j) = this%featname(n)
1687  end do
1688  end if
1689  !
1690  ! -- copy boundname into boundname_cst
1691  call this%copy_boundname()
1692  end subroutine apt_read_initial_attr
1693 
1694  !> @brief Add terms specific to advanced package transport to the explicit
1695  !! solve
1696  !!
1697  !! Explicit solve for concentration (or temperature) in advaced package
1698  !! features, which is an alternative to the iterative implicit solve.
1699  !<
1700  subroutine apt_solve(this)
1701  use constantsmodule, only: linelength
1702  ! -- dummy
1703  class(tspapttype) :: this
1704  ! -- local
1705  integer(I4B) :: n, j, igwfnode
1706  integer(I4B) :: n1, n2
1707  real(DP) :: rrate
1708  real(DP) :: ctmp
1709  real(DP) :: c1, qbnd
1710  real(DP) :: hcofval, rhsval
1711  !
1712  ! -- initialize dbuff
1713  do n = 1, this%ncv
1714  this%dbuff(n) = dzero
1715  end do
1716  !
1717  ! -- call the individual package routines to add terms specific to the
1718  ! advanced transport package
1719  call this%pak_solve()
1720  !
1721  ! -- add to mover contribution
1722  if (this%idxbudtmvr /= 0) then
1723  do j = 1, this%flowbudptr%budterm(this%idxbudtmvr)%nlist
1724  call this%apt_tmvr_term(j, n1, n2, rrate)
1725  this%dbuff(n1) = this%dbuff(n1) + rrate
1726  end do
1727  end if
1728  !
1729  ! -- add from mover contribution
1730  if (this%idxbudfmvr /= 0) then
1731  do n1 = 1, size(this%qmfrommvr)
1732  rrate = this%qmfrommvr(n1) ! Will be in terms of energy for heat transport
1733  this%dbuff(n1) = this%dbuff(n1) + rrate
1734  end do
1735  end if
1736  !
1737  ! -- go through each gwf connection and accumulate
1738  ! total mass (or energy) in dbuff mass
1739  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
1740  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
1741  this%hcof(j) = dzero
1742  this%rhs(j) = dzero
1743  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(j)
1744  qbnd = this%flowbudptr%budterm(this%idxbudgwf)%flow(j)
1745  if (qbnd <= dzero) then
1746  ctmp = this%xnewpak(n)
1747  this%rhs(j) = qbnd * ctmp * this%eqnsclfac
1748  else
1749  ctmp = this%xnew(igwfnode)
1750  this%hcof(j) = -qbnd * this%eqnsclfac
1751  end if
1752  c1 = qbnd * ctmp * this%eqnsclfac
1753  this%dbuff(n) = this%dbuff(n) + c1
1754  end do
1755  !
1756  ! -- go through each "within apt-apt" connection (e.g., lak-lak) and
1757  ! accumulate total mass (or energy) in dbuff mass
1758  if (this%idxbudfjf /= 0) then
1759  do j = 1, this%flowbudptr%budterm(this%idxbudfjf)%nlist
1760  call this%apt_fjf_term(j, n1, n2, rrate)
1761  c1 = rrate
1762  this%dbuff(n1) = this%dbuff(n1) + c1
1763  end do
1764  end if
1765  !
1766  ! -- calculate the feature concentration/temperature
1767  do n = 1, this%ncv
1768  call this%apt_stor_term(n, n1, n2, rrate, rhsval, hcofval)
1769  !
1770  ! -- at this point, dbuff has q * c for all sources, so now
1771  ! add Vold / dt * Cold
1772  this%dbuff(n) = this%dbuff(n) - rhsval
1773  !
1774  ! -- Now to calculate c, need to divide dbuff by hcofval
1775  c1 = -this%dbuff(n) / hcofval
1776  if (this%iboundpak(n) > 0) then
1777  this%xnewpak(n) = c1
1778  end if
1779  end do
1780  end subroutine apt_solve
1781 
1782  !> @brief Add terms specific to advanced package transport features to the
1783  !! explicit solve routine
1784  !!
1785  !! This routine must be overridden by the specific apt package
1786  !<
1787  subroutine pak_solve(this)
1788  ! -- dummy
1789  class(tspapttype) :: this
1790  ! -- local
1791  !
1792  ! -- this routine should never be called
1793  call store_error('Program error: pak_solve not implemented.', &
1794  terminate=.true.)
1795  end subroutine pak_solve
1796 
1797  !> @brief Accumulate constant concentration (or temperature) terms for budget
1798  !<
1799  subroutine apt_accumulate_ccterm(this, ilak, rrate, ccratin, ccratout)
1800  ! -- dummy
1801  class(tspapttype) :: this
1802  integer(I4B), intent(in) :: ilak
1803  real(DP), intent(in) :: rrate
1804  real(DP), intent(inout) :: ccratin
1805  real(DP), intent(inout) :: ccratout
1806  ! -- locals
1807  real(DP) :: q
1808  ! format
1809  ! code
1810  !
1811  if (this%iboundpak(ilak) < 0) then
1812  q = -rrate
1813  this%ccterm(ilak) = this%ccterm(ilak) + q
1814  !
1815  ! -- See if flow is into lake or out of lake.
1816  if (q < dzero) then
1817  !
1818  ! -- Flow is out of lake subtract rate from ratout.
1819  ccratout = ccratout - q
1820  else
1821  !
1822  ! -- Flow is into lake; add rate to ratin.
1823  ccratin = ccratin + q
1824  end if
1825  end if
1826  end subroutine apt_accumulate_ccterm
1827 
1828  !> @brief Define the list heading that is written to iout when PRINT_INPUT
1829  !! option is used.
1830  !<
1831  subroutine define_listlabel(this)
1832  class(tspapttype), intent(inout) :: this
1833  !
1834  ! -- create the header list label
1835  this%listlabel = trim(this%filtyp)//' NO.'
1836  if (this%dis%ndim == 3) then
1837  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
1838  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
1839  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
1840  elseif (this%dis%ndim == 2) then
1841  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
1842  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
1843  else
1844  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
1845  end if
1846  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'STRESS RATE'
1847  if (this%inamedbound == 1) then
1848  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
1849  end if
1850  end subroutine define_listlabel
1851 
1852  !> @brief Set pointers to model arrays and variables so that a package has
1853  !! access to these items.
1854  !<
1855  subroutine apt_set_pointers(this, neq, ibound, xnew, xold, flowja)
1856  class(tspapttype) :: this
1857  integer(I4B), pointer :: neq
1858  integer(I4B), dimension(:), pointer, contiguous :: ibound
1859  real(DP), dimension(:), pointer, contiguous :: xnew
1860  real(DP), dimension(:), pointer, contiguous :: xold
1861  real(DP), dimension(:), pointer, contiguous :: flowja
1862  ! -- local
1863  integer(I4B) :: istart, iend
1864  !
1865  ! -- call base BndType set_pointers
1866  call this%BndType%set_pointers(neq, ibound, xnew, xold, flowja)
1867  !
1868  ! -- Set the pointers
1869  !
1870  ! -- set package pointers
1871  if (this%imatrows /= 0) then
1872  istart = this%dis%nodes + this%ioffset + 1
1873  iend = istart + this%ncv - 1
1874  this%iboundpak => this%ibound(istart:iend)
1875  this%xnewpak => this%xnew(istart:iend)
1876  end if
1877  end subroutine apt_set_pointers
1878 
1879  !> @brief Return the feature new volume and old volume
1880  !<
1881  subroutine get_volumes(this, icv, vnew, vold, delt)
1882  ! -- modules
1883  ! -- dummy
1884  class(tspapttype) :: this
1885  integer(I4B), intent(in) :: icv
1886  real(DP), intent(inout) :: vnew, vold
1887  real(DP), intent(in) :: delt
1888  ! -- local
1889  real(DP) :: qss
1890  !
1891  ! -- get volumes
1892  vold = dzero
1893  vnew = vold
1894  if (this%idxbudsto /= 0) then
1895  qss = this%flowbudptr%budterm(this%idxbudsto)%flow(icv)
1896  vnew = this%flowbudptr%budterm(this%idxbudsto)%auxvar(1, icv)
1897  vold = vnew + qss * delt
1898  end if
1899  end subroutine get_volumes
1900 
1901  !> @brief Function to return the number of budget terms just for this package
1902  !!
1903  !! This function must be overridden.
1904  !<
1905  function pak_get_nbudterms(this) result(nbudterms)
1906  ! -- modules
1907  ! -- dummy
1908  class(tspapttype) :: this
1909  ! -- return
1910  integer(I4B) :: nbudterms
1911  ! -- local
1912  !
1913  ! -- this routine should never be called
1914  call store_error('Program error: pak_get_nbudterms not implemented.', &
1915  terminate=.true.)
1916  nbudterms = 0
1917  end function pak_get_nbudterms
1918 
1919  !> @brief Set up the budget object that stores advanced package flow terms
1920  !<
1921  subroutine apt_setup_budobj(this)
1922  ! -- modules
1923  use constantsmodule, only: lenbudtxt
1924  ! -- dummy
1925  class(tspapttype) :: this
1926  ! -- local
1927  integer(I4B) :: nbudterm
1928  integer(I4B) :: nlen
1929  integer(I4B) :: n, n1, n2
1930  integer(I4B) :: maxlist, naux
1931  integer(I4B) :: idx
1932  logical :: ordered_id1
1933  real(DP) :: q
1934  character(len=LENBUDTXT) :: bddim_opt
1935  character(len=LENBUDTXT) :: text, textt
1936  character(len=LENBUDTXT), dimension(1) :: auxtxt
1937  !
1938  ! -- initialize nbudterm
1939  nbudterm = 0
1940  !
1941  ! -- Determine if there are flow-ja-face terms
1942  nlen = 0
1943  if (this%idxbudfjf /= 0) then
1944  nlen = this%flowbudptr%budterm(this%idxbudfjf)%maxlist
1945  end if
1946  !
1947  ! -- Determine the number of budget terms associated with apt.
1948  ! These are fixed for the simulation and cannot change
1949  !
1950  ! -- add one if flow-ja-face present
1951  if (this%idxbudfjf /= 0) nbudterm = nbudterm + 1
1952  !
1953  ! -- All the APT packages have GWF, STORAGE, and CONSTANT
1954  nbudterm = nbudterm + 3
1955  !
1956  ! -- add terms for the specific package
1957  nbudterm = nbudterm + this%pak_get_nbudterms()
1958  !
1959  ! -- add for mover terms and auxiliary
1960  if (this%idxbudtmvr /= 0) nbudterm = nbudterm + 1
1961  if (this%idxbudfmvr /= 0) nbudterm = nbudterm + 1
1962  if (this%naux > 0) nbudterm = nbudterm + 1
1963  !
1964  ! -- set up budobj
1965  call budgetobject_cr(this%budobj, this%packName)
1966  !
1967  bddim_opt = this%depvarunitabbrev
1968  call this%budobj%budgetobject_df(this%ncv, nbudterm, 0, 0, &
1969  bddim_opt=bddim_opt, ibudcsv=this%ibudcsv)
1970  idx = 0
1971  !
1972  ! -- Go through and set up each budget term
1973  if (nlen > 0) then
1974  text = ' FLOW-JA-FACE'
1975  idx = idx + 1
1976  maxlist = this%flowbudptr%budterm(this%idxbudfjf)%maxlist
1977  naux = 0
1978  ordered_id1 = this%flowbudptr%budterm(this%idxbudfjf)%ordered_id1
1979  call this%budobj%budterm(idx)%initialize(text, &
1980  this%name_model, &
1981  this%packName, &
1982  this%name_model, &
1983  this%packName, &
1984  maxlist, .false., .false., &
1985  naux, ordered_id1=ordered_id1)
1986  !
1987  ! -- store outlet connectivity
1988  call this%budobj%budterm(idx)%reset(maxlist)
1989  q = dzero
1990  do n = 1, maxlist
1991  n1 = this%flowbudptr%budterm(this%idxbudfjf)%id1(n)
1992  n2 = this%flowbudptr%budterm(this%idxbudfjf)%id2(n)
1993  call this%budobj%budterm(idx)%update_term(n1, n2, q)
1994  end do
1995  end if
1996  !
1997  ! --
1998  text = ' GWF'
1999  idx = idx + 1
2000  maxlist = this%flowbudptr%budterm(this%idxbudgwf)%maxlist
2001  naux = 0
2002  call this%budobj%budterm(idx)%initialize(text, &
2003  this%name_model, &
2004  this%packName, &
2005  this%name_model, &
2006  this%name_model, &
2007  maxlist, .false., .true., &
2008  naux)
2009  call this%budobj%budterm(idx)%reset(maxlist)
2010  q = dzero
2011  do n = 1, maxlist
2012  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(n)
2013  n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(n)
2014  call this%budobj%budterm(idx)%update_term(n1, n2, q)
2015  end do
2016  !
2017  ! -- Reserve space for the package specific terms
2018  this%idxprepak = idx
2019  call this%pak_setup_budobj(idx)
2020  this%idxlastpak = idx
2021  !
2022  ! --
2023  text = ' STORAGE'
2024  idx = idx + 1
2025  maxlist = this%flowbudptr%budterm(this%idxbudsto)%maxlist
2026  naux = 1
2027  write (textt, '(a)') str_pad_left(this%depvarunit, 16)
2028  auxtxt(1) = textt ! ' MASS' or ' ENERGY'
2029  call this%budobj%budterm(idx)%initialize(text, &
2030  this%name_model, &
2031  this%packName, &
2032  this%name_model, &
2033  this%packName, &
2034  maxlist, .false., .false., &
2035  naux, auxtxt)
2036  if (this%idxbudtmvr /= 0) then
2037  !
2038  ! --
2039  text = ' TO-MVR'
2040  idx = idx + 1
2041  maxlist = this%flowbudptr%budterm(this%idxbudtmvr)%maxlist
2042  naux = 0
2043  ordered_id1 = this%flowbudptr%budterm(this%idxbudtmvr)%ordered_id1
2044  call this%budobj%budterm(idx)%initialize(text, &
2045  this%name_model, &
2046  this%packName, &
2047  this%name_model, &
2048  this%packName, &
2049  maxlist, .false., .false., &
2050  naux, ordered_id1=ordered_id1)
2051  end if
2052  if (this%idxbudfmvr /= 0) then
2053  !
2054  ! --
2055  text = ' FROM-MVR'
2056  idx = idx + 1
2057  maxlist = this%ncv
2058  naux = 0
2059  call this%budobj%budterm(idx)%initialize(text, &
2060  this%name_model, &
2061  this%packName, &
2062  this%name_model, &
2063  this%packName, &
2064  maxlist, .false., .false., &
2065  naux)
2066  end if
2067  !
2068  ! --
2069  text = ' CONSTANT'
2070  idx = idx + 1
2071  maxlist = this%ncv
2072  naux = 0
2073  call this%budobj%budterm(idx)%initialize(text, &
2074  this%name_model, &
2075  this%packName, &
2076  this%name_model, &
2077  this%packName, &
2078  maxlist, .false., .false., &
2079  naux)
2080 
2081  !
2082  ! --
2083  naux = this%naux
2084  if (naux > 0) then
2085  !
2086  ! --
2087  text = ' AUXILIARY'
2088  idx = idx + 1
2089  maxlist = this%ncv
2090  call this%budobj%budterm(idx)%initialize(text, &
2091  this%name_model, &
2092  this%packName, &
2093  this%name_model, &
2094  this%packName, &
2095  maxlist, .false., .false., &
2096  naux, this%auxname)
2097  end if
2098  !
2099  ! -- if flow for each control volume are written to the listing file
2100  if (this%iprflow /= 0) then
2101  call this%budobj%flowtable_df(this%iout)
2102  end if
2103  end subroutine apt_setup_budobj
2104 
2105  !> @brief Set up a budget object that stores an advanced package flows
2106  !!
2107  !! Individual packages set up their budget terms. Must be overridden.
2108  !<
2109  subroutine pak_setup_budobj(this, idx)
2110  ! -- modules
2111  ! -- dummy
2112  class(tspapttype) :: this
2113  integer(I4B), intent(inout) :: idx
2114  ! -- local
2115  !
2116  ! -- this routine should never be called
2117  call store_error('Program error: pak_setup_budobj not implemented.', &
2118  terminate=.true.)
2119  end subroutine pak_setup_budobj
2120 
2121  !> @brief Copy flow terms into this%budobj
2122  !<
2123  subroutine apt_fill_budobj(this, x, flowja)
2124  ! -- modules
2125  use tdismodule, only: delt
2126  ! -- dummy
2127  class(tspapttype) :: this
2128  real(DP), dimension(:), intent(in) :: x
2129  real(DP), dimension(:), contiguous, intent(inout) :: flowja
2130  ! -- local
2131  integer(I4B) :: naux
2132  real(DP), dimension(:), allocatable :: auxvartmp
2133  integer(I4B) :: i, j, n1, n2
2134  integer(I4B) :: idx
2135  integer(I4B) :: nlen
2136  integer(I4B) :: nlist
2137  integer(I4B) :: igwfnode
2138  real(DP) :: q
2139  real(DP) :: v0, v1
2140  real(DP) :: ccratin, ccratout
2141  ! -- formats
2142  !
2143  ! -- initialize counter
2144  idx = 0
2145  !
2146  ! -- initialize ccterm, which is used to sum up all mass (or energy) flows
2147  ! into a constant concentration (or temperature) cell
2148  ccratin = dzero
2149  ccratout = dzero
2150  do n1 = 1, this%ncv
2151  this%ccterm(n1) = dzero
2152  end do
2153  !
2154  ! -- FLOW JA FACE
2155  nlen = 0
2156  if (this%idxbudfjf /= 0) then
2157  nlen = this%flowbudptr%budterm(this%idxbudfjf)%maxlist
2158  end if
2159  if (nlen > 0) then
2160  idx = idx + 1
2161  nlist = this%flowbudptr%budterm(this%idxbudfjf)%maxlist
2162  call this%budobj%budterm(idx)%reset(nlist)
2163  q = dzero
2164  do j = 1, nlist
2165  call this%apt_fjf_term(j, n1, n2, q)
2166  call this%budobj%budterm(idx)%update_term(n1, n2, q)
2167  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
2168  end do
2169  end if
2170  !
2171  ! -- GWF (LEAKAGE)
2172  idx = idx + 1
2173  call this%budobj%budterm(idx)%reset(this%maxbound)
2174  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
2175  q = dzero
2176  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
2177  if (this%iboundpak(n1) /= 0) then
2178  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(j)
2179  q = this%hcof(j) * x(igwfnode) - this%rhs(j)
2180  q = -q ! flip sign so relative to advanced package feature
2181  end if
2182  call this%budobj%budterm(idx)%update_term(n1, igwfnode, q)
2183  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
2184  end do
2185  !
2186  ! -- skip individual package terms for now and process them last
2187  ! -- in case they depend on the other terms (as for uze)
2188  idx = this%idxlastpak
2189  !
2190  ! -- STORAGE
2191  idx = idx + 1
2192  call this%budobj%budterm(idx)%reset(this%ncv)
2193  allocate (auxvartmp(1))
2194  do n1 = 1, this%ncv
2195  call this%get_volumes(n1, v1, v0, delt)
2196  auxvartmp(1) = v1 * this%xnewpak(n1) ! Note: When GWE is added, check if this needs a factor of eqnsclfac
2197  q = this%qsto(n1)
2198  call this%budobj%budterm(idx)%update_term(n1, n1, q, auxvartmp)
2199  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
2200  end do
2201  deallocate (auxvartmp)
2202  !
2203  ! -- TO MOVER
2204  if (this%idxbudtmvr /= 0) then
2205  idx = idx + 1
2206  nlist = this%flowbudptr%budterm(this%idxbudtmvr)%nlist
2207  call this%budobj%budterm(idx)%reset(nlist)
2208  do j = 1, nlist
2209  call this%apt_tmvr_term(j, n1, n2, q)
2210  call this%budobj%budterm(idx)%update_term(n1, n2, q)
2211  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
2212  end do
2213  end if
2214  !
2215  ! -- FROM MOVER
2216  if (this%idxbudfmvr /= 0) then
2217  idx = idx + 1
2218  nlist = this%ncv
2219  call this%budobj%budterm(idx)%reset(nlist)
2220  do j = 1, nlist
2221  call this%apt_fmvr_term(j, n1, n2, q)
2222  call this%budobj%budterm(idx)%update_term(n1, n1, q)
2223  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
2224  end do
2225  end if
2226  !
2227  ! -- CONSTANT FLOW
2228  idx = idx + 1
2229  call this%budobj%budterm(idx)%reset(this%ncv)
2230  do n1 = 1, this%ncv
2231  q = this%ccterm(n1)
2232  call this%budobj%budterm(idx)%update_term(n1, n1, q)
2233  end do
2234  !
2235  ! -- AUXILIARY VARIABLES
2236  naux = this%naux
2237  if (naux > 0) then
2238  idx = idx + 1
2239  allocate (auxvartmp(naux))
2240  call this%budobj%budterm(idx)%reset(this%ncv)
2241  do n1 = 1, this%ncv
2242  q = dzero
2243  do i = 1, naux
2244  auxvartmp(i) = this%lauxvar(i, n1)
2245  end do
2246  call this%budobj%budterm(idx)%update_term(n1, n1, q, auxvartmp)
2247  end do
2248  deallocate (auxvartmp)
2249  end if
2250  !
2251  ! -- individual package terms processed last
2252  idx = this%idxprepak
2253  call this%pak_fill_budobj(idx, x, flowja, ccratin, ccratout)
2254  !
2255  ! --Terms are filled, now accumulate them for this time step
2256  call this%budobj%accumulate_terms()
2257  end subroutine apt_fill_budobj
2258 
2259  !> @brief Copy flow terms into this%budobj, must be overridden
2260  !<
2261  subroutine pak_fill_budobj(this, idx, x, flowja, ccratin, ccratout)
2262  ! -- modules
2263  ! -- dummy
2264  class(tspapttype) :: this
2265  integer(I4B), intent(inout) :: idx
2266  real(DP), dimension(:), intent(in) :: x
2267  real(DP), dimension(:), contiguous, intent(inout) :: flowja
2268  real(DP), intent(inout) :: ccratin
2269  real(DP), intent(inout) :: ccratout
2270  ! -- local
2271  ! -- formats
2272  !
2273  ! -- this routine should never be called
2274  call store_error('Program error: pak_fill_budobj not implemented.', &
2275  terminate=.true.)
2276  end subroutine pak_fill_budobj
2277 
2278  !> @brief Account for mass or energy storage in advanced package features
2279  !<
2280  subroutine apt_stor_term(this, ientry, n1, n2, rrate, &
2281  rhsval, hcofval)
2282  use tdismodule, only: delt
2283  class(tspapttype) :: this
2284  integer(I4B), intent(in) :: ientry
2285  integer(I4B), intent(inout) :: n1
2286  integer(I4B), intent(inout) :: n2
2287  real(DP), intent(inout), optional :: rrate
2288  real(DP), intent(inout), optional :: rhsval
2289  real(DP), intent(inout), optional :: hcofval
2290  real(DP) :: v0, v1
2291  real(DP) :: c0, c1
2292  !
2293  n1 = ientry
2294  n2 = ientry
2295  call this%get_volumes(n1, v1, v0, delt)
2296  c0 = this%xoldpak(n1)
2297  c1 = this%xnewpak(n1)
2298  if (present(rrate)) then
2299  rrate = (-c1 * v1 / delt + c0 * v0 / delt) * this%eqnsclfac
2300  end if
2301  if (present(rhsval)) rhsval = -c0 * v0 * this%eqnsclfac / delt
2302  if (present(hcofval)) hcofval = -v1 * this%eqnsclfac / delt
2303  end subroutine apt_stor_term
2304 
2305  !> @brief Account for mass or energy transferred to the MVR package
2306  !<
2307  subroutine apt_tmvr_term(this, ientry, n1, n2, rrate, &
2308  rhsval, hcofval)
2309  ! -- modules
2310  ! -- dummy
2311  class(tspapttype) :: this
2312  integer(I4B), intent(in) :: ientry
2313  integer(I4B), intent(inout) :: n1
2314  integer(I4B), intent(inout) :: n2
2315  real(DP), intent(inout), optional :: rrate
2316  real(DP), intent(inout), optional :: rhsval
2317  real(DP), intent(inout), optional :: hcofval
2318  ! -- local
2319  real(DP) :: qbnd
2320  real(DP) :: ctmp
2321  !
2322  ! -- Calculate MVR-related terms
2323  n1 = this%flowbudptr%budterm(this%idxbudtmvr)%id1(ientry)
2324  n2 = this%flowbudptr%budterm(this%idxbudtmvr)%id2(ientry)
2325  qbnd = this%flowbudptr%budterm(this%idxbudtmvr)%flow(ientry)
2326  ctmp = this%xnewpak(n1)
2327  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
2328  if (present(rhsval)) rhsval = dzero
2329  if (present(hcofval)) hcofval = qbnd * this%eqnsclfac
2330  end subroutine apt_tmvr_term
2331 
2332  !> @brief Account for mass or energy transferred to this package from the
2333  !! MVR package
2334  !<
2335  subroutine apt_fmvr_term(this, ientry, n1, n2, rrate, &
2336  rhsval, hcofval)
2337  ! -- modules
2338  ! -- dummy
2339  class(tspapttype) :: this
2340  integer(I4B), intent(in) :: ientry
2341  integer(I4B), intent(inout) :: n1
2342  integer(I4B), intent(inout) :: n2
2343  real(DP), intent(inout), optional :: rrate
2344  real(DP), intent(inout), optional :: rhsval
2345  real(DP), intent(inout), optional :: hcofval
2346  !
2347  ! -- Calculate MVR-related terms
2348  n1 = ientry
2349  n2 = n1
2350  if (present(rrate)) rrate = this%qmfrommvr(n1) ! NOTE: When bringing in GWE, ensure this is in terms of energy. Might need to apply eqnsclfac here.
2351  if (present(rhsval)) rhsval = this%qmfrommvr(n1)
2352  if (present(hcofval)) hcofval = dzero
2353  end subroutine apt_fmvr_term
2354 
2355  !> @brief Go through each "within apt-apt" connection (e.g., lkt-lkt, or
2356  !! sft-sft) and accumulate total mass (or energy) in dbuff mass
2357  !<
2358  subroutine apt_fjf_term(this, ientry, n1, n2, rrate, &
2359  rhsval, hcofval)
2360  ! -- modules
2361  ! -- dummy
2362  class(tspapttype) :: this
2363  integer(I4B), intent(in) :: ientry
2364  integer(I4B), intent(inout) :: n1
2365  integer(I4B), intent(inout) :: n2
2366  real(DP), intent(inout), optional :: rrate
2367  real(DP), intent(inout), optional :: rhsval
2368  real(DP), intent(inout), optional :: hcofval
2369  ! -- local
2370  real(DP) :: qbnd
2371  real(DP) :: ctmp
2372  !
2373  n1 = this%flowbudptr%budterm(this%idxbudfjf)%id1(ientry)
2374  n2 = this%flowbudptr%budterm(this%idxbudfjf)%id2(ientry)
2375  qbnd = this%flowbudptr%budterm(this%idxbudfjf)%flow(ientry)
2376  if (qbnd <= 0) then
2377  ctmp = this%xnewpak(n1)
2378  else
2379  ctmp = this%xnewpak(n2)
2380  end if
2381  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
2382  if (present(rhsval)) rhsval = -rrate * this%eqnsclfac
2383  if (present(hcofval)) hcofval = dzero
2384  end subroutine apt_fjf_term
2385 
2386  !> @brief Copy concentrations (or temperatures) into flow package aux
2387  !! variable
2388  !<
2389  subroutine apt_copy2flowp(this)
2390  ! -- modules
2391  ! -- dummy
2392  class(tspapttype) :: this
2393  ! -- local
2394  integer(I4B) :: n, j
2395  !
2396  ! -- copy
2397  if (this%iauxfpconc /= 0) then
2398  !
2399  ! -- go through each apt-gwf connection
2400  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
2401  !
2402  ! -- set n to feature number and process if active feature
2403  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
2404  this%flowpackagebnd%auxvar(this%iauxfpconc, j) = this%xnewpak(n)
2405  end do
2406  end if
2407  end subroutine apt_copy2flowp
2408 
2409  !> @brief Determine whether an obs type is supported
2410  !!
2411  !! This function:
2412  !! - returns true if APT package supports named observation.
2413  !! - overrides BndType%bnd_obs_supported()
2414  !<
2415  logical function apt_obs_supported(this)
2416  ! -- modules
2417  ! -- dummy
2418  class(tspapttype) :: this
2419  !
2420  ! -- Set to true
2421  apt_obs_supported = .true.
2422  end function apt_obs_supported
2423 
2424  !> @brief Define observation type
2425  !!
2426  !! This routine:
2427  !! - stores observation types supported by APT package.
2428  !! - overrides BndType%bnd_df_obs
2429  !<
2430  subroutine apt_df_obs(this)
2431  ! -- modules
2432  ! -- dummy
2433  class(tspapttype) :: this
2434  ! -- local
2435  !
2436  ! -- call additional specific observations for lkt, sft, mwt, and uzt
2437  call this%pak_df_obs()
2438  end subroutine apt_df_obs
2439 
2440  !> @brief Define apt observation type
2441  !!
2442  !! This routine:
2443  !! - stores observations supported by the APT package
2444  !! - must be overridden by child class
2445  subroutine pak_df_obs(this)
2446  ! -- modules
2447  ! -- dummy
2448  class(tspapttype) :: this
2449  ! -- local
2450  !
2451  ! -- this routine should never be called
2452  call store_error('Program error: pak_df_obs not implemented.', &
2453  terminate=.true.)
2454  end subroutine pak_df_obs
2455 
2456  !> @brief Process package specific obs
2457  !!
2458  !! Method to process specific observations for this package.
2459  !<
2460  subroutine pak_rp_obs(this, obsrv, found)
2461  ! -- dummy
2462  class(tspapttype), intent(inout) :: this !< package class
2463  type(observetype), intent(inout) :: obsrv !< observation object
2464  logical, intent(inout) :: found !< indicate whether observation was found
2465  ! -- local
2466  !
2467  ! -- this routine should never be called
2468  call store_error('Program error: pak_rp_obs not implemented.', &
2469  terminate=.true.)
2470  end subroutine pak_rp_obs
2471 
2472  !> @brief Prepare observation
2473  !!
2474  !! Find the indices for this observation assuming they are indexed by
2475  !! feature number
2476  !<
2477  subroutine rp_obs_byfeature(this, obsrv)
2478  class(tspapttype), intent(inout) :: this !< object
2479  type(observetype), intent(inout) :: obsrv !< observation
2480  integer(I4B) :: nn1
2481  integer(I4B) :: j
2482  logical :: jfound
2483  character(len=*), parameter :: fmterr = &
2484  "('Boundary ', a, ' for observation ', a, &
2485  &' is invalid in package ', a)"
2486  nn1 = obsrv%NodeNumber
2487  if (nn1 == namedboundflag) then
2488  jfound = .false.
2489  do j = 1, this%ncv
2490  if (this%featname(j) == obsrv%FeatureName) then
2491  jfound = .true.
2492  call obsrv%AddObsIndex(j)
2493  end if
2494  end do
2495  if (.not. jfound) then
2496  write (errmsg, fmterr) trim(obsrv%FeatureName), trim(obsrv%Name), &
2497  trim(this%packName)
2498  call store_error(errmsg)
2499  end if
2500  else
2501  !
2502  ! -- ensure nn1 is > 0 and < ncv
2503  if (nn1 < 0 .or. nn1 > this%ncv) then
2504  write (errmsg, '(7a, i0, a, i0, a)') &
2505  'Observation ', trim(obsrv%Name), ' of type ', &
2506  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2507  trim(this%packName), ' was assigned ID = ', nn1, &
2508  '. ID must be >= 1 and <= ', this%ncv, '.'
2509  call store_error(errmsg)
2510  end if
2511  call obsrv%AddObsIndex(nn1)
2512  end if
2513  end subroutine rp_obs_byfeature
2514 
2515  !> @brief Prepare observation
2516  !!
2517  !! Find the indices for this observation assuming they are first indexed
2518  !! by feature number and secondly by a connection number
2519  !<
2520  subroutine rp_obs_budterm(this, obsrv, budterm)
2521  class(tspapttype), intent(inout) :: this !< object
2522  type(observetype), intent(inout) :: obsrv !< observation
2523  type(budgettermtype), intent(in) :: budterm !< budget term
2524  integer(I4B) :: nn1
2525  integer(I4B) :: iconn
2526  integer(I4B) :: icv
2527  integer(I4B) :: idx
2528  integer(I4B) :: j
2529  logical :: jfound
2530  character(len=*), parameter :: fmterr = &
2531  "('Boundary ', a, ' for observation ', a, &
2532  &' is invalid in package ', a)"
2533  nn1 = obsrv%NodeNumber
2534  if (nn1 == namedboundflag) then
2535  jfound = .false.
2536  do j = 1, budterm%nlist
2537  icv = budterm%id1(j)
2538  if (this%featname(icv) == obsrv%FeatureName) then
2539  jfound = .true.
2540  call obsrv%AddObsIndex(j)
2541  end if
2542  end do
2543  if (.not. jfound) then
2544  write (errmsg, fmterr) trim(obsrv%FeatureName), trim(obsrv%Name), &
2545  trim(this%packName)
2546  call store_error(errmsg)
2547  end if
2548  else
2549  !
2550  ! -- ensure nn1 is > 0 and < ncv
2551  if (nn1 < 0 .or. nn1 > this%ncv) then
2552  write (errmsg, '(7a, i0, a, i0, a)') &
2553  'Observation ', trim(obsrv%Name), ' of type ', &
2554  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2555  trim(this%packName), ' was assigned ID = ', nn1, &
2556  '. ID must be >= 1 and <= ', this%ncv, '.'
2557  call store_error(errmsg)
2558  end if
2559  iconn = obsrv%NodeNumber2
2560  do j = 1, budterm%nlist
2561  if (budterm%id1(j) == nn1) then
2562  ! -- Look for the first occurrence of nn1, then set indxbnds
2563  ! to the iconn record after that
2564  idx = j + iconn - 1
2565  call obsrv%AddObsIndex(idx)
2566  exit
2567  end if
2568  end do
2569  if (idx < 1 .or. idx > budterm%nlist) then
2570  write (errmsg, '(7a, i0, a, i0, a)') &
2571  'Observation ', trim(obsrv%Name), ' of type ', &
2572  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2573  trim(this%packName), ' specifies iconn = ', iconn, &
2574  ', but this is not a valid connection for ID ', nn1, '.'
2575  call store_error(errmsg)
2576  else if (budterm%id1(idx) /= nn1) then
2577  write (errmsg, '(7a, i0, a, i0, a)') &
2578  'Observation ', trim(obsrv%Name), ' of type ', &
2579  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2580  trim(this%packName), ' specifies iconn = ', iconn, &
2581  ', but this is not a valid connection for ID ', nn1, '.'
2582  call store_error(errmsg)
2583  end if
2584  end if
2585  end subroutine rp_obs_budterm
2586 
2587  !> @brief Prepare observation
2588  !!
2589  !! Find the indices for this observation assuming they are first indexed
2590  !! by a feature number and secondly by a second feature number
2591  !<
2592  subroutine rp_obs_flowjaface(this, obsrv, budterm)
2593  class(tspapttype), intent(inout) :: this !< object
2594  type(observetype), intent(inout) :: obsrv !< observation
2595  type(budgettermtype), intent(in) :: budterm !< budget term
2596  integer(I4B) :: nn1
2597  integer(I4B) :: nn2
2598  integer(I4B) :: icv
2599  integer(I4B) :: j
2600  logical :: jfound
2601  character(len=*), parameter :: fmterr = &
2602  "('Boundary ', a, ' for observation ', a, &
2603  &' is invalid in package ', a)"
2604  nn1 = obsrv%NodeNumber
2605  if (nn1 == namedboundflag) then
2606  jfound = .false.
2607  do j = 1, budterm%nlist
2608  icv = budterm%id1(j)
2609  if (this%featname(icv) == obsrv%FeatureName) then
2610  jfound = .true.
2611  call obsrv%AddObsIndex(j)
2612  end if
2613  end do
2614  if (.not. jfound) then
2615  write (errmsg, fmterr) trim(obsrv%FeatureName), trim(obsrv%Name), &
2616  trim(this%packName)
2617  call store_error(errmsg)
2618  end if
2619  else
2620  !
2621  ! -- ensure nn1 is > 0 and < ncv
2622  if (nn1 < 0 .or. nn1 > this%ncv) then
2623  write (errmsg, '(7a, i0, a, i0, a)') &
2624  'Observation ', trim(obsrv%Name), ' of type ', &
2625  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2626  trim(this%packName), ' was assigned ID = ', nn1, &
2627  '. ID must be >= 1 and <= ', this%ncv, '.'
2628  call store_error(errmsg)
2629  end if
2630  nn2 = obsrv%NodeNumber2
2631  !
2632  ! -- ensure nn2 is > 0 and < ncv
2633  if (nn2 < 0 .or. nn2 > this%ncv) then
2634  write (errmsg, '(7a, i0, a, i0, a)') &
2635  'Observation ', trim(obsrv%Name), ' of type ', &
2636  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2637  trim(this%packName), ' was assigned ID2 = ', nn2, &
2638  '. ID must be >= 1 and <= ', this%ncv, '.'
2639  call store_error(errmsg)
2640  end if
2641  ! -- Look for nn1 and nn2 in id1 and id2
2642  jfound = .false.
2643  do j = 1, budterm%nlist
2644  if (budterm%id1(j) == nn1 .and. budterm%id2(j) == nn2) then
2645  call obsrv%AddObsIndex(j)
2646  jfound = .true.
2647  end if
2648  end do
2649  if (.not. jfound) then
2650  write (errmsg, '(7a, i0, a, i0, a)') &
2651  'Observation ', trim(obsrv%Name), ' of type ', &
2652  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2653  trim(this%packName), &
2654  ' specifies a connection between feature ', nn1, &
2655  ' feature ', nn2, ', but these features are not connected.'
2656  call store_error(errmsg)
2657  end if
2658  end if
2659  end subroutine rp_obs_flowjaface
2660 
2661  !> @brief Read and prepare apt-related observations
2662  !!
2663  !! Method to process specific observations for an apt package
2664  !<
2665  subroutine apt_rp_obs(this)
2666  ! -- modules
2667  use tdismodule, only: kper
2668  ! -- dummy
2669  class(tspapttype), intent(inout) :: this
2670  ! -- local
2671  integer(I4B) :: i
2672  logical :: found
2673  class(observetype), pointer :: obsrv => null()
2674  !
2675  if (kper == 1) then
2676  do i = 1, this%obs%npakobs
2677  obsrv => this%obs%pakobs(i)%obsrv
2678  select case (obsrv%ObsTypeId)
2679  case ('CONCENTRATION', 'TEMPERATURE')
2680  call this%rp_obs_byfeature(obsrv)
2681  !
2682  ! -- catch non-cumulative observation assigned to observation defined
2683  ! by a boundname that is assigned to more than one element
2684  if (obsrv%indxbnds_count > 1) then
2685  write (errmsg, '(a, a, a, a)') &
2686  trim(adjustl(this%depvartype))// &
2687  ' for observation', trim(adjustl(obsrv%Name)), &
2688  ' must be assigned to a feature with a unique boundname.'
2689  call store_error(errmsg)
2690  end if
2691  case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE', 'MWE', 'UZE')
2692  call this%rp_obs_budterm(obsrv, &
2693  this%flowbudptr%budterm(this%idxbudgwf))
2694  case ('FLOW-JA-FACE')
2695  if (this%idxbudfjf > 0) then
2696  call this%rp_obs_flowjaface(obsrv, &
2697  this%flowbudptr%budterm(this%idxbudfjf))
2698  else
2699  write (errmsg, '(7a)') &
2700  'Observation ', trim(obsrv%Name), ' of type ', &
2701  trim(adjustl(obsrv%ObsTypeId)), ' in package ', &
2702  trim(this%packName), &
2703  ' cannot be processed because there are no flow connections.'
2704  call store_error(errmsg)
2705  end if
2706  case ('STORAGE')
2707  call this%rp_obs_byfeature(obsrv)
2708  case ('CONSTANT')
2709  call this%rp_obs_byfeature(obsrv)
2710  case ('FROM-MVR')
2711  call this%rp_obs_byfeature(obsrv)
2712  case default
2713  !
2714  ! -- check the child package for any specific obs
2715  found = .false.
2716  call this%pak_rp_obs(obsrv, found)
2717  !
2718  ! -- if none found then terminate with an error
2719  if (.not. found) then
2720  errmsg = 'Unrecognized observation type "'// &
2721  trim(obsrv%ObsTypeId)//'" for '// &
2722  trim(adjustl(this%text))//' package '// &
2723  trim(this%packName)
2724  call store_error(errmsg, terminate=.true.)
2725  end if
2726  end select
2727 
2728  end do
2729  !
2730  ! -- check for errors
2731  if (count_errors() > 0) then
2732  call store_error_unit(this%obs%inunitobs)
2733  end if
2734  end if
2735  end subroutine apt_rp_obs
2736 
2737  !> @brief Calculate observation values
2738  !!
2739  !! Routine calculates observations common to SFT/LKT/MWT/UZT
2740  !! (or SFE/LKE/MWE/UZE) for as many TspAptType observations that are common
2741  !! among the advanced transport packages
2742  !<
2743  subroutine apt_bd_obs(this)
2744  ! -- modules
2745  ! -- dummy
2746  class(tspapttype) :: this
2747  ! -- local
2748  integer(I4B) :: i
2749  integer(I4B) :: igwfnode
2750  integer(I4B) :: j
2751  integer(I4B) :: jj
2752  integer(I4B) :: n
2753  integer(I4B) :: n1
2754  integer(I4B) :: n2
2755  real(DP) :: v
2756  type(observetype), pointer :: obsrv => null()
2757  logical :: found
2758  !
2759  ! -- Write simulated values for all Advanced Package observations
2760  if (this%obs%npakobs > 0) then
2761  call this%obs%obs_bd_clear()
2762  do i = 1, this%obs%npakobs
2763  obsrv => this%obs%pakobs(i)%obsrv
2764  do j = 1, obsrv%indxbnds_count
2765  v = dnodata
2766  jj = obsrv%indxbnds(j)
2767  select case (obsrv%ObsTypeId)
2768  case ('CONCENTRATION', 'TEMPERATURE')
2769  if (this%iboundpak(jj) /= 0) then
2770  v = this%xnewpak(jj)
2771  end if
2772  case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE', 'MWE', 'UZE')
2773  n = this%flowbudptr%budterm(this%idxbudgwf)%id1(jj)
2774  if (this%iboundpak(n) /= 0) then
2775  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(jj)
2776  v = this%hcof(jj) * this%xnew(igwfnode) - this%rhs(jj)
2777  v = -v
2778  end if
2779  case ('FLOW-JA-FACE')
2780  n = this%flowbudptr%budterm(this%idxbudfjf)%id1(jj)
2781  if (this%iboundpak(n) /= 0) then
2782  call this%apt_fjf_term(jj, n1, n2, v)
2783  end if
2784  case ('STORAGE')
2785  if (this%iboundpak(jj) /= 0) then
2786  v = this%qsto(jj)
2787  end if
2788  case ('CONSTANT')
2789  if (this%iboundpak(jj) /= 0) then
2790  v = this%ccterm(jj)
2791  end if
2792  case ('FROM-MVR')
2793  if (this%iboundpak(jj) /= 0 .and. this%idxbudfmvr > 0) then
2794  call this%apt_fmvr_term(jj, n1, n2, v)
2795  end if
2796  case ('TO-MVR')
2797  if (this%idxbudtmvr > 0) then
2798  n = this%flowbudptr%budterm(this%idxbudtmvr)%id1(jj)
2799  if (this%iboundpak(n) /= 0) then
2800  call this%apt_tmvr_term(jj, n1, n2, v)
2801  end if
2802  end if
2803  case default
2804  found = .false.
2805  !
2806  ! -- check the child package for any specific obs
2807  call this%pak_bd_obs(obsrv%ObsTypeId, jj, v, found)
2808  !
2809  ! -- if none found then terminate with an error
2810  if (.not. found) then
2811  errmsg = 'Unrecognized observation type "'// &
2812  trim(obsrv%ObsTypeId)//'" for '// &
2813  trim(adjustl(this%text))//' package '// &
2814  trim(this%packName)
2815  call store_error(errmsg, terminate=.true.)
2816  end if
2817  end select
2818  call this%obs%SaveOneSimval(obsrv, v)
2819  end do
2820  end do
2821  !
2822  ! -- write summary of error messages
2823  if (count_errors() > 0) then
2824  call store_error_unit(this%obs%inunitobs)
2825  end if
2826  end if
2827  end subroutine apt_bd_obs
2828 
2829  !> @brief Check if observation exists in an advanced package
2830  !<
2831  subroutine pak_bd_obs(this, obstypeid, jj, v, found)
2832  ! -- dummy
2833  class(tspapttype), intent(inout) :: this
2834  character(len=*), intent(in) :: obstypeid
2835  integer(I4B), intent(in) :: jj
2836  real(DP), intent(inout) :: v
2837  logical, intent(inout) :: found
2838  ! -- local
2839  !
2840  ! -- set found = .false. because obstypeid is not known
2841  found = .false.
2842  end subroutine pak_bd_obs
2843 
2844  !> @brief Process observation IDs for an advanced package
2845  !!
2846  !! Method to process observation ID strings for an APT package.
2847  !! This processor is only for observation types that support ID1
2848  !! and not ID2.
2849  !<
2850  subroutine apt_process_obsid(obsrv, dis, inunitobs, iout)
2851  ! -- dummy variables
2852  type(observetype), intent(inout) :: obsrv !< Observation object
2853  class(disbasetype), intent(in) :: dis !< Discretization object
2854  integer(I4B), intent(in) :: inunitobs !< file unit number for the package observation file
2855  integer(I4B), intent(in) :: iout !< model listing file unit number
2856  ! -- local variables
2857  integer(I4B) :: nn1
2858  integer(I4B) :: icol
2859  integer(I4B) :: istart
2860  integer(I4B) :: istop
2861  character(len=LINELENGTH) :: string
2862  character(len=LENBOUNDNAME) :: bndname
2863  !
2864  ! -- initialize local variables
2865  string = obsrv%IDstring
2866  !
2867  ! -- Extract reach number from string and store it.
2868  ! If 1st item is not an integer(I4B), it should be a
2869  ! boundary name--deal with it.
2870  icol = 1
2871  !
2872  ! -- get reach number or boundary name
2873  call extract_idnum_or_bndname(string, icol, istart, istop, nn1, bndname)
2874  if (nn1 == namedboundflag) then
2875  obsrv%FeatureName = bndname
2876  end if
2877  !
2878  ! -- store reach number (NodeNumber)
2879  obsrv%NodeNumber = nn1
2880  !
2881  ! -- store NodeNumber2 as 1 so that this can be used
2882  ! as the iconn value for SFT. This works for SFT
2883  ! because there is only one reach per GWT connection.
2884  obsrv%NodeNumber2 = 1
2885  end subroutine apt_process_obsid
2886 
2887  !> @brief Process observation IDs for a package
2888  !!
2889  !! Method to process observation ID strings for an APT package. This
2890  !! processor is for the case where if ID1 is an integer then ID2 must be
2891  !! provided.
2892  !<
2893  subroutine apt_process_obsid12(obsrv, dis, inunitobs, iout)
2894  ! -- dummy variables
2895  type(observetype), intent(inout) :: obsrv !< Observation object
2896  class(disbasetype), intent(in) :: dis !< Discretization object
2897  integer(I4B), intent(in) :: inunitobs !< file unit number for the package observation file
2898  integer(I4B), intent(in) :: iout !< model listing file unit number
2899  ! -- local variables
2900  integer(I4B) :: nn1
2901  integer(I4B) :: iconn
2902  integer(I4B) :: icol
2903  integer(I4B) :: istart
2904  integer(I4B) :: istop
2905  character(len=LINELENGTH) :: string
2906  character(len=LENBOUNDNAME) :: bndname
2907  !
2908  ! -- initialize local variables
2909  string = obsrv%IDstring
2910  !
2911  ! -- Extract reach number from string and store it.
2912  ! If 1st item is not an integer(I4B), it should be a
2913  ! boundary name--deal with it.
2914  icol = 1
2915  !
2916  ! -- get reach number or boundary name
2917  call extract_idnum_or_bndname(string, icol, istart, istop, nn1, bndname)
2918  if (nn1 == namedboundflag) then
2919  obsrv%FeatureName = bndname
2920  else
2921  call extract_idnum_or_bndname(string, icol, istart, istop, iconn, bndname)
2922  if (len_trim(bndname) < 1 .and. iconn < 0) then
2923  write (errmsg, '(a,1x,a,a,1x,a,1x,a)') &
2924  'For observation type', trim(adjustl(obsrv%ObsTypeId)), &
2925  ', ID given as an integer and not as boundname,', &
2926  'but ID2 is missing. Either change ID to valid', &
2927  'boundname or supply valid entry for ID2.'
2928  call store_error(errmsg)
2929  end if
2930  obsrv%NodeNumber2 = iconn
2931  end if
2932  !
2933  ! -- store reach number (NodeNumber)
2934  obsrv%NodeNumber = nn1
2935  end subroutine apt_process_obsid12
2936 
2937  !> @brief Setup a table object an advanced package
2938  !!
2939  !! Set up the table object that is used to write the apt concentration
2940  !! (or temperature) data. The terms listed here must correspond in the
2941  !! apt_ot method.
2942  !<
2943  subroutine apt_setup_tableobj(this)
2944  ! -- modules
2946  ! -- dummy
2947  class(tspapttype) :: this
2948  ! -- local
2949  integer(I4B) :: nterms
2950  character(len=LINELENGTH) :: title
2951  character(len=LINELENGTH) :: text_temp
2952  !
2953  ! -- setup well head table
2954  if (this%iprconc > 0) then
2955  !
2956  ! -- Determine the number of head table columns
2957  nterms = 2
2958  if (this%inamedbound == 1) nterms = nterms + 1
2959  !
2960  ! -- set up table title
2961  title = trim(adjustl(this%text))//' PACKAGE ('// &
2962  trim(adjustl(this%packName))// &
2963  ') '//trim(adjustl(this%depvartype))// &
2964  &' FOR EACH CONTROL VOLUME'
2965  !
2966  ! -- set up dv tableobj
2967  call table_cr(this%dvtab, this%packName, title)
2968  call this%dvtab%table_df(this%ncv, nterms, this%iout, &
2969  transient=.true.)
2970  !
2971  ! -- Go through and set up table budget term
2972  if (this%inamedbound == 1) then
2973  text_temp = 'NAME'
2974  call this%dvtab%initialize_column(text_temp, 20, alignment=tableft)
2975  end if
2976  !
2977  ! -- feature number
2978  text_temp = 'NUMBER'
2979  call this%dvtab%initialize_column(text_temp, 10, alignment=tabcenter)
2980  !
2981  ! -- feature conc
2982  text_temp = this%depvartype(1:4)
2983  call this%dvtab%initialize_column(text_temp, 12, alignment=tabcenter)
2984  end if
2985  end subroutine apt_setup_tableobj
2986 
2987 end module tspaptmodule
This module contains the base boundary package.
This module contains the BudgetModule.
Definition: Budget.f90:20
subroutine, public budget_cr(this, name_model)
@ brief Create a new budget object
Definition: Budget.f90:84
subroutine, public budgetobject_cr(this, name)
Create a new budget object.
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 dhdry
real dry cell constant
Definition: Constants.f90:94
@ tabcenter
centered table column
Definition: Constants.f90:172
@ tabright
right justified table column
Definition: Constants.f90:173
@ tableft
left justified table column
Definition: Constants.f90:171
@ tabucstring
upper case string table data
Definition: Constants.f90:180
@ tabstring
string table data
Definition: Constants.f90:179
@ tabreal
real table data
Definition: Constants.f90:182
@ tabinteger
integer table data
Definition: Constants.f90:181
integer(i4b), parameter lenpackagename
maximum length of the package name
Definition: Constants.f90:23
integer(i4b), parameter namedboundflag
named bound flag
Definition: Constants.f90:49
real(dp), parameter dnodata
real no data constant
Definition: Constants.f90:95
real(dp), parameter dhnoflo
real no flow constant
Definition: Constants.f90:93
real(dp), parameter dep20
real constant 1e20
Definition: Constants.f90:91
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
integer(i4b), parameter lenftype
maximum length of a package type (DIS, WEL, OC, etc.)
Definition: Constants.f90:39
integer(i4b), parameter lenauxname
maximum length of a aux variable
Definition: Constants.f90:35
integer(i4b), parameter lenboundname
maximum length of a bound name
Definition: Constants.f90:36
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter maxcharlen
maximum length of char string
Definition: Constants.f90:47
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, public assign_iounit(iounit, errunit, description)
@ brief assign io unit number
subroutine, public extract_idnum_or_bndname(line, icol, istart, istop, idnum, bndname)
Starting at position icol, define string as line(istart:istop).
integer(i4b) function, public getunit()
Get a free unit number.
character(len=max(len_trim(str), width)) function, public str_pad_left(str, width)
Function for string manipulation.
subroutine, public ulasav(buf, text, kstp, kper, pertim, totim, ncol, nrow, ilay, ichn)
Save 1 layer array on disk.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
subroutine, public urword(line, icol, istart, istop, ncode, n, r, iout, in)
Extract a word from a string.
This module defines variable data types.
Definition: kind.f90:8
This module contains the derived types ObserveType and ObsDataType.
Definition: Observe.f90:15
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
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_unit(iunit, terminate)
Store the file unit number.
Definition: Sim.f90:168
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
integer(i4b) ifailedstepretry
current retry for this time step
subroutine, public table_cr(this, name, title)
Definition: Table.f90:87
real(dp), pointer, public pertim
time relative to start of stress period
Definition: tdis.f90:30
real(dp), pointer, public totim
time relative to start of simulation
Definition: tdis.f90:32
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
integer(i4b), pointer, public nper
number of stress period
Definition: tdis.f90:21
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 apt_ot_bdsummary(this, kstp, kper, iout, ibudfl)
Print advanced package transport dependent variables.
Definition: tsp-apt.f90:1043
subroutine apt_cfupdate(this)
Advanced package transport routine.
Definition: tsp-apt.f90:893
subroutine pak_setup_budobj(this, idx)
Set up a budget object that stores an advanced package flows.
Definition: tsp-apt.f90:2110
subroutine apt_ar(this)
Advanced package transport allocate and read (ar) routine.
Definition: tsp-apt.f90:303
subroutine pak_fill_budobj(this, idx, x, flowja, ccratin, ccratout)
Copy flow terms into thisbudobj, must be overridden.
Definition: tsp-apt.f90:2262
subroutine rp_obs_flowjaface(this, obsrv, budterm)
Prepare observation.
Definition: tsp-apt.f90:2593
subroutine pak_set_stressperiod(this, itemno, keyword, found)
Advanced package transport set stress period routine.
Definition: tsp-apt.f90:592
subroutine pak_bd_obs(this, obstypeid, jj, v, found)
Check if observation exists in an advanced package.
Definition: tsp-apt.f90:2832
subroutine apt_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
Advanced package transport fill coefficient (fc) method.
Definition: tsp-apt.f90:764
subroutine pak_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
Advanced package transport fill coefficient (fc) method.
Definition: tsp-apt.f90:873
subroutine apt_fjf_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Go through each "within apt-apt" connection (e.g., lkt-lkt, or sft-sft) and accumulate total mass (or...
Definition: tsp-apt.f90:2360
subroutine rp_obs_budterm(this, obsrv, budterm)
Prepare observation.
Definition: tsp-apt.f90:2521
subroutine allocate_scalars(this)
@ brief Allocate scalars
Definition: tsp-apt.f90:1060
subroutine apt_read_initial_attr(this)
Read the initial parameters for an advanced package.
Definition: tsp-apt.f90:1652
subroutine apt_setup_budobj(this)
Set up the budget object that stores advanced package flow terms.
Definition: tsp-apt.f90:1922
subroutine apt_allocate_index_arrays(this)
@ brief Allocate index arrays
Definition: tsp-apt.f90:1116
subroutine apt_rp_obs(this)
Read and prepare apt-related observations.
Definition: tsp-apt.f90:2666
character(len=lenvarname) text
Definition: tsp-apt.f90:63
subroutine apt_ot_package_flows(this, icbcfl, ibudfl)
Save advanced package flows routine.
Definition: tsp-apt.f90:963
real(dp) function, dimension(:), pointer, contiguous get_mvr_depvar(this)
Advanced package transport utility function.
Definition: tsp-apt.f90:632
subroutine get_volumes(this, icv, vnew, vold, delt)
Return the feature new volume and old volume.
Definition: tsp-apt.f90:1882
subroutine apt_allocate_arrays(this)
@ brief Allocate arrays
Definition: tsp-apt.f90:1173
integer(i4b) function pak_get_nbudterms(this)
Function to return the number of budget terms just for this package.
Definition: tsp-apt.f90:1906
subroutine apt_set_stressperiod(this, itemno)
Advanced package transport set stress period routine.
Definition: tsp-apt.f90:495
subroutine apt_df_obs(this)
Define observation type.
Definition: tsp-apt.f90:2431
character(len=lenftype) ftype
Definition: tsp-apt.f90:62
subroutine apt_fmvr_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Account for mass or energy transferred to this package from the MVR package.
Definition: tsp-apt.f90:2337
subroutine apt_da(this)
@ brief Deallocate memory
Definition: tsp-apt.f90:1228
subroutine apt_read_dimensions(this)
Determine dimensions for this advanced package.
Definition: tsp-apt.f90:1409
subroutine apt_mc(this, moffset, matrix_sln)
Advanced package transport map package connections to matrix.
Definition: tsp-apt.f90:242
subroutine apt_fill_budobj(this, x, flowja)
Copy flow terms into thisbudobj.
Definition: tsp-apt.f90:2124
subroutine pak_rp_obs(this, obsrv, found)
Process package specific obs.
Definition: tsp-apt.f90:2461
logical function apt_obs_supported(this)
Determine whether an obs type is supported.
Definition: tsp-apt.f90:2416
subroutine apt_read_cvs(this)
Read feature information for this advanced package.
Definition: tsp-apt.f90:1475
subroutine apt_bd_obs(this)
Calculate observation values.
Definition: tsp-apt.f90:2744
subroutine apt_solve(this)
Add terms specific to advanced package transport to the explicit solve.
Definition: tsp-apt.f90:1701
subroutine apt_tmvr_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Account for mass or energy transferred to the MVR package.
Definition: tsp-apt.f90:2309
subroutine apt_set_pointers(this, neq, ibound, xnew, xold, flowja)
Set pointers to model arrays and variables so that a package has access to these items.
Definition: tsp-apt.f90:1856
subroutine apt_accumulate_ccterm(this, ilak, rrate, ccratin, ccratout)
Accumulate constant concentration (or temperature) terms for budget.
Definition: tsp-apt.f90:1800
subroutine apt_setup_tableobj(this)
Setup a table object an advanced package.
Definition: tsp-apt.f90:2944
subroutine find_apt_package(this)
Find corresponding advanced package transport package.
Definition: tsp-apt.f90:1302
subroutine apt_rp(this)
Advanced package transport read and prepare (rp) routine.
Definition: tsp-apt.f90:370
subroutine apt_stor_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Account for mass or energy storage in advanced package features.
Definition: tsp-apt.f90:2282
subroutine apt_options(this, option, found)
Set options specific to the TspAptType.
Definition: tsp-apt.f90:1318
subroutine rp_obs_byfeature(this, obsrv)
Prepare observation.
Definition: tsp-apt.f90:2478
subroutine apt_copy2flowp(this)
Copy concentrations (or temperatures) into flow package aux variable.
Definition: tsp-apt.f90:2390
subroutine apt_cq(this, x, flowja, iadv)
Advanced package transport calculate flows (cq) routine.
Definition: tsp-apt.f90:923
subroutine apt_ac(this, moffset, sparse)
Add package connection to matrix.
Definition: tsp-apt.f90:195
subroutine, public apt_process_obsid(obsrv, dis, inunitobs, iout)
Process observation IDs for an advanced package.
Definition: tsp-apt.f90:2851
subroutine apt_ad(this)
Advanced package transport routine.
Definition: tsp-apt.f90:645
integer(i4b) function apt_check_valid(this, itemno)
Advanced package transport routine.
Definition: tsp-apt.f90:612
subroutine apt_ad_chk(this)
Definition: tsp-apt.f90:484
subroutine define_listlabel(this)
Define the list heading that is written to iout when PRINT_INPUT option is used.
Definition: tsp-apt.f90:1832
subroutine pak_solve(this)
Add terms specific to advanced package transport features to the explicit solve routine.
Definition: tsp-apt.f90:1788
subroutine apt_fc_nonexpanded(this, rhs, ia, idxglo, matrix_sln)
Advanced package transport fill coefficient (fc) method.
Definition: tsp-apt.f90:735
subroutine apt_reset(this)
Override bnd reset for custom mover logic.
Definition: tsp-apt.f90:702
subroutine apt_fc(this, rhs, ia, idxglo, matrix_sln)
Definition: tsp-apt.f90:712
subroutine apt_ot_dv(this, idvsave, idvprint)
Definition: tsp-apt.f90:987
subroutine pak_df_obs(this)
Define apt observation type.
Definition: tsp-apt.f90:2446
subroutine, public apt_process_obsid12(obsrv, dis, inunitobs, iout)
Process observation IDs for a package.
Definition: tsp-apt.f90:2894
@ brief BndType