MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
gwf-chd.f90
Go to the documentation of this file.
1 module chdmodule
2  !
3  use kindmodule, only: dp, i4b
6  use simvariablesmodule, only: errmsg
10  use bndmodule, only: bndtype
11  use bndextmodule, only: bndexttype
12  use observemodule, only: observetype
16  !
17  implicit none
18  !
19  private
20  public :: chd_create, chdtype
21  !
22  character(len=LENFTYPE) :: ftype = 'CHD'
23  character(len=LENPACKAGENAME) :: text = ' CHD'
24  !
25  type, extends(bndexttype) :: chdtype
26  real(dp), dimension(:), pointer, contiguous :: head => null() !< constant head array
27  real(dp), dimension(:), pointer, contiguous :: ratechdin => null() !simulated flows into constant head (excluding other chds)
28  real(dp), dimension(:), pointer, contiguous :: ratechdout => null() !simulated flows out of constant head (excluding to other chds)
29  contains
30  procedure :: bnd_rp => chd_rp
31  procedure :: bnd_ad => chd_ad
32  procedure :: bnd_ck => chd_ck
33  procedure :: bnd_fc => chd_fc
34  procedure :: bnd_cq => chd_cq
35  procedure :: bnd_bd => chd_bd
36  procedure :: bnd_da => chd_da
37  procedure :: allocate_arrays => chd_allocate_arrays
38  procedure :: define_listlabel
39  procedure :: bound_value => chd_bound_value
40  procedure :: head_mult
41  ! -- methods for observations
42  procedure, public :: bnd_obs_supported => chd_obs_supported
43  procedure, public :: bnd_df_obs => chd_df_obs
44  !
45  procedure, private :: calc_chd_rate
46  end type chdtype
47 
48 contains
49 
50  !> @brief Create a new constant head package
51  !!
52  !! Routine points packobj to the newly created package
53  !<
54  subroutine chd_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
55  mempath)
56  ! -- dummy
57  class(bndtype), pointer :: packobj
58  integer(I4B), intent(in) :: id
59  integer(I4B), intent(in) :: ibcnum
60  integer(I4B), intent(in) :: inunit
61  integer(I4B), intent(in) :: iout
62  character(len=*), intent(in) :: namemodel
63  character(len=*), intent(in) :: pakname
64  character(len=*), intent(in) :: mempath
65  ! -- local
66  type(chdtype), pointer :: chdobj
67  !
68  ! -- allocate the object and assign values to object variables
69  allocate (chdobj)
70  packobj => chdobj
71  !
72  ! -- create name and memory path
73  call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath)
74  packobj%text = text
75  !
76  ! -- allocate scalars
77  call chdobj%allocate_scalars()
78  !
79  ! -- initialize package
80  call packobj%pack_initialize()
81  !
82  ! -- store values
83  packobj%inunit = inunit
84  packobj%iout = iout
85  packobj%id = id
86  packobj%ibcnum = ibcnum
87  packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
88  end subroutine chd_create
89 
90  !> @brief Allocate arrays specific to the constant head package
91  !<
92  subroutine chd_allocate_arrays(this, nodelist, auxvar)
93  ! -- modules
95  ! -- dummy
96  class(chdtype) :: this
97  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
98  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
99  ! -- local
100  integer(I4B) :: i
101  !
102  ! -- call standard BndType allocate scalars
103  call this%BndExtType%allocate_arrays(nodelist, auxvar)
104  !
105  ! -- allocate ratechdex
106  call mem_allocate(this%ratechdin, this%maxbound, 'RATECHDIN', this%memoryPath)
107  call mem_allocate(this%ratechdout, this%maxbound, 'RATECHDOUT', &
108  this%memoryPath)
109  do i = 1, this%maxbound
110  this%ratechdin(i) = dzero
111  this%ratechdout(i) = dzero
112  end do
113  !
114  ! -- set constant head array input context pointer
115  call mem_setptr(this%head, 'HEAD', this%input_mempath)
116  !
117  ! -- checkin constant head array input context pointer
118  call mem_checkin(this%head, 'HEAD', this%memoryPath, &
119  'HEAD', this%input_mempath)
120  end subroutine chd_allocate_arrays
121 
122  !> @brief Constant concentration/temperature read and prepare (rp) routine
123  !<
124  subroutine chd_rp(this)
125  !
126  use tdismodule, only: kper
127  ! -- dummy
128  class(chdtype), intent(inout) :: this
129  ! -- local
130  character(len=30) :: nodestr
131  integer(I4B) :: i, node, ibd, ierr
132  !
133  if (this%iper /= kper) return
134  !
135  ! -- Reset previous CHDs to active cell
136  do i = 1, this%nbound
137  node = this%nodelist(i)
138  this%ibound(node) = this%ibcnum
139  end do
140  !
141  ! -- Call the parent class read and prepare
142  call this%BndExtType%bnd_rp()
143  !
144  ! -- Set ibound to -(ibcnum + 1) for constant head cells
145  ierr = 0
146  do i = 1, this%nbound
147  node = this%nodelist(i)
148  ibd = this%ibound(node)
149  if (ibd < 0) then
150  call this%dis%noder_to_string(node, nodestr)
151  write (errmsg, '(3a)') &
152  'Cell is already a constant head (', trim(adjustl(nodestr)), ').'
153  call store_error(errmsg)
154  ierr = ierr + 1
155  else
156  this%ibound(node) = -this%ibcnum
157  end if
158  end do
159  !
160  ! -- Stop if errors detected
161  if (ierr > 0) then
162  call store_error_filename(this%input_fname)
163  end if
164  !
165  ! -- Write the list to iout if requested
166  if (this%iprpak /= 0) then
167  call this%write_list()
168  end if
169  end subroutine chd_rp
170 
171  !> @brief Constant head package advance routine
172  !!
173  !! Add package connections to matrix
174  !<
175  subroutine chd_ad(this)
176  ! -- modules
177  ! -- dummy
178  class(chdtype) :: this
179  ! -- local
180  integer(I4B) :: i, node
181  real(DP) :: hb
182  ! -- formats
183  !
184  ! -- Process each entry in the specified-head cell list
185  do i = 1, this%nbound
186  node = this%nodelist(i)
187  hb = this%head_mult(i)
188  !
189  this%xnew(node) = hb
190  this%xold(node) = this%xnew(node)
191  end do
192  !
193  ! -- For each observation, push simulated value and corresponding
194  ! simulation time from "current" to "preceding" and reset
195  ! "current" value.
196  call this%obs%obs_ad()
197  end subroutine chd_ad
198 
199  !> @brief Check constant concentration/temperature boundary condition data
200  !<
201  subroutine chd_ck(this)
202  ! -- modules
203  ! -- dummy
204  class(chdtype), intent(inout) :: this
205  ! -- local
206  character(len=30) :: nodestr
207  integer(I4B) :: i
208  integer(I4B) :: node
209  real(DP) :: bt
210  ! -- formats
211  character(len=*), parameter :: fmtchderr = &
212  "('CHD BOUNDARY ',i0,' HEAD (',g0,') IS LESS THAN CELL &
213  &BOTTOM (',g0,')',' FOR CELL ',a)"
214  !
215  ! -- check stress period data
216  do i = 1, this%nbound
217  node = this%nodelist(i)
218  bt = this%dis%bot(node)
219  ! -- accumulate errors
220  if (this%head_mult(i) < bt .and. this%icelltype(node) /= 0) then
221  call this%dis%noder_to_string(node, nodestr)
222  write (errmsg, fmt=fmtchderr) i, this%head_mult(i), bt, trim(nodestr)
223  call store_error(errmsg)
224  end if
225  end do
226  !
227  ! write summary of chd package error messages
228  if (count_errors() > 0) then
229  call store_error_filename(this%input_fname)
230  end if
231  end subroutine chd_ck
232 
233  !> @brief Override bnd_fc and do nothing
234  !!
235  !! For constant head boundary type, the call to bnd_fc
236  !! needs to be overwritten to do nothing
237  !<
238  subroutine chd_fc(this, rhs, ia, idxglo, matrix_sln)
239  ! -- dummy
240  class(chdtype) :: this
241  real(DP), dimension(:), intent(inout) :: rhs
242  integer(I4B), dimension(:), intent(in) :: ia
243  integer(I4B), dimension(:), intent(in) :: idxglo
244  class(matrixbasetype), pointer :: matrix_sln
245  ! -- local
246  end subroutine chd_fc
247 
248  !> @brief Calculate flow associated with constant head boundary
249  !!
250  !! This method overrides bnd_cq()
251  !<
252  subroutine chd_cq(this, x, flowja, iadv)
253  class(chdtype), intent(inout) :: this
254  real(DP), dimension(:), intent(in) :: x
255  real(DP), dimension(:), contiguous, intent(inout) :: flowja
256  integer(I4B), optional, intent(in) :: iadv
257 
258  ! NB: the rate calculation cannot be done until chd_bd below
259 
260  end subroutine chd_cq
261 
262  !> @brief Calculate the CHD cell rates, to be called
263  !< after all updates to the model flowja are done
264  subroutine calc_chd_rate(this)
265  ! -- modules
266  ! -- dummy
267  class(chdtype), intent(inout) :: this
268  ! -- local
269  integer(I4B) :: i
270  integer(I4B) :: ipos
271  integer(I4B) :: node
272  integer(I4B) :: n2
273  integer(I4B) :: idiag
274  real(DP) :: rate
275  real(DP) :: ratein, rateout
276  real(DP) :: q
277  !
278  ! -- If no boundaries, skip flow calculations.
279  if (this%nbound > 0) then
280  !
281  ! -- Loop through each boundary calculating flow.
282  do i = 1, this%nbound
283  node = this%nodelist(i)
284  idiag = this%dis%con%ia(node)
285  rate = dzero
286  ratein = dzero
287  rateout = dzero
288  !
289  ! -- Calculate the flow rate into the cell.
290  do ipos = this%dis%con%ia(node) + 1, &
291  this%dis%con%ia(node + 1) - 1
292  q = this%flowja(ipos)
293  rate = rate - q
294  ! -- only accumulate chin and chout for active
295  ! connected cells
296  n2 = this%dis%con%ja(ipos)
297  if (this%ibound(n2) > 0) then
298  if (q < dzero) then
299  ratein = ratein - q
300  else
301  rateout = rateout + q
302  end if
303  end if
304  end do
305  !
306  ! -- For chd, store total flow in rhs so it is available for other
307  ! calculations
308  this%rhs(i) = -rate
309  this%hcof(i) = dzero
310  !
311  ! -- Save simulated value to simvals array.
312  this%simvals(i) = rate
313  this%ratechdin(i) = ratein
314  this%ratechdout(i) = rateout
315  this%flowja(idiag) = this%flowja(idiag) + rate
316  !
317  end do
318  !
319  end if
320  end subroutine calc_chd_rate
321 
322  !> @brief Add package ratin/ratout to model budget
323  !<
324  subroutine chd_bd(this, model_budget)
325  ! -- add package ratin/ratout to model budget
326  use tdismodule, only: delt
328  class(chdtype) :: this
329  type(budgettype), intent(inout) :: model_budget
330  real(DP) :: ratin
331  real(DP) :: ratout
332  real(DP) :: dum
333  integer(I4B) :: isuppress_output
334 
335  ! For CHDs at an exchange, under some conditions
336  ! (XT3D), the model flowja into the cell is not
337  ! finalized until after exg_cq. So we calculate
338  ! the CHD rate here
339  call this%calc_chd_rate()
340 
341  isuppress_output = 0
342  call rate_accumulator(this%ratechdin(1:this%nbound), ratin, dum)
343  call rate_accumulator(this%ratechdout(1:this%nbound), ratout, dum)
344  call model_budget%addentry(ratin, ratout, delt, this%text, &
345  isuppress_output, this%packName)
346  end subroutine chd_bd
347 
348  !> @brief Deallocate memory
349  !<
350  subroutine chd_da(this)
351  ! -- modules
353  ! -- dummy
354  class(chdtype) :: this
355  !
356  ! -- Deallocate parent package
357  call this%BndExtType%bnd_da()
358  !
359  ! -- arrays
360  call mem_deallocate(this%ratechdin)
361  call mem_deallocate(this%ratechdout)
362  call mem_deallocate(this%head, 'HEAD', this%memoryPath)
363  end subroutine chd_da
364 
365  !> @brief Define the list heading that is written to iout when PRINT_INPUT
366  !! option is used.
367  !<
368  subroutine define_listlabel(this)
369  class(chdtype), intent(inout) :: this
370  !
371  ! -- create the header list label
372  this%listlabel = trim(this%filtyp)//' NO.'
373  if (this%dis%ndim == 3) then
374  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
375  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
376  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
377  elseif (this%dis%ndim == 2) then
378  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
379  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
380  else
381  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
382  end if
383  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'HEAD'
384  if (this%inamedbound == 1) then
385  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
386  end if
387  end subroutine define_listlabel
388 
389  ! -- Procedures related to observations
390 
391  !> @brief Overrides bnd_obs_supported from bndType class
392  !!
393  !! Return true since CHD package supports observations
394  !<
395  logical function chd_obs_supported(this)
396  implicit none
397  !
398  class(chdtype) :: this
399  !
400  chd_obs_supported = .true.
401  end function chd_obs_supported
402 
403  !> @brief Overrides bnd_df_obs from bndType class
404  !!
405  !! (1) Store observation type supported by CHD package and (2) override
406  !! BndType%bnd_df_obs
407  !<
408  subroutine chd_df_obs(this)
409  implicit none
410  ! -- dummy
411  class(chdtype) :: this
412  ! -- local
413  integer(I4B) :: indx
414  !
415  call this%obs%StoreObsType('chd', .true., indx)
416  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
417  end subroutine chd_df_obs
418 
419  !> @brief Apply auxiliary multiplier to specified head if appropriate
420  !<
421  function head_mult(this, row) result(head)
422  ! -- modules
423  use constantsmodule, only: dzero
424  ! -- dummy variables
425  class(chdtype), intent(inout) :: this !< BndExtType object
426  integer(I4B), intent(in) :: row
427  ! -- result
428  real(dp) :: head
429  !
430  if (this%iauxmultcol > 0) then
431  head = this%head(row) * this%auxvar(this%iauxmultcol, row)
432  else
433  head = this%head(row)
434  end if
435  end function head_mult
436 
437  !> @ brief Return a bound value
438  !!
439  !! Return a bound value associated with an ncolbnd index
440  !! and row.
441  !<
442  function chd_bound_value(this, col, row) result(bndval)
443  ! -- modules
444  use constantsmodule, only: dzero
445  ! -- dummy variables
446  class(chdtype), intent(inout) :: this !< BndType object
447  integer(I4B), intent(in) :: col
448  integer(I4B), intent(in) :: row
449  ! -- result
450  real(dp) :: bndval
451  !
452  select case (col)
453  case (1)
454  bndval = this%head_mult(row)
455  case default
456  errmsg = 'Programming error. CHD bound value requested column '&
457  &'outside range of ncolbnd (1).'
458  call store_error(errmsg)
459  call store_error_filename(this%input_fname)
460  end select
461  end function chd_bound_value
462 
463 end module chdmodule
This module contains the extended boundary package.
This module contains the base boundary package.
This module contains the BudgetModule.
Definition: Budget.f90:20
subroutine, public rate_accumulator(flow, rin, rout)
@ brief Rate accumulator subroutine
Definition: Budget.f90:632
subroutine chd_ck(this)
Check constant concentration/temperature boundary condition data.
Definition: gwf-chd.f90:202
character(len=lenpackagename) text
Definition: gwf-chd.f90:23
subroutine calc_chd_rate(this)
Calculate the CHD cell rates, to be called.
Definition: gwf-chd.f90:265
subroutine chd_da(this)
Deallocate memory.
Definition: gwf-chd.f90:351
subroutine define_listlabel(this)
Define the list heading that is written to iout when PRINT_INPUT option is used.
Definition: gwf-chd.f90:369
subroutine chd_cq(this, x, flowja, iadv)
Calculate flow associated with constant head boundary.
Definition: gwf-chd.f90:253
character(len=lenftype) ftype
Definition: gwf-chd.f90:22
real(dp) function chd_bound_value(this, col, row)
@ brief Return a bound value
Definition: gwf-chd.f90:443
subroutine chd_df_obs(this)
Overrides bnd_df_obs from bndType class.
Definition: gwf-chd.f90:409
subroutine chd_rp(this)
Constant concentration/temperature read and prepare (rp) routine.
Definition: gwf-chd.f90:125
subroutine chd_ad(this)
Constant head package advance routine.
Definition: gwf-chd.f90:176
subroutine, public chd_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, mempath)
Create a new constant head package.
Definition: gwf-chd.f90:56
real(dp) function head_mult(this, row)
Apply auxiliary multiplier to specified head if appropriate.
Definition: gwf-chd.f90:422
subroutine chd_allocate_arrays(this, nodelist, auxvar)
Allocate arrays specific to the constant head package.
Definition: gwf-chd.f90:93
logical function chd_obs_supported(this)
Overrides bnd_obs_supported from bndType class.
Definition: gwf-chd.f90:396
subroutine chd_bd(this, model_budget)
Add package ratin/ratout to model budget.
Definition: gwf-chd.f90:325
subroutine chd_fc(this, rhs, ia, idxglo, matrix_sln)
Override bnd_fc and do nothing.
Definition: gwf-chd.f90:239
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
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
integer(i4b), parameter lenftype
maximum length of a package type (DIS, WEL, OC, etc.)
Definition: Constants.f90:39
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
This module defines variable data types.
Definition: kind.f90:8
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
This module contains the derived types ObserveType and ObsDataType.
Definition: Observe.f90:15
This module contains the derived type ObsType.
Definition: Obs.f90:127
subroutine, public defaultobsidprocessor(obsrv, dis, inunitobs, iout)
@ brief Process IDstring provided for each observation
Definition: Obs.f90:246
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
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
type(timeserieslinktype) function, pointer, public gettimeserieslinkfromlist(list, indx)
Get time series link from a list.
@ brief BndType
Derived type for the Budget object.
Definition: Budget.f90:39