MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
gwe-sfe.f90
Go to the documentation of this file.
1 ! -- Stream Energy Transport Module
2 ! -- todo: Temperature decay?
3 ! -- todo: save the sfe temperature into the sfr aux variable? (perhaps needed for GWT-GWE exchanges)
4 ! -- todo: calculate the sfr VISC aux variable using temperature?
5 !
6 ! SFR flows (sfrbudptr) index var SFE term Transport Type
7 !---------------------------------------------------------------------------------
8 
9 ! -- terms from SFR that will be handled by parent APT Package
10 ! FLOW-JA-FACE idxbudfjf FLOW-JA-FACE cv2cv
11 ! GWF (aux FLOW-AREA) idxbudgwf GWF cv2gwf
12 ! STORAGE (aux VOLUME) idxbudsto none used for cv volumes
13 ! FROM-MVR idxbudfmvr FROM-MVR rhow * cpw * q * tmpext = this%qfrommvr(:)
14 ! TO-MVR idxbudtmvr TO-MVR rhow * cpw * q * tfeat
15 
16 ! -- terms from SFR that will be handled by SFE
17 ! RAINFALL idxbudrain RAINFALL rhow * cpw * q * train
18 ! EVAPORATION idxbudevap EVAPORATION rhow * latheatvap * q
19 ! RUNOFF idxbudroff RUNOFF rhow * cpw * q * troff
20 ! EXT-INFLOW idxbudiflw EXT-INFLOW rhow * cpw * q * tiflw
21 ! EXT-OUTFLOW idxbudoutf EXT-OUTFLOW rhow * cpw * q * tfeat
22 
23 ! -- terms not associated with SFR
24 ! STRMBD-COND none STRMBD-COND ctherm * (tcell - tfeat) (ctherm is a thermal conductance for the streambed)
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) dE/dt
32 ! none none AUXILIARY none
33 ! none none CONSTANT accumulate
34 !
35 !
37 
38  use kindmodule, only: dp, i4b
40  use simvariablesmodule, only: errmsg
42  use bndmodule, only: bndtype, getbndfromlist
43  use tspfmimodule, only: tspfmitype
44  use sfrmodule, only: sfrtype
45  use observemodule, only: observetype
50  !
51  implicit none
52  !
53  private
54  public :: sfe_create
55  !
56  character(len=*), parameter :: ftype = 'SFE'
57  character(len=*), parameter :: flowtype = 'SFR'
58  character(len=16) :: text = ' SFE'
59 
60  type, extends(tspapttype) :: gwesfetype
61 
62  type(gweinputdatatype), pointer :: gwecommon => null() !< pointer to shared gwe data used by multiple packages but set in mst
63 
64  integer(I4B), pointer :: idxbudrain => null() !< index of rainfall terms in flowbudptr
65  integer(I4B), pointer :: idxbudevap => null() !< index of evaporation terms in flowbudptr
66  integer(I4B), pointer :: idxbudroff => null() !< index of runoff terms in flowbudptr
67  integer(I4B), pointer :: idxbudiflw => null() !< index of inflow terms in flowbudptr
68  integer(I4B), pointer :: idxbudoutf => null() !< index of outflow terms in flowbudptr
69 
70  real(dp), dimension(:), pointer, contiguous :: temprain => null() !< rainfall temperature
71  real(dp), dimension(:), pointer, contiguous :: tempevap => null() !< evaporation temperature
72  real(dp), dimension(:), pointer, contiguous :: temproff => null() !< runoff temperature
73  real(dp), dimension(:), pointer, contiguous :: tempiflw => null() !< inflow temperature
74  real(dp), dimension(:), pointer, contiguous :: ktf => null() !< thermal conductivity between the sfe and groundwater cell
75  real(dp), dimension(:), pointer, contiguous :: rfeatthk => null() !< thickness of streambed material through which thermal conduction occurs
76 
77  contains
78 
79  procedure :: bnd_da => sfe_da
80  procedure :: allocate_scalars
81  procedure :: apt_allocate_arrays => sfe_allocate_arrays
82  procedure :: find_apt_package => find_sfe_package
83  procedure :: pak_fc_expanded => sfe_fc_expanded
84  procedure :: pak_solve => sfe_solve
85  procedure :: pak_get_nbudterms => sfe_get_nbudterms
86  procedure :: pak_setup_budobj => sfe_setup_budobj
87  procedure :: pak_fill_budobj => sfe_fill_budobj
88  procedure :: sfe_rain_term
89  procedure :: sfe_evap_term
90  procedure :: sfe_roff_term
91  procedure :: sfe_iflw_term
92  procedure :: sfe_outf_term
93  procedure, private :: sfe_sbcd_term
94  procedure :: pak_df_obs => sfe_df_obs
95  procedure :: pak_rp_obs => sfe_rp_obs
96  procedure :: pak_bd_obs => sfe_bd_obs
97  procedure :: pak_set_stressperiod => sfe_set_stressperiod
98  procedure :: apt_read_cvs => sfe_read_cvs
99 
100  end type gwesfetype
101 
102 contains
103 
104  !> @brief Create a new sfe package
105  !<
106  subroutine sfe_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
107  fmi, eqnsclfac, gwecommon, dvt, dvu, dvua)
108  ! -- dummy
109  class(bndtype), pointer :: packobj
110  integer(I4B), intent(in) :: id
111  integer(I4B), intent(in) :: ibcnum
112  integer(I4B), intent(in) :: inunit
113  integer(I4B), intent(in) :: iout
114  character(len=*), intent(in) :: namemodel
115  character(len=*), intent(in) :: pakname
116  type(tspfmitype), pointer :: fmi
117  real(dp), intent(in), pointer :: eqnsclfac !< Governing equation scale factor
118  type(gweinputdatatype), intent(in), target :: gwecommon !< Shared data container for use by multiple GWE packages
119  character(len=*), intent(in) :: dvt !< For GWE, set to "TEMPERATURE" in TspAptType
120  character(len=*), intent(in) :: dvu !< For GWE, set to "energy" in TspAptType
121  character(len=*), intent(in) :: dvua !< For GWE, set to "E" in TspAptType
122  ! -- local
123  type(gwesfetype), pointer :: sfeobj
124  !
125  ! -- Allocate the object and assign values to object variables
126  allocate (sfeobj)
127  packobj => sfeobj
128  !
129  ! -- Create name and memory path
130  call packobj%set_names(ibcnum, namemodel, pakname, ftype)
131  packobj%text = text
132  !
133  ! -- Allocate scalars
134  call sfeobj%allocate_scalars()
135  !
136  ! -- Initialize package
137  call packobj%pack_initialize()
138  !
139  packobj%inunit = inunit
140  packobj%iout = iout
141  packobj%id = id
142  packobj%ibcnum = ibcnum
143  packobj%ncolbnd = 1
144  packobj%iscloc = 1
145  !
146  ! -- Store pointer to flow model interface. When the GwfGwt exchange is
147  ! created, it sets fmi%bndlist so that the GWT model has access to all
148  ! the flow packages
149  sfeobj%fmi => fmi
150  !
151  ! -- Store pointer to governing equation scale factor
152  sfeobj%eqnsclfac => eqnsclfac
153  !
154  ! -- Store pointer to shared data module for accessing cpw, rhow
155  ! for the budget calculations, and for accessing the latent heat of
156  ! vaporization for evaporative cooling.
157  sfeobj%gwecommon => gwecommon
158  !
159  ! -- Set labels that will be used in generalized APT class
160  sfeobj%depvartype = dvt
161  sfeobj%depvarunit = dvu
162  sfeobj%depvarunitabbrev = dvua
163  end subroutine sfe_create
164 
165  !> @brief Find corresponding sfe package
166  !<
167  subroutine find_sfe_package(this)
168  ! -- modules
170  ! -- dummy
171  class(gwesfetype) :: this
172  ! -- local
173  character(len=LINELENGTH) :: errmsg
174  class(bndtype), pointer :: packobj
175  integer(I4B) :: ip, icount
176  integer(I4B) :: nbudterm
177  logical :: found
178  !
179  ! -- Initialize found to false, and error later if flow package cannot
180  ! be found
181  found = .false.
182  !
183  ! -- If user is specifying flows in a binary budget file, then set up
184  ! the budget file reader, otherwise set a pointer to the flow package
185  ! budobj
186  if (this%fmi%flows_from_file) then
187  call this%fmi%set_aptbudobj_pointer(this%flowpackagename, this%flowbudptr)
188  if (associated(this%flowbudptr)) found = .true.
189  !
190  else
191  if (associated(this%fmi%gwfbndlist)) then
192  ! -- Look through gwfbndlist for a flow package with the same name as
193  ! this transport package name
194  do ip = 1, this%fmi%gwfbndlist%Count()
195  packobj => getbndfromlist(this%fmi%gwfbndlist, ip)
196  if (packobj%packName == this%flowpackagename) then
197  found = .true.
198  !
199  ! -- Store BndType pointer to packobj, and then
200  ! use the select type to point to the budobj in flow package
201  this%flowpackagebnd => packobj
202  select type (packobj)
203  type is (sfrtype)
204  this%flowbudptr => packobj%budobj
205  end select
206  end if
207  if (found) exit
208  end do
209  end if
210  end if
211  !
212  ! -- Error if flow package not found
213  if (.not. found) then
214  write (errmsg, '(a)') 'Could not find flow package with name '&
215  &//trim(adjustl(this%flowpackagename))//'.'
216  call store_error(errmsg)
217  call this%parser%StoreErrorUnit()
218  end if
219  !
220  ! -- Allocate space for idxbudssm, which indicates whether this is a
221  ! special budget term or one that is a general source and sink
222  nbudterm = this%flowbudptr%nbudterm
223  call mem_allocate(this%idxbudssm, nbudterm, 'IDXBUDSSM', this%memoryPath)
224  !
225  ! -- Process budget terms and identify special budget terms
226  write (this%iout, '(/, a, a)') &
227  'PROCESSING '//ftype//' INFORMATION FOR ', this%packName
228  write (this%iout, '(a)') ' IDENTIFYING FLOW TERMS IN '//flowtype//' PACKAGE'
229  write (this%iout, '(a, i0)') &
230  ' NUMBER OF '//flowtype//' = ', this%flowbudptr%ncv
231  icount = 1
232  do ip = 1, this%flowbudptr%nbudterm
233  select case (trim(adjustl(this%flowbudptr%budterm(ip)%flowtype)))
234  case ('FLOW-JA-FACE')
235  this%idxbudfjf = ip
236  this%idxbudssm(ip) = 0
237  case ('GWF')
238  this%idxbudgwf = ip
239  this%idxbudssm(ip) = 0
240  case ('STORAGE')
241  this%idxbudsto = ip
242  this%idxbudssm(ip) = 0
243  case ('RAINFALL')
244  this%idxbudrain = ip
245  this%idxbudssm(ip) = 0
246  case ('EVAPORATION')
247  this%idxbudevap = ip
248  this%idxbudssm(ip) = 0
249  case ('RUNOFF')
250  this%idxbudroff = ip
251  this%idxbudssm(ip) = 0
252  case ('EXT-INFLOW')
253  this%idxbudiflw = ip
254  this%idxbudssm(ip) = 0
255  case ('EXT-OUTFLOW')
256  this%idxbudoutf = ip
257  this%idxbudssm(ip) = 0
258  case ('TO-MVR')
259  this%idxbudtmvr = ip
260  this%idxbudssm(ip) = 0
261  case ('FROM-MVR')
262  this%idxbudfmvr = ip
263  this%idxbudssm(ip) = 0
264  case ('AUXILIARY')
265  this%idxbudaux = ip
266  this%idxbudssm(ip) = 0
267  case default
268  !
269  ! -- Set idxbudssm equal to a column index for where the temperatures
270  ! are stored in the concbud(nbudssm, ncv) array
271  this%idxbudssm(ip) = icount
272  icount = icount + 1
273  end select
274  write (this%iout, '(a, i0, " = ", a,/, a, i0)') &
275  ' TERM ', ip, trim(adjustl(this%flowbudptr%budterm(ip)%flowtype)), &
276  ' MAX NO. OF ENTRIES = ', this%flowbudptr%budterm(ip)%maxlist
277  end do
278  write (this%iout, '(a, //)') 'DONE PROCESSING '//ftype//' INFORMATION'
279  end subroutine find_sfe_package
280 
281  !> @brief Add matrix terms related to SFE
282  !!
283  !! This will be called from TspAptType%apt_fc_expanded()
284  !! in order to add matrix terms specifically for SFE
285  !<
286  subroutine sfe_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
287  ! -- dummy
288  class(gwesfetype) :: this
289  real(DP), dimension(:), intent(inout) :: rhs
290  integer(I4B), dimension(:), intent(in) :: ia
291  integer(I4B), dimension(:), intent(in) :: idxglo
292  class(matrixbasetype), pointer :: matrix_sln
293  ! -- local
294  integer(I4B) :: j, n1, n2
295  integer(I4B) :: iloc
296  integer(I4B) :: iposd, iposoffd
297  integer(I4B) :: ipossymd, ipossymoffd
298  real(DP) :: rrate
299  real(DP) :: rhsval
300  real(DP) :: hcofval
301  !
302  ! -- Add rainfall contribution
303  if (this%idxbudrain /= 0) then
304  do j = 1, this%flowbudptr%budterm(this%idxbudrain)%nlist
305  call this%sfe_rain_term(j, n1, n2, rrate, rhsval, hcofval)
306  iloc = this%idxlocnode(n1)
307  iposd = this%idxpakdiag(n1)
308  call matrix_sln%add_value_pos(iposd, hcofval)
309  rhs(iloc) = rhs(iloc) + rhsval
310  end do
311  end if
312  !
313  ! -- Add evaporation contribution
314  if (this%idxbudevap /= 0) then
315  do j = 1, this%flowbudptr%budterm(this%idxbudevap)%nlist
316  call this%sfe_evap_term(j, n1, n2, rrate, rhsval, hcofval)
317  iloc = this%idxlocnode(n1)
318  iposd = this%idxpakdiag(n1)
319  call matrix_sln%add_value_pos(iposd, hcofval)
320  rhs(iloc) = rhs(iloc) + rhsval
321  end do
322  end if
323  !
324  ! -- Add runoff contribution
325  if (this%idxbudroff /= 0) then
326  do j = 1, this%flowbudptr%budterm(this%idxbudroff)%nlist
327  call this%sfe_roff_term(j, n1, n2, rrate, rhsval, hcofval)
328  iloc = this%idxlocnode(n1)
329  iposd = this%idxpakdiag(n1)
330  call matrix_sln%add_value_pos(iposd, hcofval)
331  rhs(iloc) = rhs(iloc) + rhsval
332  end do
333  end if
334  !
335  ! -- Add inflow contribution
336  if (this%idxbudiflw /= 0) then
337  do j = 1, this%flowbudptr%budterm(this%idxbudiflw)%nlist
338  call this%sfe_iflw_term(j, n1, n2, rrate, rhsval, hcofval)
339  iloc = this%idxlocnode(n1)
340  iposd = this%idxpakdiag(n1)
341  call matrix_sln%add_value_pos(iposd, hcofval)
342  rhs(iloc) = rhs(iloc) + rhsval
343  end do
344  end if
345  !
346  ! -- Add outflow contribution
347  if (this%idxbudoutf /= 0) then
348  do j = 1, this%flowbudptr%budterm(this%idxbudoutf)%nlist
349  call this%sfe_outf_term(j, n1, n2, rrate, rhsval, hcofval)
350  iloc = this%idxlocnode(n1)
351  iposd = this%idxpakdiag(n1)
352  call matrix_sln%add_value_pos(iposd, hcofval)
353  rhs(iloc) = rhs(iloc) + rhsval
354  end do
355  end if
356  !
357  ! -- Add streambed conduction contribution
358  do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
359  !
360  ! -- Set n to feature number and process if active feature
361  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
362  if (this%iboundpak(n1) /= 0) then
363  !
364  call this%sfe_sbcd_term(j, n1, n2, rrate, rhsval, hcofval)
365  !
366  ! -- Add to sfe row
367  iposd = this%idxdglo(j)
368  iposoffd = this%idxoffdglo(j)
369  call matrix_sln%add_value_pos(iposd, -hcofval)
370  call matrix_sln%add_value_pos(iposoffd, hcofval)
371  !
372  ! -- Add to gwe row for sfe connection
373  ipossymd = this%idxsymdglo(j)
374  ipossymoffd = this%idxsymoffdglo(j)
375  call matrix_sln%add_value_pos(ipossymd, -hcofval)
376  call matrix_sln%add_value_pos(ipossymoffd, hcofval)
377  end if
378  end do
379  end subroutine sfe_fc_expanded
380 
381  !> @ brief Add terms specific to sfr to the explicit sfe solve
382  !<
383  subroutine sfe_solve(this)
384  ! -- dummy
385  class(gwesfetype) :: this
386  ! -- local
387  integer(I4B) :: j
388  integer(I4B) :: n1, n2
389  real(DP) :: rrate
390  !
391  ! -- Add rainfall contribution
392  if (this%idxbudrain /= 0) then
393  do j = 1, this%flowbudptr%budterm(this%idxbudrain)%nlist
394  call this%sfe_rain_term(j, n1, n2, rrate)
395  this%dbuff(n1) = this%dbuff(n1) + rrate
396  end do
397  end if
398  !
399  ! -- Add evaporation contribution
400  if (this%idxbudevap /= 0) then
401  do j = 1, this%flowbudptr%budterm(this%idxbudevap)%nlist
402  call this%sfe_evap_term(j, n1, n2, rrate)
403  this%dbuff(n1) = this%dbuff(n1) + rrate
404  end do
405  end if
406  !
407  ! -- Add runoff contribution
408  if (this%idxbudroff /= 0) then
409  do j = 1, this%flowbudptr%budterm(this%idxbudroff)%nlist
410  call this%sfe_roff_term(j, n1, n2, rrate)
411  this%dbuff(n1) = this%dbuff(n1) + rrate
412  end do
413  end if
414  !
415  ! -- Add inflow contribution
416  if (this%idxbudiflw /= 0) then
417  do j = 1, this%flowbudptr%budterm(this%idxbudiflw)%nlist
418  call this%sfe_iflw_term(j, n1, n2, rrate)
419  this%dbuff(n1) = this%dbuff(n1) + rrate
420  end do
421  end if
422  !
423  ! -- Add outflow contribution
424  if (this%idxbudoutf /= 0) then
425  do j = 1, this%flowbudptr%budterm(this%idxbudoutf)%nlist
426  call this%sfe_outf_term(j, n1, n2, rrate)
427  this%dbuff(n1) = this%dbuff(n1) + rrate
428  end do
429  end if
430  !
431  ! Note: explicit streambed conduction terms???
432  end subroutine sfe_solve
433 
434  !> @brief Function to return the number of budget terms just for this package.
435  !!
436  !! This overrides a function in the parent class.
437  !<
438  function sfe_get_nbudterms(this) result(nbudterms)
439  ! -- dummy
440  class(gwesfetype) :: this
441  ! -- return
442  integer(I4B) :: nbudterms
443  !
444  ! -- Number of budget terms is 6:
445  ! 1. rainfall
446  ! 2. evaporation
447  ! 3. runoff
448  ! 4. ext-inflow
449  ! 5. ext-outflow
450  ! 6. strmbd-cond
451  nbudterms = 6
452  end function sfe_get_nbudterms
453 
454  !> @brief Set up the budget object that stores all the sfe flows
455  !<
456  subroutine sfe_setup_budobj(this, idx)
457  ! -- modules
458  use constantsmodule, only: lenbudtxt
459  ! -- dummy
460  class(gwesfetype) :: this
461  integer(I4B), intent(inout) :: idx
462  ! -- local
463  integer(I4B) :: n, n1, n2
464  integer(I4B) :: maxlist, naux
465  real(DP) :: q
466  character(len=LENBUDTXT) :: text
467  !
468  ! --
469  text = ' RAINFALL'
470  idx = idx + 1
471  maxlist = this%flowbudptr%budterm(this%idxbudrain)%maxlist
472  naux = 0
473  call this%budobj%budterm(idx)%initialize(text, &
474  this%name_model, &
475  this%packName, &
476  this%name_model, &
477  this%packName, &
478  maxlist, .false., .false., &
479  naux)
480  !
481  ! --
482  text = ' EVAPORATION'
483  idx = idx + 1
484  maxlist = this%flowbudptr%budterm(this%idxbudevap)%maxlist
485  naux = 0
486  call this%budobj%budterm(idx)%initialize(text, &
487  this%name_model, &
488  this%packName, &
489  this%name_model, &
490  this%packName, &
491  maxlist, .false., .false., &
492  naux)
493  !
494  ! --
495  text = ' RUNOFF'
496  idx = idx + 1
497  maxlist = this%flowbudptr%budterm(this%idxbudroff)%maxlist
498  naux = 0
499  call this%budobj%budterm(idx)%initialize(text, &
500  this%name_model, &
501  this%packName, &
502  this%name_model, &
503  this%packName, &
504  maxlist, .false., .false., &
505  naux)
506  !
507  ! --
508  text = ' EXT-INFLOW'
509  idx = idx + 1
510  maxlist = this%flowbudptr%budterm(this%idxbudiflw)%maxlist
511  naux = 0
512  call this%budobj%budterm(idx)%initialize(text, &
513  this%name_model, &
514  this%packName, &
515  this%name_model, &
516  this%packName, &
517  maxlist, .false., .false., &
518  naux)
519  !
520  ! --
521  text = ' EXT-OUTFLOW'
522  idx = idx + 1
523  maxlist = this%flowbudptr%budterm(this%idxbudoutf)%maxlist
524  naux = 0
525  call this%budobj%budterm(idx)%initialize(text, &
526  this%name_model, &
527  this%packName, &
528  this%name_model, &
529  this%packName, &
530  maxlist, .false., .false., &
531  naux)
532  !
533  ! -- Conduction through the wetted streambed
534  text = ' STRMBD-COND'
535  idx = idx + 1
536  maxlist = this%flowbudptr%budterm(this%idxbudgwf)%maxlist
537  naux = 0
538  call this%budobj%budterm(idx)%initialize(text, &
539  this%name_model, &
540  this%packName, &
541  this%name_model, &
542  this%packName, &
543  maxlist, .false., .false., &
544  naux)
545  call this%budobj%budterm(idx)%reset(maxlist)
546  q = dzero
547  do n = 1, maxlist
548  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(n)
549  n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(n)
550  call this%budobj%budterm(idx)%update_term(n1, n2, q)
551  end do
552  end subroutine sfe_setup_budobj
553 
554  !> @brief Copy flow terms into this%budobj
555  !<
556  subroutine sfe_fill_budobj(this, idx, x, flowja, ccratin, ccratout)
557  ! -- dummy
558  class(gwesfetype) :: this
559  integer(I4B), intent(inout) :: idx
560  real(DP), dimension(:), intent(in) :: x
561  real(DP), dimension(:), contiguous, intent(inout) :: flowja
562  real(DP), intent(inout) :: ccratin
563  real(DP), intent(inout) :: ccratout
564  ! -- local
565  integer(I4B) :: j, n1, n2
566  integer(I4B) :: igwfnode
567  integer(I4B) :: nlist
568  integer(I4B) :: idiag
569  real(DP) :: q
570  !
571  ! -- Rain
572  idx = idx + 1
573  nlist = this%flowbudptr%budterm(this%idxbudrain)%nlist
574  call this%budobj%budterm(idx)%reset(nlist)
575  do j = 1, nlist
576  call this%sfe_rain_term(j, n1, n2, q)
577  call this%budobj%budterm(idx)%update_term(n1, n2, q)
578  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
579  end do
580  !
581  ! -- Evaporation
582  idx = idx + 1
583  nlist = this%flowbudptr%budterm(this%idxbudevap)%nlist
584  call this%budobj%budterm(idx)%reset(nlist)
585  do j = 1, nlist
586  call this%sfe_evap_term(j, n1, n2, q)
587  call this%budobj%budterm(idx)%update_term(n1, n2, q)
588  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
589  end do
590  !
591  ! -- Runoff
592  idx = idx + 1
593  nlist = this%flowbudptr%budterm(this%idxbudroff)%nlist
594  call this%budobj%budterm(idx)%reset(nlist)
595  do j = 1, nlist
596  call this%sfe_roff_term(j, n1, n2, q)
597  call this%budobj%budterm(idx)%update_term(n1, n2, q)
598  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
599  end do
600  !
601  ! -- Ext-inflow
602  idx = idx + 1
603  nlist = this%flowbudptr%budterm(this%idxbudiflw)%nlist
604  call this%budobj%budterm(idx)%reset(nlist)
605  do j = 1, nlist
606  call this%sfe_iflw_term(j, n1, n2, q)
607  call this%budobj%budterm(idx)%update_term(n1, n2, q)
608  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
609  end do
610  !
611  ! -- Ext-outflow
612  idx = idx + 1
613  nlist = this%flowbudptr%budterm(this%idxbudoutf)%nlist
614  call this%budobj%budterm(idx)%reset(nlist)
615  do j = 1, nlist
616  call this%sfe_outf_term(j, n1, n2, q)
617  call this%budobj%budterm(idx)%update_term(n1, n2, q)
618  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
619  end do
620  !
621  ! -- Strmbd-cond. The idxbudgwf index is used since it contains the
622  ! sfe/cell mapping information
623  idx = idx + 1
624  nlist = this%flowbudptr%budterm(this%idxbudgwf)%nlist
625  call this%budobj%budterm(idx)%reset(nlist)
626  do j = 1, nlist
627  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(j)
628  if (this%iboundpak(n1) /= 0) then
629  ! -- use igwfnode instead of n2 for consistency with usage in apt;
630  ! helps highlight that cell number is sought and used
631  call this%sfe_sbcd_term(j, n1, igwfnode, q)
632  call this%budobj%budterm(idx)%update_term(n1, igwfnode, q)
633  call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
634  ! -- Contribution to gwe cell budget, which is different than the
635  ! other terms (i.e., ext-inflow, evap, etc.)
636  this%simvals(n1) = this%simvals(n1) - q
637  idiag = this%dis%con%ia(igwfnode)
638  flowja(idiag) = flowja(idiag) - q
639  end if
640  end do
641  end subroutine sfe_fill_budobj
642 
643  !> @brief Allocate scalars specific to the streamflow energy transport (SFE)
644  !! package.
645  !<
646  subroutine allocate_scalars(this)
647  ! -- modules
649  ! -- dummy
650  class(gwesfetype) :: this
651  !
652  ! -- Allocate scalars in TspAptType
653  call this%TspAptType%allocate_scalars()
654  !
655  ! -- Allocate
656  call mem_allocate(this%idxbudrain, 'IDXBUDRAIN', this%memoryPath)
657  call mem_allocate(this%idxbudevap, 'IDXBUDEVAP', this%memoryPath)
658  call mem_allocate(this%idxbudroff, 'IDXBUDROFF', this%memoryPath)
659  call mem_allocate(this%idxbudiflw, 'IDXBUDIFLW', this%memoryPath)
660  call mem_allocate(this%idxbudoutf, 'IDXBUDOUTF', this%memoryPath)
661  !
662  ! -- Initialize
663  this%idxbudrain = 0
664  this%idxbudevap = 0
665  this%idxbudroff = 0
666  this%idxbudiflw = 0
667  this%idxbudoutf = 0
668  end subroutine allocate_scalars
669 
670  !> @brief Allocate arrays specific to the streamflow energy transport (SFE)
671  !! package.
672  !<
673  subroutine sfe_allocate_arrays(this)
674  ! -- modules
676  ! -- dummy
677  class(gwesfetype), intent(inout) :: this
678  ! -- local
679  integer(I4B) :: n
680  !
681  ! -- Time series
682  call mem_allocate(this%temprain, this%ncv, 'TEMPRAIN', this%memoryPath)
683  call mem_allocate(this%tempevap, this%ncv, 'TEMPEVAP', this%memoryPath)
684  call mem_allocate(this%temproff, this%ncv, 'TEMPROFF', this%memoryPath)
685  call mem_allocate(this%tempiflw, this%ncv, 'TEMPIFLW', this%memoryPath)
686  !
687  ! -- Call standard TspAptType allocate arrays
688  call this%TspAptType%apt_allocate_arrays()
689  !
690  ! -- Initialize
691  do n = 1, this%ncv
692  this%temprain(n) = dzero
693  this%tempevap(n) = dzero
694  this%temproff(n) = dzero
695  this%tempiflw(n) = dzero
696  end do
697  end subroutine sfe_allocate_arrays
698 
699  !> @brief Deallocate memory
700  !<
701  subroutine sfe_da(this)
702  ! -- modules
704  ! -- dummy
705  class(gwesfetype) :: this
706  !
707  ! -- Deallocate scalars
708  call mem_deallocate(this%idxbudrain)
709  call mem_deallocate(this%idxbudevap)
710  call mem_deallocate(this%idxbudroff)
711  call mem_deallocate(this%idxbudiflw)
712  call mem_deallocate(this%idxbudoutf)
713  !
714  ! -- Deallocate time series
715  call mem_deallocate(this%temprain)
716  call mem_deallocate(this%tempevap)
717  call mem_deallocate(this%temproff)
718  call mem_deallocate(this%tempiflw)
719  !
720  ! -- Deallocate arrays
721  call mem_deallocate(this%ktf)
722  call mem_deallocate(this%rfeatthk)
723  !
724  ! -- Deallocate scalars in TspAptType
725  call this%TspAptType%bnd_da()
726  end subroutine sfe_da
727 
728  !> @brief Rain term
729  !<
730  subroutine sfe_rain_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
731  ! -- dummy
732  class(gwesfetype) :: this
733  integer(I4B), intent(in) :: ientry
734  integer(I4B), intent(inout) :: n1
735  integer(I4B), intent(inout) :: n2
736  real(DP), intent(inout), optional :: rrate
737  real(DP), intent(inout), optional :: rhsval
738  real(DP), intent(inout), optional :: hcofval
739  ! -- local
740  real(DP) :: qbnd
741  real(DP) :: ctmp
742  !
743  n1 = this%flowbudptr%budterm(this%idxbudrain)%id1(ientry)
744  n2 = this%flowbudptr%budterm(this%idxbudrain)%id2(ientry)
745  qbnd = this%flowbudptr%budterm(this%idxbudrain)%flow(ientry)
746  ctmp = this%temprain(n1)
747  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
748  if (present(rhsval)) rhsval = -rrate
749  if (present(hcofval)) hcofval = dzero
750  end subroutine sfe_rain_term
751 
752  !> @brief Evaporative term
753  !<
754  subroutine sfe_evap_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
755  ! -- dummy
756  class(gwesfetype) :: this
757  integer(I4B), intent(in) :: ientry
758  integer(I4B), intent(inout) :: n1
759  integer(I4B), intent(inout) :: n2
760  real(DP), intent(inout), optional :: rrate
761  real(DP), intent(inout), optional :: rhsval
762  real(DP), intent(inout), optional :: hcofval
763  ! -- local
764  real(DP) :: qbnd
765  real(DP) :: heatlat
766  !
767  n1 = this%flowbudptr%budterm(this%idxbudevap)%id1(ientry)
768  n2 = this%flowbudptr%budterm(this%idxbudevap)%id2(ientry)
769  ! -- note that qbnd is negative for evap
770  qbnd = this%flowbudptr%budterm(this%idxbudevap)%flow(ientry)
771  heatlat = this%gwecommon%gwerhow * this%gwecommon%gwelatheatvap
772  if (present(rrate)) rrate = qbnd * heatlat
773  if (present(rhsval)) rhsval = -rrate
774  if (present(hcofval)) hcofval = dzero
775  end subroutine sfe_evap_term
776 
777  !> @brief Runoff term
778  !<
779  subroutine sfe_roff_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
780  ! -- dummy
781  class(gwesfetype) :: this
782  integer(I4B), intent(in) :: ientry
783  integer(I4B), intent(inout) :: n1
784  integer(I4B), intent(inout) :: n2
785  real(DP), intent(inout), optional :: rrate
786  real(DP), intent(inout), optional :: rhsval
787  real(DP), intent(inout), optional :: hcofval
788  ! -- local
789  real(DP) :: qbnd
790  real(DP) :: ctmp
791  !
792  n1 = this%flowbudptr%budterm(this%idxbudroff)%id1(ientry)
793  n2 = this%flowbudptr%budterm(this%idxbudroff)%id2(ientry)
794  qbnd = this%flowbudptr%budterm(this%idxbudroff)%flow(ientry)
795  ctmp = this%temproff(n1)
796  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
797  if (present(rhsval)) rhsval = -rrate
798  if (present(hcofval)) hcofval = dzero
799  end subroutine sfe_roff_term
800 
801  !> @brief Inflow Term
802  !!
803  !! Accounts for energy added via streamflow entering into a stream channel;
804  !! for example, energy entering the model domain via a specified flow in a
805  !! stream channel.
806  !<
807  subroutine sfe_iflw_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
808  ! -- dummy
809  class(gwesfetype) :: this
810  integer(I4B), intent(in) :: ientry
811  integer(I4B), intent(inout) :: n1
812  integer(I4B), intent(inout) :: n2
813  real(DP), intent(inout), optional :: rrate
814  real(DP), intent(inout), optional :: rhsval
815  real(DP), intent(inout), optional :: hcofval
816  ! -- local
817  real(DP) :: qbnd
818  real(DP) :: ctmp
819  !
820  n1 = this%flowbudptr%budterm(this%idxbudiflw)%id1(ientry)
821  n2 = this%flowbudptr%budterm(this%idxbudiflw)%id2(ientry)
822  qbnd = this%flowbudptr%budterm(this%idxbudiflw)%flow(ientry)
823  ctmp = this%tempiflw(n1)
824  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
825  if (present(rhsval)) rhsval = -rrate
826  if (present(hcofval)) hcofval = dzero
827  end subroutine sfe_iflw_term
828 
829  !> @brief Outflow term
830  !!
831  !! Accounts for the energy leaving a stream channel, for example, energy exiting the
832  !! model domain via a flow in a stream channel flowing out of the active domain.
833  !<
834  subroutine sfe_outf_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
835  ! -- dummy
836  class(gwesfetype) :: this
837  integer(I4B), intent(in) :: ientry
838  integer(I4B), intent(inout) :: n1
839  integer(I4B), intent(inout) :: n2
840  real(DP), intent(inout), optional :: rrate
841  real(DP), intent(inout), optional :: rhsval
842  real(DP), intent(inout), optional :: hcofval
843  ! -- local
844  real(DP) :: qbnd
845  real(DP) :: ctmp
846  !
847  n1 = this%flowbudptr%budterm(this%idxbudoutf)%id1(ientry)
848  n2 = this%flowbudptr%budterm(this%idxbudoutf)%id2(ientry)
849  qbnd = this%flowbudptr%budterm(this%idxbudoutf)%flow(ientry)
850  ctmp = this%xnewpak(n1)
851  if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac
852  if (present(rhsval)) rhsval = dzero
853  if (present(hcofval)) hcofval = qbnd * this%eqnsclfac
854  end subroutine sfe_outf_term
855 
856  !> @brief Streambed conduction term
857  !!
858  !! Accounts for the energy entering or leaving a stream channel by
859  !! thermal conduction through the streambed.
860  !<
861  subroutine sfe_sbcd_term(this, ientry, n1, igwfnode, rrate, rhsval, hcofval)
862  ! -- dummy
863  class(gwesfetype) :: this
864  integer(I4B), intent(in) :: ientry
865  integer(I4B), intent(inout) :: n1
866  integer(I4B), intent(inout) :: igwfnode
867  real(DP), intent(inout), optional :: rrate
868  real(DP), intent(inout), optional :: rhsval
869  real(DP), intent(inout), optional :: hcofval
870  ! -- local
871  integer(I4B) :: auxpos
872  real(DP) :: wa !< wetted area
873  real(DP) :: ktf !< thermal conductivity of streambed material
874  real(DP) :: s !< thickness of conductive streambed material
875  real(DP) :: ctherm
876  !
877  rrate = dzero
878  n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(ientry)
879  ! -- use igwfnode instead of n2 for consistency with usage in apt;
880  ! helps highlight that cell number is sought and used
881  igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(ientry)
882  ! -- For now, there is only 1 aux variable under 'GWF'
883  auxpos = this%flowbudptr%budterm(this%idxbudgwf)%naux
884  wa = this%flowbudptr%budterm(this%idxbudgwf)%auxvar(auxpos, ientry)
885  ktf = this%ktf(n1)
886  s = this%rfeatthk(n1)
887  ctherm = ktf * wa / s
888  ! -- this%xnew available b/c set in parent class (TspAptType) using
889  ! routine set_pointers from the "grandparent" class BndType
890  if (present(rrate)) rrate = ctherm * (this%xnew(igwfnode) - this%xnewpak(n1))
891  if (present(rhsval)) rhsval = dzero
892  if (present(hcofval)) hcofval = ctherm
893  end subroutine sfe_sbcd_term
894 
895  !> @brief Observations
896  !!
897  !! Store the observation type supported by the APT package and override
898  !! BndType%bnd_df_obs
899  !<
900  subroutine sfe_df_obs(this)
901  ! -- modules
902  ! -- dummy
903  class(gwesfetype) :: this
904  ! -- local
905  integer(I4B) :: indx
906  !
907  ! -- Store obs type and assign procedure pointer
908  ! for temperature observation type.
909  call this%obs%StoreObsType('temperature', .false., indx)
910  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
911  !
912  ! -- Store obs type and assign procedure pointer
913  ! for flow between reaches.
914  call this%obs%StoreObsType('flow-ja-face', .true., indx)
915  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid12
916  !
917  ! -- Store obs type and assign procedure pointer
918  ! for from-mvr observation type.
919  call this%obs%StoreObsType('from-mvr', .true., indx)
920  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
921  !
922  ! -- Store obs type and assign procedure pointer
923  ! for to-mvr observation type.
924  call this%obs%StoreObsType('to-mvr', .true., indx)
925  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
926  !
927  ! -- Store obs type and assign procedure pointer
928  ! for storage observation type.
929  call this%obs%StoreObsType('storage', .true., indx)
930  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
931  !
932  ! -- Store obs type and assign procedure pointer
933  ! for constant observation type.
934  call this%obs%StoreObsType('constant', .true., indx)
935  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
936  !
937  ! -- Store obs type and assign procedure pointer
938  ! for observation type: sfe
939  call this%obs%StoreObsType('sfe', .true., indx)
940  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
941  !
942  ! -- Store obs type and assign procedure pointer
943  ! for rainfall observation type.
944  call this%obs%StoreObsType('rainfall', .true., indx)
945  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
946  !
947  ! -- Store obs type and assign procedure pointer
948  ! for evaporation observation type.
949  call this%obs%StoreObsType('evaporation', .true., indx)
950  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
951  !
952  ! -- Store obs type and assign procedure pointer
953  ! for runoff observation type.
954  call this%obs%StoreObsType('runoff', .true., indx)
955  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
956  !
957  ! -- Store obs type and assign procedure pointer
958  ! for inflow observation type.
959  call this%obs%StoreObsType('ext-inflow', .true., indx)
960  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
961  !
962  ! -- Store obs type and assign procedure pointer
963  ! for ext-outflow observation type.
964  call this%obs%StoreObsType('ext-outflow', .true., indx)
965  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
966  !
967  ! -- Store obs type and assign procedure pointer
968  ! for strmbd-cond observation type.
969  call this%obs%StoreObsType('strmbd-cond', .true., indx)
970  this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsid
971  end subroutine sfe_df_obs
972 
973  !> @brief Process package specific obs
974  !!
975  !! Method to process specific observations for this package.
976  !<
977  subroutine sfe_rp_obs(this, obsrv, found)
978  ! -- dummy
979  class(gwesfetype), intent(inout) :: this !< package class
980  type(observetype), intent(inout) :: obsrv !< observation object
981  logical, intent(inout) :: found !< indicate whether observation was found
982  ! -- local
983  !
984  found = .true.
985  select case (obsrv%ObsTypeId)
986  case ('RAINFALL')
987  call this%rp_obs_byfeature(obsrv)
988  case ('EVAPORATION')
989  call this%rp_obs_byfeature(obsrv)
990  case ('RUNOFF')
991  call this%rp_obs_byfeature(obsrv)
992  case ('EXT-INFLOW')
993  call this%rp_obs_byfeature(obsrv)
994  case ('EXT-OUTFLOW')
995  call this%rp_obs_byfeature(obsrv)
996  case ('TO-MVR')
997  call this%rp_obs_byfeature(obsrv)
998  case ('STRMBD-COND')
999  call this%rp_obs_byfeature(obsrv)
1000  case default
1001  found = .false.
1002  end select
1003  end subroutine sfe_rp_obs
1004 
1005  !> @brief Calculate observation value and pass it back to APT
1006  !<
1007  subroutine sfe_bd_obs(this, obstypeid, jj, v, found)
1008  ! -- dummy
1009  class(gwesfetype), intent(inout) :: this
1010  character(len=*), intent(in) :: obstypeid
1011  real(DP), intent(inout) :: v
1012  integer(I4B), intent(in) :: jj
1013  logical, intent(inout) :: found
1014  ! -- local
1015  integer(I4B) :: n1, n2
1016  !
1017  found = .true.
1018  select case (obstypeid)
1019  case ('RAINFALL')
1020  if (this%iboundpak(jj) /= 0) then
1021  call this%sfe_rain_term(jj, n1, n2, v)
1022  end if
1023  case ('EVAPORATION')
1024  if (this%iboundpak(jj) /= 0) then
1025  call this%sfe_evap_term(jj, n1, n2, v)
1026  end if
1027  case ('RUNOFF')
1028  if (this%iboundpak(jj) /= 0) then
1029  call this%sfe_roff_term(jj, n1, n2, v)
1030  end if
1031  case ('EXT-INFLOW')
1032  if (this%iboundpak(jj) /= 0) then
1033  call this%sfe_iflw_term(jj, n1, n2, v)
1034  end if
1035  case ('EXT-OUTFLOW')
1036  if (this%iboundpak(jj) /= 0) then
1037  call this%sfe_outf_term(jj, n1, n2, v)
1038  end if
1039  case ('STRMBD-COND')
1040  if (this%iboundpak(jj) /= 0) then
1041  call this%sfe_sbcd_term(jj, n1, n2, v)
1042  end if
1043  case default
1044  found = .false.
1045  end select
1046  end subroutine sfe_bd_obs
1047 
1048  !> @brief Sets the stress period attributes for keyword use.
1049  !<
1050  subroutine sfe_set_stressperiod(this, itemno, keyword, found)
1051  ! -- modules
1053  ! -- dummy
1054  class(gwesfetype), intent(inout) :: this
1055  integer(I4B), intent(in) :: itemno
1056  character(len=*), intent(in) :: keyword
1057  logical, intent(inout) :: found
1058  ! -- local
1059  character(len=LINELENGTH) :: text
1060  integer(I4B) :: ierr
1061  integer(I4B) :: jj
1062  real(DP), pointer :: bndElem => null()
1063  !
1064  ! RAINFALL <rainfall>
1065  ! EVAPORATION <evaporation>
1066  ! RUNOFF <runoff>
1067  ! INFLOW <inflow>
1068  ! WITHDRAWAL <withdrawal>
1069  !
1070  found = .true.
1071  select case (keyword)
1072  case ('RAINFALL')
1073  ierr = this%apt_check_valid(itemno)
1074  if (ierr /= 0) then
1075  goto 999
1076  end if
1077  call this%parser%GetString(text)
1078  jj = 1
1079  bndelem => this%temprain(itemno)
1080  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
1081  this%packName, 'BND', this%tsManager, &
1082  this%iprpak, 'RAINFALL')
1083  case ('EVAPORATION')
1084  ierr = this%apt_check_valid(itemno)
1085  if (ierr /= 0) then
1086  goto 999
1087  end if
1088  call this%parser%GetString(text)
1089  jj = 1
1090  bndelem => this%tempevap(itemno)
1091  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
1092  this%packName, 'BND', this%tsManager, &
1093  this%iprpak, 'EVAPORATION')
1094  case ('RUNOFF')
1095  ierr = this%apt_check_valid(itemno)
1096  if (ierr /= 0) then
1097  goto 999
1098  end if
1099  call this%parser%GetString(text)
1100  jj = 1
1101  bndelem => this%temproff(itemno)
1102  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
1103  this%packName, 'BND', this%tsManager, &
1104  this%iprpak, 'RUNOFF')
1105  case ('INFLOW')
1106  ierr = this%apt_check_valid(itemno)
1107  if (ierr /= 0) then
1108  goto 999
1109  end if
1110  call this%parser%GetString(text)
1111  jj = 1
1112  bndelem => this%tempiflw(itemno)
1113  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
1114  this%packName, 'BND', this%tsManager, &
1115  this%iprpak, 'INFLOW')
1116  case default
1117  !
1118  ! -- Keyword not recognized so return to caller with found = .false.
1119  found = .false.
1120  end select
1121  !
1122 999 continue
1123  end subroutine sfe_set_stressperiod
1124 
1125  !> @brief Read feature information for this advanced package
1126  !<
1127  subroutine sfe_read_cvs(this)
1128  ! -- modules
1131  ! -- dummy
1132  class(gwesfetype), intent(inout) :: this
1133  ! -- local
1134  character(len=LINELENGTH) :: text
1135  character(len=LENBOUNDNAME) :: bndName, bndNameTemp
1136  character(len=9) :: cno
1137  character(len=50), dimension(:), allocatable :: caux
1138  integer(I4B) :: ierr
1139  logical :: isfound, endOfBlock
1140  integer(I4B) :: n
1141  integer(I4B) :: ii, jj
1142  integer(I4B) :: iaux
1143  integer(I4B) :: itmp
1144  integer(I4B) :: nlak
1145  integer(I4B) :: nconn
1146  integer(I4B), dimension(:), pointer, contiguous :: nboundchk
1147  real(DP), pointer :: bndElem => null()
1148  !
1149  ! -- initialize itmp
1150  itmp = 0
1151  !
1152  ! -- allocate apt data
1153  call mem_allocate(this%strt, this%ncv, 'STRT', this%memoryPath)
1154  call mem_allocate(this%ktf, this%ncv, 'KTF', this%memoryPath)
1155  call mem_allocate(this%rfeatthk, this%ncv, 'RFEATTHK', this%memoryPath)
1156  call mem_allocate(this%lauxvar, this%naux, this%ncv, 'LAUXVAR', &
1157  this%memoryPath)
1158  !
1159  ! -- stream boundary and temperatures
1160  if (this%imatrows == 0) then
1161  call mem_allocate(this%iboundpak, this%ncv, 'IBOUND', this%memoryPath)
1162  call mem_allocate(this%xnewpak, this%ncv, 'XNEWPAK', this%memoryPath)
1163  end if
1164  call mem_allocate(this%xoldpak, this%ncv, 'XOLDPAK', this%memoryPath)
1165  !
1166  ! -- allocate character storage not managed by the memory manager
1167  allocate (this%featname(this%ncv)) ! ditch after boundnames allocated??
1168  !allocate(this%status(this%ncv))
1169  !
1170  do n = 1, this%ncv
1171  this%strt(n) = dep20
1172  this%ktf(n) = dzero
1173  this%rfeatthk(n) = dzero
1174  this%lauxvar(:, n) = dzero
1175  this%xoldpak(n) = dep20
1176  if (this%imatrows == 0) then
1177  this%iboundpak(n) = 1
1178  this%xnewpak(n) = dep20
1179  end if
1180  end do
1181  !
1182  ! -- allocate local storage for aux variables
1183  if (this%naux > 0) then
1184  allocate (caux(this%naux))
1185  end if
1186  !
1187  ! -- allocate and initialize temporary variables
1188  allocate (nboundchk(this%ncv))
1189  do n = 1, this%ncv
1190  nboundchk(n) = 0
1191  end do
1192  !
1193  ! -- get packagedata block
1194  call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
1195  supportopenclose=.true.)
1196  !
1197  ! -- parse locations block if detected
1198  if (isfound) then
1199  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
1200  ' PACKAGEDATA'
1201  nlak = 0
1202  nconn = 0
1203  do
1204  call this%parser%GetNextLine(endofblock)
1205  if (endofblock) exit
1206  n = this%parser%GetInteger()
1207 
1208  if (n < 1 .or. n > this%ncv) then
1209  write (errmsg, '(a,1x,i6)') &
1210  'Itemno must be > 0 and <= ', this%ncv
1211  call store_error(errmsg)
1212  cycle
1213  end if
1214  !
1215  ! -- increment nboundchk
1216  nboundchk(n) = nboundchk(n) + 1
1217  !
1218  ! -- strt
1219  this%strt(n) = this%parser%GetDouble()
1220  !
1221  ! -- read additional thermal conductivity terms
1222  this%ktf(n) = this%parser%GetDouble()
1223  this%rfeatthk(n) = this%parser%GetDouble()
1224  if (this%rfeatthk(n) <= dzero) then
1225  write (errmsg, '(4x,a)') &
1226  '****ERROR. Specified thickness used for thermal &
1227  &conduction MUST BE > 0 else divide by zero error occurs'
1228  call store_error(errmsg)
1229  cycle
1230  end if
1231  !
1232  ! -- get aux data
1233  do iaux = 1, this%naux
1234  call this%parser%GetString(caux(iaux))
1235  end do
1236 
1237  ! -- set default bndName
1238  write (cno, '(i9.9)') n
1239  bndname = 'Feature'//cno
1240 
1241  ! -- featname
1242  if (this%inamedbound /= 0) then
1243  call this%parser%GetStringCaps(bndnametemp)
1244  if (bndnametemp /= '') then
1245  bndname = bndnametemp
1246  end if
1247  end if
1248  this%featname(n) = bndname
1249 
1250  ! -- fill time series aware data
1251  ! -- fill aux data
1252  do jj = 1, this%naux
1253  text = caux(jj)
1254  ii = n
1255  bndelem => this%lauxvar(jj, ii)
1256  call read_value_or_time_series_adv(text, ii, jj, bndelem, &
1257  this%packName, 'AUX', &
1258  this%tsManager, this%iprpak, &
1259  this%auxname(jj))
1260  end do
1261  !
1262  nlak = nlak + 1
1263  end do
1264  !
1265  ! -- check for duplicate or missing lakes
1266  do n = 1, this%ncv
1267  if (nboundchk(n) == 0) then
1268  write (errmsg, '(a,1x,i0)') 'No data specified for feature', n
1269  call store_error(errmsg)
1270  else if (nboundchk(n) > 1) then
1271  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1272  'Data for feature', n, 'specified', nboundchk(n), 'times'
1273  call store_error(errmsg)
1274  end if
1275  end do
1276  !
1277  write (this%iout, '(1x,a)') &
1278  'END OF '//trim(adjustl(this%text))//' PACKAGEDATA'
1279  else
1280  call store_error('Required packagedata block not found.')
1281  end if
1282  !
1283  ! -- terminate if any errors were detected
1284  if (count_errors() > 0) then
1285  call this%parser%StoreErrorUnit()
1286  end if
1287  !
1288  ! -- deallocate local storage for aux variables
1289  if (this%naux > 0) then
1290  deallocate (caux)
1291  end if
1292  !
1293  ! -- deallocate local storage for nboundchk
1294  deallocate (nboundchk)
1295  end subroutine sfe_read_cvs
1296 
1297 end module gwesfemodule
This module contains the base boundary package.
class(bndtype) function, pointer, public getbndfromlist(list, idx)
Get boundary from package list.
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
real(dp), parameter dep20
real constant 1e20
Definition: Constants.f90:91
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 lenbudtxt
maximum length of a budget component names
Definition: Constants.f90:37
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
subroutine sfe_bd_obs(this, obstypeid, jj, v, found)
Calculate observation value and pass it back to APT.
Definition: gwe-sfe.f90:1008
character(len= *), parameter flowtype
Definition: gwe-sfe.f90:57
subroutine sfe_df_obs(this)
Observations.
Definition: gwe-sfe.f90:901
subroutine allocate_scalars(this)
Allocate scalars specific to the streamflow energy transport (SFE) package.
Definition: gwe-sfe.f90:647
subroutine sfe_setup_budobj(this, idx)
Set up the budget object that stores all the sfe flows.
Definition: gwe-sfe.f90:457
integer(i4b) function sfe_get_nbudterms(this)
Function to return the number of budget terms just for this package.
Definition: gwe-sfe.f90:439
subroutine sfe_outf_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Outflow term.
Definition: gwe-sfe.f90:835
subroutine sfe_iflw_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Inflow Term.
Definition: gwe-sfe.f90:808
subroutine sfe_set_stressperiod(this, itemno, keyword, found)
Sets the stress period attributes for keyword use.
Definition: gwe-sfe.f90:1051
character(len=16) text
Definition: gwe-sfe.f90:58
subroutine sfe_rp_obs(this, obsrv, found)
Process package specific obs.
Definition: gwe-sfe.f90:978
subroutine sfe_fill_budobj(this, idx, x, flowja, ccratin, ccratout)
Copy flow terms into thisbudobj.
Definition: gwe-sfe.f90:557
subroutine, public sfe_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, fmi, eqnsclfac, gwecommon, dvt, dvu, dvua)
Create a new sfe package.
Definition: gwe-sfe.f90:108
subroutine sfe_sbcd_term(this, ientry, n1, igwfnode, rrate, rhsval, hcofval)
Streambed conduction term.
Definition: gwe-sfe.f90:862
subroutine sfe_rain_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Rain term.
Definition: gwe-sfe.f90:731
subroutine sfe_evap_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Evaporative term.
Definition: gwe-sfe.f90:755
character(len= *), parameter ftype
Definition: gwe-sfe.f90:56
subroutine sfe_solve(this)
@ brief Add terms specific to sfr to the explicit sfe solve
Definition: gwe-sfe.f90:384
subroutine sfe_allocate_arrays(this)
Allocate arrays specific to the streamflow energy transport (SFE) package.
Definition: gwe-sfe.f90:674
subroutine sfe_roff_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
Runoff term.
Definition: gwe-sfe.f90:780
subroutine sfe_read_cvs(this)
Read feature information for this advanced package.
Definition: gwe-sfe.f90:1128
subroutine sfe_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
Add matrix terms related to SFE.
Definition: gwe-sfe.f90:287
subroutine sfe_da(this)
Deallocate memory.
Definition: gwe-sfe.f90:702
subroutine find_sfe_package(this)
Find corresponding sfe package.
Definition: gwe-sfe.f90:168
This module defines variable data types.
Definition: kind.f90:8
This module contains the derived types ObserveType and ObsDataType.
Definition: Observe.f90:15
This module contains the SFR package methods.
Definition: gwf-sfr.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
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
subroutine, public read_value_or_time_series_adv(textInput, ii, jj, bndElem, pkgName, auxOrBnd, tsManager, iprpak, varName)
Call this subroutine from advanced packages to define timeseries link for a variable (varName).
subroutine, public apt_process_obsid(obsrv, dis, inunitobs, iout)
Process observation IDs for an advanced package.
Definition: tsp-apt.f90:2827
subroutine, public apt_process_obsid12(obsrv, dis, inunitobs, iout)
Process observation IDs for a package.
Definition: tsp-apt.f90:2870
@ brief BndType
Data for sharing among multiple packages. Originally read in from.