MODFLOW 6  version 6.7.0.dev0
USGS Modular Hydrologic Model
All Classes Namespaces Files Functions Variables Enumerator Pages
gwf-ghb.f90
Go to the documentation of this file.
1 module ghbmodule
2  use kindmodule, only: dp, i4b
4  use simvariablesmodule, only: errmsg
7  use bndmodule, only: bndtype
8  use bndextmodule, only: bndexttype
11  !
12  implicit none
13  !
14  private
15  public :: ghb_create
16  public :: ghbtype
17  !
18  character(len=LENFTYPE) :: ftype = 'GHB'
19  character(len=LENPACKAGENAME) :: text = ' GHB'
20  !
21  type, extends(bndexttype) :: ghbtype
22  real(dp), dimension(:), pointer, contiguous :: bhead => null() !< GHB boundary head
23  real(dp), dimension(:), pointer, contiguous :: cond => null() !< GHB hydraulic conductance
24  contains
25  procedure :: allocate_arrays => ghb_allocate_arrays
26  procedure :: source_options => ghb_options
27  procedure :: log_ghb_options
28  procedure :: bnd_rp => ghb_rp
29  procedure :: bnd_ck => ghb_ck
30  procedure :: bnd_cf => ghb_cf
31  procedure :: bnd_fc => ghb_fc
32  procedure :: bnd_da => ghb_da
33  procedure :: define_listlabel
34  procedure :: bound_value => ghb_bound_value
35  procedure :: cond_mult
36  ! -- methods for observations
37  procedure, public :: bnd_obs_supported => ghb_obs_supported
38  procedure, public :: bnd_df_obs => ghb_df_obs
39  procedure, public :: ghb_store_user_cond
40  end type ghbtype
41 
42 contains
43 
44  !> @brief Create a New Ghb Package and point bndobj to the new package
45  !<
46  subroutine ghb_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
47  mempath)
48  ! -- dummy
49  class(bndtype), pointer :: packobj
50  integer(I4B), intent(in) :: id
51  integer(I4B), intent(in) :: ibcnum
52  integer(I4B), intent(in) :: inunit
53  integer(I4B), intent(in) :: iout
54  character(len=*), intent(in) :: namemodel
55  character(len=*), intent(in) :: pakname
56  character(len=*), intent(in) :: mempath
57  ! -- local
58  type(ghbtype), pointer :: ghbobj
59  !
60  ! -- allocate the object and assign values to object variables
61  allocate (ghbobj)
62  packobj => ghbobj
63  !
64  ! -- create name and memory path
65  call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath)
66  packobj%text = text
67  !
68  ! -- allocate scalars
69  call packobj%allocate_scalars()
70  !
71  ! -- initialize package
72  call packobj%pack_initialize()
73  !
74  packobj%inunit = inunit
75  packobj%iout = iout
76  packobj%id = id
77  packobj%ibcnum = ibcnum
78  packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
79  end subroutine ghb_create
80 
81  !> @brief Deallocate memory
82  !<
83  subroutine ghb_da(this)
84  ! -- modules
86  ! -- dummy
87  class(ghbtype) :: this
88  !
89  ! -- Deallocate parent package
90  call this%BndExtType%bnd_da()
91  !
92  ! -- arrays
93  call mem_deallocate(this%bhead, 'BHEAD', this%memoryPath)
94  call mem_deallocate(this%cond, 'COND', this%memoryPath)
95  end subroutine ghb_da
96 
97  !> @brief Set options specific to GhbType
98  !<
99  subroutine ghb_options(this)
100  ! -- modules
104  ! -- dummy
105  class(ghbtype), intent(inout) :: this
106  ! -- local
107  type(gwfghbparamfoundtype) :: found
108  !
109  ! -- source base class options
110  call this%BndExtType%source_options()
111  !
112  ! -- source options from input context
113  call mem_set_value(this%imover, 'MOVER', this%input_mempath, found%mover)
114  !
115  ! -- log ghb specific options
116  call this%log_ghb_options(found)
117  end subroutine ghb_options
118 
119  !> @brief Log options specific to GhbType
120  !<
121  subroutine log_ghb_options(this, found)
122  ! -- modules
124  ! -- dummy
125  class(ghbtype), intent(inout) :: this !< BndExtType object
126  type(gwfghbparamfoundtype), intent(in) :: found
127  !
128  ! -- log found options
129  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text)) &
130  //' OPTIONS'
131  !
132  if (found%mover) then
133  write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
134  end if
135  !
136  ! -- close logging block
137  write (this%iout, '(1x,a)') &
138  'END OF '//trim(adjustl(this%text))//' OPTIONS'
139  end subroutine log_ghb_options
140 
141  !> @brief Allocate arrays
142  !<
143  subroutine ghb_allocate_arrays(this, nodelist, auxvar)
144  ! -- modules
146  ! -- dummy
147  class(ghbtype) :: this
148  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
149  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
150  !
151  ! -- call base type allocate arrays
152  call this%BndExtType%allocate_arrays(nodelist, auxvar)
153  !
154  ! -- set ghb input context pointers
155  call mem_setptr(this%bhead, 'BHEAD', this%input_mempath)
156  call mem_setptr(this%cond, 'COND', this%input_mempath)
157  !
158  ! --checkin ghb input context pointers
159  call mem_checkin(this%bhead, 'BHEAD', this%memoryPath, &
160  'BHEAD', this%input_mempath)
161  call mem_checkin(this%cond, 'COND', this%memoryPath, &
162  'COND', this%input_mempath)
163  end subroutine ghb_allocate_arrays
164 
165  !> @brief Read and prepare
166  !<
167  subroutine ghb_rp(this)
168  ! -- modules
169  use tdismodule, only: kper
170  ! -- dummy
171  class(ghbtype), intent(inout) :: this
172  !
173  if (this%iper /= kper) return
174  !
175  ! -- Call the parent class read and prepare
176  call this%BndExtType%bnd_rp()
177  !
178  ! -- store user cond
179  if (this%ivsc == 1) then
180  call this%ghb_store_user_cond()
181  end if
182  !
183  ! -- Write the list to iout if requested
184  if (this%iprpak /= 0) then
185  call this%write_list()
186  end if
187  end subroutine ghb_rp
188 
189  !> @brief Check ghb boundary condition data
190  !<
191  subroutine ghb_ck(this)
192  ! -- modules
193  use constantsmodule, only: linelength
195  ! -- dummy
196  class(ghbtype), intent(inout) :: this
197  ! -- local
198  character(len=LINELENGTH) :: errmsg
199  integer(I4B) :: i
200  integer(I4B) :: node
201  real(DP) :: bt
202  ! -- formats
203  character(len=*), parameter :: fmtghberr = &
204  "('GHB BOUNDARY (',i0,') HEAD (',f10.3,') IS LESS THAN CELL &
205  &BOTTOM (',f10.3,')')"
206  character(len=*), parameter :: fmtcondmulterr = &
207  "('GHB BOUNDARY (',i0,') CONDUCTANCE MULTIPLIER (',g10.3,') IS &
208  &LESS THAN ZERO')"
209  character(len=*), parameter :: fmtconderr = &
210  "('GHB BOUNDARY (',i0,') CONDUCTANCE (',g10.3,') IS LESS THAN &
211  &ZERO')"
212  !
213  ! -- check stress period data
214  do i = 1, this%nbound
215  node = this%nodelist(i)
216  bt = this%dis%bot(node)
217  ! -- accumulate errors
218  if (this%bhead(i) < bt .and. this%icelltype(node) /= 0) then
219  write (errmsg, fmt=fmtghberr) i, this%bhead(i), bt
220  call store_error(errmsg)
221  end if
222  if (this%iauxmultcol > 0) then
223  if (this%auxvar(this%iauxmultcol, i) < dzero) then
224  write (errmsg, fmt=fmtcondmulterr) &
225  i, this%auxvar(this%iauxmultcol, i)
226  call store_error(errmsg)
227  end if
228  end if
229  if (this%cond(i) < dzero) then
230  write (errmsg, fmt=fmtconderr) i, this%cond(i)
231  call store_error(errmsg)
232  end if
233  end do
234  !
235  !write summary of ghb package error messages
236  if (count_errors() > 0) then
237  call store_error_unit(this%inunit)
238  end if
239  end subroutine ghb_ck
240 
241  !> @brief Formulate the HCOF and RHS terms
242  !!
243  !! Skip if no GHBs
244  !<
245  subroutine ghb_cf(this)
246  ! -- dummy
247  class(ghbtype) :: this
248  ! -- local
249  integer(I4B) :: i, node
250  !
251  ! -- Return if no ghbs
252  if (this%nbound .eq. 0) return
253  !
254  ! -- Calculate hcof and rhs for each ghb entry
255  do i = 1, this%nbound
256  node = this%nodelist(i)
257  if (this%ibound(node) .le. 0) then
258  this%hcof(i) = dzero
259  this%rhs(i) = dzero
260  cycle
261  end if
262  this%hcof(i) = -this%cond_mult(i)
263  this%rhs(i) = -this%cond_mult(i) * this%bhead(i)
264  end do
265  end subroutine ghb_cf
266 
267  !> @brief Copy rhs and hcof into solution rhs and amat
268  !<
269  subroutine ghb_fc(this, rhs, ia, idxglo, matrix_sln)
270  ! -- dummy
271  class(ghbtype) :: this
272  real(DP), dimension(:), intent(inout) :: rhs
273  integer(I4B), dimension(:), intent(in) :: ia
274  integer(I4B), dimension(:), intent(in) :: idxglo
275  class(matrixbasetype), pointer :: matrix_sln
276  ! -- local
277  integer(I4B) :: i, n, ipos
278  real(DP) :: cond, bhead, qghb
279  !
280  ! -- pakmvrobj fc
281  if (this%imover == 1) then
282  call this%pakmvrobj%fc()
283  end if
284  !
285  ! -- Copy package rhs and hcof into solution rhs and amat
286  do i = 1, this%nbound
287  n = this%nodelist(i)
288  rhs(n) = rhs(n) + this%rhs(i)
289  ipos = ia(n)
290  call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))
291  !
292  ! -- If mover is active and this boundary is discharging,
293  ! store available water (as positive value).
294  bhead = this%bhead(i)
295  if (this%imover == 1 .and. this%xnew(n) > bhead) then
296  cond = this%cond_mult(i)
297  qghb = cond * (this%xnew(n) - bhead)
298  call this%pakmvrobj%accumulate_qformvr(i, qghb)
299  end if
300  end do
301  end subroutine ghb_fc
302 
303  !> @brief Define the list heading that is written to iout when PRINT_INPUT
304  !! option is used
305  !<
306  subroutine define_listlabel(this)
307  ! -- dummy
308  class(ghbtype), intent(inout) :: this
309  !
310  ! -- create the header list label
311  this%listlabel = trim(this%filtyp)//' NO.'
312  if (this%dis%ndim == 3) then
313  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
314  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
315  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
316  elseif (this%dis%ndim == 2) then
317  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
318  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
319  else
320  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
321  end if
322  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'STAGE'
323  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'CONDUCTANCE'
324  if (this%inamedbound == 1) then
325  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
326  end if
327  end subroutine define_listlabel
328 
329  ! -- Procedures related to observations
330 
331  !> @brief Return true because GHB package supports observations
332  !!
333  !! Overrides BndType%bnd_obs_supported()
334  !<
335  logical function ghb_obs_supported(this)
336  implicit none
337  ! -- dummy
338  class(ghbtype) :: this
339  !
340  ghb_obs_supported = .true.
341  end function ghb_obs_supported
342 
343  !> @brief Store observation type supported by GHB package
344  !!
345  !! Overrides BndType%bnd_df_obs
346  !<
347  subroutine ghb_df_obs(this)
348  implicit none
349  ! -- dummy
350  class(ghbtype) :: this
351  ! -- local
352  integer(I4B) :: indx
353  !
354  call this%obs%StoreObsType('ghb', .true., indx)
355  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
356  !
357  ! -- Store obs type and assign procedure pointer
358  ! for to-mvr observation type.
359  call this%obs%StoreObsType('to-mvr', .true., indx)
360  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
361  end subroutine ghb_df_obs
362 
363  !> @brief Store user-specified conductance for GHB boundary type
364  !<
365  subroutine ghb_store_user_cond(this)
366  ! -- dummy
367  class(ghbtype), intent(inout) :: this !< BndExtType object
368  ! -- local
369  integer(I4B) :: n
370  !
371  ! -- store backup copy of conductance values
372  do n = 1, this%nbound
373  this%condinput(n) = this%cond_mult(n)
374  end do
375  end subroutine ghb_store_user_cond
376 
377  !> @brief Apply multiplier to GHB conductance if option AUXMULTCOL is used
378  !<
379  function cond_mult(this, row) result(cond)
380  ! -- modules
381  use constantsmodule, only: dzero
382  ! -- dummy variables
383  class(ghbtype), intent(inout) :: this !< BndExtType object
384  integer(I4B), intent(in) :: row
385  ! -- result
386  real(dp) :: cond
387  !
388  if (this%iauxmultcol > 0) then
389  cond = this%cond(row) * this%auxvar(this%iauxmultcol, row)
390  else
391  cond = this%cond(row)
392  end if
393  end function cond_mult
394 
395  !> @brief Return requested boundary value
396  !<
397  function ghb_bound_value(this, col, row) result(bndval)
398  ! -- modules
399  use constantsmodule, only: dzero
400  ! -- dummy
401  class(ghbtype), intent(inout) :: this !< BndExtType object
402  integer(I4B), intent(in) :: col
403  integer(I4B), intent(in) :: row
404  ! -- result
405  real(dp) :: bndval
406  !
407  select case (col)
408  case (1)
409  bndval = this%bhead(row)
410  case (2)
411  bndval = this%cond_mult(row)
412  case default
413  errmsg = 'Programming error. GHB bound value requested column '&
414  &'outside range of ncolbnd (2).'
415  call store_error(errmsg)
416  call store_error_filename(this%input_fname)
417  end select
418  end function ghb_bound_value
419 
420 end module ghbmodule
This module contains the extended boundary package.
This module contains the base boundary package.
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 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
subroutine ghb_cf(this)
Formulate the HCOF and RHS terms.
Definition: gwf-ghb.f90:246
subroutine ghb_allocate_arrays(this, nodelist, auxvar)
Allocate arrays.
Definition: gwf-ghb.f90:144
subroutine log_ghb_options(this, found)
Log options specific to GhbType.
Definition: gwf-ghb.f90:122
subroutine ghb_ck(this)
Check ghb boundary condition data.
Definition: gwf-ghb.f90:192
subroutine ghb_da(this)
Deallocate memory.
Definition: gwf-ghb.f90:84
subroutine ghb_fc(this, rhs, ia, idxglo, matrix_sln)
Copy rhs and hcof into solution rhs and amat.
Definition: gwf-ghb.f90:270
subroutine ghb_df_obs(this)
Store observation type supported by GHB package.
Definition: gwf-ghb.f90:348
character(len=lenpackagename) text
Definition: gwf-ghb.f90:19
real(dp) function ghb_bound_value(this, col, row)
Return requested boundary value.
Definition: gwf-ghb.f90:398
real(dp) function cond_mult(this, row)
Apply multiplier to GHB conductance if option AUXMULTCOL is used.
Definition: gwf-ghb.f90:380
character(len=lenftype) ftype
Definition: gwf-ghb.f90:18
subroutine ghb_store_user_cond(this)
Store user-specified conductance for GHB boundary type.
Definition: gwf-ghb.f90:366
logical function ghb_obs_supported(this)
Return true because GHB package supports observations.
Definition: gwf-ghb.f90:336
subroutine ghb_rp(this)
Read and prepare.
Definition: gwf-ghb.f90:168
subroutine, public ghb_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, mempath)
Create a New Ghb Package and point bndobj to the new package.
Definition: gwf-ghb.f90:48
subroutine ghb_options(this)
Set options specific to GhbType.
Definition: gwf-ghb.f90:100
subroutine define_listlabel(this)
Define the list heading that is written to iout when PRINT_INPUT option is used.
Definition: gwf-ghb.f90:307
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 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
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), pointer, public kper
current stress period number
Definition: tdis.f90:23
@ brief BndType
This class is used to store a single deferred-length character string. It was designed to work in an ...
Definition: CharString.f90:23