MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
DiscretizationBase.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b, lgp
6  dis, disv, disu, &
7  dis2d, disv2d, disu2d, &
9 
13  use simvariablesmodule, only: errmsg
14  use simmodule, only: count_errors, store_error, &
20  use tdismodule, only: kstp, kper, pertim, totim, delt
23 
24  implicit none
25 
26  private
27  public :: disbasetype
28  public :: dis_transform_xy
29 
30  type :: disbasetype
31  character(len=LENMEMPATH) :: memorypath !< path for memory allocation
32  character(len=LENMEMPATH) :: input_mempath = '' !< input context mempath
33  character(len=LENMODELNAME), pointer :: name_model => null() !< name of the model
34  character(len=LINELENGTH), pointer :: input_fname => null() !< input file name
35  character(len=LINELENGTH), pointer :: output_fname => null() !< output file name
36  integer(I4B), pointer :: inunit => null() !< unit number for input file
37  integer(I4B), pointer :: iout => null() !< unit number for output file
38  integer(I4B), pointer :: nodes => null() !< number of nodes in solution
39  integer(I4B), pointer :: nodesuser => null() !< number of user nodes (same as nodes for disu grid)
40  integer(I4B), pointer :: nja => null() !< number of connections plus number of nodes
41  integer(I4B), pointer :: njas => null() !< (nja-nodes)/2
42  integer(I4B), pointer :: lenuni => null() !< length unit
43  integer(I4B), pointer :: ndim => null() !< number of spatial model dimensions (1 for disu grid)
44  integer(I4B), pointer :: icondir => null() !< flag indicating if grid has enough info to calculate connection vectors
45  integer(I4B), pointer :: nogrb => null() !< don't write binary grid file
46  real(dp), dimension(:), pointer, contiguous :: xc => null() !< x-coordinate of the cell center
47  real(dp), dimension(:), pointer, contiguous :: yc => null() !< y-coordinate of the cell center
48  real(dp), pointer :: yorigin => null() !< y-position of the lower-left grid corner (default is 0.)
49  real(dp), pointer :: xorigin => null() !< x-position of the lower-left grid corner (default is 0.)
50  real(dp), pointer :: angrot => null() !< counter-clockwise rotation angle of the lower-left corner (default is 0.0)
51  integer(I4B), dimension(:), pointer, contiguous :: mshape => null() !< shape of the model; (nodes) for DisBaseType
52  real(dp), dimension(:), pointer, contiguous :: top => null() !< (size:nodes) cell top elevation
53  real(dp), dimension(:), pointer, contiguous :: bot => null() !< (size:nodes) cell bottom elevation
54  real(dp), dimension(:), pointer, contiguous :: area => null() !< (size:nodes) cell area, in plan view
55  type(connectionstype), pointer :: con => null() !< connections object
56  type(blockparsertype) :: parser !< object to read blocks
57  real(dp), dimension(:), pointer, contiguous :: dbuff => null() !< helper double array of size nodesuser
58  integer(I4B), dimension(:), pointer, contiguous :: ibuff => null() !< helper int array of size nodesuser
59  integer(I4B), dimension(:), pointer, contiguous :: nodereduced => null() !< (size:nodesuser)contains reduced nodenumber (size 0 if not reduced); -1 means vertical pass through, 0 is idomain = 0
60  integer(I4B), dimension(:), pointer, contiguous :: nodeuser => null() !< (size:nodes) given a reduced nodenumber, provide the user nodenumber (size 0 if not reduced)
61  contains
62  procedure :: dis_df
63  procedure :: dis_ac
64  procedure :: dis_mc
65  procedure :: dis_ar
66  procedure :: dis_da
67  ! -- helper functions
68  !
69  ! -- get_nodenumber is an overloaded integer function that will always
70  ! return the reduced nodenumber. For all grids, get_nodenumber can
71  ! be passed the user nodenumber. For some other grids, it can also
72  ! be passed an index. For dis3d the index is k, i, j, and for
73  ! disv the index is k, n.
74  generic :: get_nodenumber => get_nodenumber_idx1, &
77  procedure :: get_nodenumber_idx1
78  procedure :: get_nodenumber_idx2
79  procedure :: get_nodenumber_idx3
80  procedure :: get_nodeuser
81  procedure :: nodeu_to_string
82  procedure :: nodeu_to_array
83  procedure :: nodeu_from_string
84  procedure :: nodeu_from_cellid
85  procedure :: noder_from_string
86  procedure :: noder_from_cellid
87  procedure :: connection_normal
88  procedure :: connection_vector
89  procedure :: get_dis_type
90  procedure :: get_dis_enum
91  procedure :: supports_layers
92  procedure :: allocate_scalars
93  procedure :: allocate_arrays
94  procedure :: get_ncpl
95  procedure :: get_cell_volume
96  procedure :: get_polyverts
97  procedure :: get_npolyverts
98  procedure :: get_max_npolyverts
99  procedure :: write_grb
100  !
101  procedure :: read_int_array
102  procedure :: read_dbl_array
103  generic, public :: read_grid_array => read_int_array, read_dbl_array
104  procedure, public :: read_layer_array
105  procedure :: fill_int_array
106  procedure :: fill_dbl_array
107  generic, public :: fill_grid_array => fill_int_array, fill_dbl_array
108  procedure, public :: read_list
109  !
110  procedure, public :: record_array
111  procedure, public :: record_connection_array
112  procedure, public :: noder_to_string
113  procedure, public :: noder_to_array
114  procedure, public :: record_srcdst_list_header
115  procedure, private :: record_srcdst_list_entry
116  generic, public :: record_mf6_list_entry => record_srcdst_list_entry
117  procedure, public :: nlarray_to_nodelist
118  procedure, public :: highest_active
119  procedure, public :: highest_saturated
120  procedure, public :: get_area
121  procedure, public :: get_area_factor
122  procedure, public :: get_flow_width
123  procedure, public :: is_3d
124  procedure, public :: is_2d
125  procedure, public :: is_1d
126  end type disbasetype
127 
128 contains
129 
130  !> @brief Define the discretization
131  subroutine dis_df(this)
132  class(disbasetype) :: this
133  call store_error('Programmer error: dis_df must be overridden', &
134  terminate=.true.)
135  end subroutine dis_df
136 
137  !> @brief Add connections to sparse cell connectivity matrix
138  subroutine dis_ac(this, moffset, sparse)
139  ! -- modules
140  use sparsemodule, only: sparsematrix
141  ! -- dummy
142  class(disbasetype) :: this
143  integer(I4B), intent(in) :: moffset
144  type(sparsematrix), intent(inout) :: sparse
145  ! -- local
146  integer(I4B) :: i, j, ipos, iglo, jglo
147  !
148  do i = 1, this%nodes
149  do ipos = this%con%ia(i), this%con%ia(i + 1) - 1
150  j = this%con%ja(ipos)
151  iglo = i + moffset
152  jglo = j + moffset
153  call sparse%addconnection(iglo, jglo, 1)
154  end do
155  end do
156  end subroutine dis_ac
157 
158  !> @brief Map cell connections in the numerical solution coefficient matrix.
159  subroutine dis_mc(this, moffset, idxglo, matrix_sln)
160  ! -- dummy
161  class(disbasetype) :: this
162  integer(I4B), intent(in) :: moffset
163  integer(I4B), dimension(:), intent(inout) :: idxglo
164  class(matrixbasetype), pointer :: matrix_sln
165  ! -- local
166  integer(I4B) :: i, j, ipos, iglo, jglo
167  !
168  do i = 1, this%nodes
169  iglo = i + moffset
170  do ipos = this%con%ia(i), this%con%ia(i + 1) - 1
171  j = this%con%ja(ipos)
172  jglo = j + moffset
173  idxglo(ipos) = matrix_sln%get_position(iglo, jglo)
174  end do
175  end do
176  end subroutine dis_mc
177 
178  !> @brief Allocate and setup variables, and write binary grid file.
179  subroutine dis_ar(this, icelltype)
180  ! -- dummy
181  class(disbasetype) :: this
182  integer(I4B), dimension(:), intent(in) :: icelltype
183  ! -- local
184  integer(I4B), dimension(:), allocatable :: ict
185  integer(I4B) :: nu, nr
186  !
187  ! -- Expand icelltype to full grid; fill with 0 if cell is excluded
188  allocate (ict(this%nodesuser))
189  do nu = 1, this%nodesuser
190  nr = this%get_nodenumber(nu, 0)
191  if (nr > 0) then
192  ict(nu) = icelltype(nr)
193  else
194  ict(nu) = 0
195  end if
196  end do
197  !
198  if (this%nogrb == 0) call this%write_grb(ict)
199  end subroutine dis_ar
200 
201  !> @brief Write a binary grid file
202  subroutine write_grb(this, icelltype)
203  class(disbasetype) :: this
204  integer(I4B), dimension(:), intent(in) :: icelltype
205  call store_error('Programmer error: write_grb must be overridden', &
206  terminate=.true.)
207  end subroutine write_grb
208 
209  !> @brier Deallocate variables
210  subroutine dis_da(this)
211  ! -- modules
213  ! -- dummy
214  class(disbasetype) :: this
215  !
216  ! -- Strings
217  deallocate (this%name_model)
218  deallocate (this%input_fname)
219  deallocate (this%output_fname)
220  !
221  ! -- Scalars
222  call mem_deallocate(this%inunit)
223  call mem_deallocate(this%iout)
224  call mem_deallocate(this%nodes)
225  call mem_deallocate(this%nodesuser)
226  call mem_deallocate(this%ndim)
227  call mem_deallocate(this%icondir)
228  call mem_deallocate(this%nogrb)
229  call mem_deallocate(this%xorigin)
230  call mem_deallocate(this%yorigin)
231  call mem_deallocate(this%angrot)
232  call mem_deallocate(this%nja)
233  call mem_deallocate(this%njas)
234  call mem_deallocate(this%lenuni)
235  !
236  ! -- Arrays
237  call mem_deallocate(this%mshape)
238  call mem_deallocate(this%xc)
239  call mem_deallocate(this%yc)
240  call mem_deallocate(this%top)
241  call mem_deallocate(this%bot)
242  call mem_deallocate(this%area)
243  call mem_deallocate(this%dbuff)
244  call mem_deallocate(this%ibuff)
245  !
246  ! -- Connections
247  call this%con%con_da()
248  deallocate (this%con)
249  end subroutine dis_da
250 
251  !> @brief Convert a user nodenumber to a string (nodenumber), (k,j), or (k,i,j)
252  subroutine nodeu_to_string(this, nodeu, str)
253  class(disbasetype) :: this
254  integer(I4B), intent(in) :: nodeu
255  character(len=*), intent(inout) :: str
256 
257  call store_error('Programmer error: nodeu_to_string must be overridden', &
258  terminate=.true.)
259  end subroutine nodeu_to_string
260 
261  !> @brief Convert a user nodenumber to an array (nodenumber), (k,j), or (k,i,j)
262  subroutine nodeu_to_array(this, nodeu, arr)
263  class(disbasetype) :: this
264  integer(I4B), intent(in) :: nodeu
265  integer(I4B), dimension(:), intent(inout) :: arr
266 
267  call store_error('Programmer error: nodeu_to_array must be overridden', &
268  terminate=.true.)
269  end subroutine nodeu_to_array
270 
271  !> @brief Convert a reduced nodenumber to a user node number
272  function get_nodeuser(this, noder) result(nodenumber)
273  class(disbasetype) :: this
274  integer(I4B), intent(in) :: noder
275  integer(I4B) :: nodenumber
276 
277  if (this%nodes < this%nodesuser) then
278  nodenumber = this%nodeuser(noder)
279  else
280  nodenumber = noder
281  end if
282  end function get_nodeuser
283 
284  function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
285  class(disbasetype), intent(in) :: this
286  integer(I4B), intent(in) :: nodeu
287  integer(I4B), intent(in) :: icheck
288  integer(I4B) :: nodenumber
289 
290  nodenumber = 0
291  call store_error('Programmer error: get_nodenumber_idx1 must be overridden', &
292  terminate=.true.)
293  end function get_nodenumber_idx1
294 
295  function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber)
296  class(disbasetype), intent(in) :: this
297  integer(I4B), intent(in) :: k, j
298  integer(I4B), intent(in) :: icheck
299  integer(I4B) :: nodenumber
300 
301  nodenumber = 0
302  call store_error('Programmer error: get_nodenumber_idx2 must be overridden', &
303  terminate=.true.)
304  end function get_nodenumber_idx2
305 
306  function get_nodenumber_idx3(this, k, i, j, icheck) result(nodenumber)
307  class(disbasetype), intent(in) :: this
308  integer(I4B), intent(in) :: k, i, j
309  integer(I4B), intent(in) :: icheck
310  integer(I4B) :: nodenumber
311 
312  nodenumber = 0
313  call store_error('Programmer error: get_nodenumber_idx3 must be overridden', &
314  terminate=.true.)
315  end function get_nodenumber_idx3
316 
317  !> @brief Get normal vector components between the cell and a given neighbor.
318  !! The normal points outward from the shared face between noden and nodem.
319  subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, &
320  ipos)
321  class(disbasetype) :: this
322  integer(I4B), intent(in) :: noden !< cell (reduced nn)
323  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
324  integer(I4B), intent(in) :: ihc !< horizontal connection flag
325  real(DP), intent(inout) :: xcomp
326  real(DP), intent(inout) :: ycomp
327  real(DP), intent(inout) :: zcomp
328  integer(I4B), intent(in) :: ipos
329 
330  call store_error('Programmer error: connection_normal must be overridden', &
331  terminate=.true.)
332  end subroutine connection_normal
333 
334  !> @brief Get unit vector components between the cell and a given neighbor.
335  !! Saturation must be provided to compute cell center vertical coordinates.
336  !! Also return the straight-line connection length.
337  subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, &
338  xcomp, ycomp, zcomp, conlen)
339  class(disbasetype) :: this
340  integer(I4B), intent(in) :: noden !< cell (reduced nn)
341  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
342  logical, intent(in) :: nozee
343  real(DP), intent(in) :: satn
344  real(DP), intent(in) :: satm
345  integer(I4B), intent(in) :: ihc !< horizontal connection flag
346  real(DP), intent(inout) :: xcomp
347  real(DP), intent(inout) :: ycomp
348  real(DP), intent(inout) :: zcomp
349  real(DP), intent(inout) :: conlen
350 
351  call store_error('Programmer error: connection_vector must be overridden', &
352  terminate=.true.)
353  end subroutine connection_vector
354 
355  !> @brief Get global (x, y) coordinates from cell-local coordinates.
356  subroutine dis_transform_xy(x, y, xorigin, yorigin, angrot, xglo, yglo)
357  real(dp), intent(in) :: x !< the cell-x coordinate to transform
358  real(dp), intent(in) :: y !< the cell-y coordinate to transform
359  real(dp), intent(in) :: xorigin !< the cell-y coordinate to transform
360  real(dp), intent(in) :: yorigin !< the cell-y coordinate to transform
361  real(dp), intent(in) :: angrot !< the cell-y coordinate to transform
362  real(dp), intent(out) :: xglo !< the global cell-x coordinate
363  real(dp), intent(out) :: yglo !< the global cell-y coordinate
364  ! local
365  real(dp) :: ang
366 
367  xglo = x
368  yglo = y
369 
370  ! first _rotate_ to 'real world'
371  ang = angrot * dpio180
372  if (ang /= dzero) then
373  xglo = x * cos(ang) - y * sin(ang)
374  yglo = x * sin(ang) + y * cos(ang)
375  end if
376 
377  ! then _translate_
378  xglo = xglo + xorigin
379  yglo = yglo + yorigin
380  end subroutine dis_transform_xy
381 
382  !> @brief Get the discretization type (DIS, DISV, or DISU)
383  subroutine get_dis_type(this, dis_type)
384  class(disbasetype), intent(in) :: this
385  character(len=*), intent(out) :: dis_type
386 
387  dis_type = "Not implemented"
388  call store_error('Programmer error: get_dis_type must be overridden', &
389  terminate=.true.)
390  end subroutine get_dis_type
391 
392  !> @brief Get the discretization type enumeration
393  function get_dis_enum(this) result(dis_enum)
394  use constantsmodule, only: disundef
395  class(disbasetype), intent(in) :: this
396  integer(I4B) :: dis_enum
397 
398  dis_enum = disundef
399  call store_error('Programmer error: get_dis_enum must be overridden', &
400  terminate=.true.)
401  end function get_dis_enum
402 
403  !> @brief Allocate and initialize scalar variables
404  subroutine allocate_scalars(this, name_model, input_mempath)
405  ! -- dummy
406  class(disbasetype) :: this
407  character(len=*), intent(in) :: name_model
408  character(len=*), intent(in) :: input_mempath
409  logical(LGP) :: found
410  !
411  ! -- Create memory path
412  this%memoryPath = create_mem_path(name_model, 'DIS')
413  !
414  ! -- Allocate
415  allocate (this%name_model)
416  allocate (this%input_fname)
417  allocate (this%output_fname)
418  !
419  call mem_allocate(this%inunit, 'INUNIT', this%memoryPath)
420  call mem_allocate(this%iout, 'IOUT', this%memoryPath)
421  call mem_allocate(this%nodes, 'NODES', this%memoryPath)
422  call mem_allocate(this%nodesuser, 'NODESUSER', this%memoryPath)
423  call mem_allocate(this%ndim, 'NDIM', this%memoryPath)
424  call mem_allocate(this%icondir, 'ICONDIR', this%memoryPath)
425  call mem_allocate(this%nogrb, 'NOGRB', this%memoryPath)
426  call mem_allocate(this%xorigin, 'XORIGIN', this%memoryPath)
427  call mem_allocate(this%yorigin, 'YORIGIN', this%memoryPath)
428  call mem_allocate(this%angrot, 'ANGROT', this%memoryPath)
429  call mem_allocate(this%nja, 'NJA', this%memoryPath)
430  call mem_allocate(this%njas, 'NJAS', this%memoryPath)
431  call mem_allocate(this%lenuni, 'LENUNI', this%memoryPath)
432  !
433  ! -- Initialize
434  this%name_model = name_model
435  this%input_mempath = input_mempath
436  this%input_fname = ''
437  this%output_fname = ''
438  this%inunit = 0
439  this%iout = 0
440  this%nodes = 0
441  this%nodesuser = 0
442  this%ndim = 1
443  this%icondir = 1
444  this%nogrb = 0
445  this%xorigin = dzero
446  this%yorigin = dzero
447  this%angrot = dzero
448  this%nja = 0
449  this%njas = 0
450  this%lenuni = 0
451  !
452  ! -- update input and output filenames
453  call mem_set_value(this%input_fname, 'INPUT_FNAME', &
454  this%input_mempath, found)
455  call mem_set_value(this%output_fname, 'GRB6_FILENAME', &
456  this%input_mempath, found)
457  if (.not. found) then
458  this%output_fname = trim(this%input_fname)//'.grb'
459  end if
460  end subroutine allocate_scalars
461 
462  !> @brief Allocate and initialize arrays
463  subroutine allocate_arrays(this)
464  class(disbasetype) :: this
465  integer :: isize
466  !
467  ! -- Allocate
468  call mem_allocate(this%mshape, this%ndim, 'MSHAPE', this%memoryPath)
469  call mem_allocate(this%xc, this%nodes, 'XC', this%memoryPath)
470  call mem_allocate(this%yc, this%nodes, 'YC', this%memoryPath)
471  call mem_allocate(this%top, this%nodes, 'TOP', this%memoryPath)
472  call mem_allocate(this%bot, this%nodes, 'BOT', this%memoryPath)
473  call mem_allocate(this%area, this%nodes, 'AREA', this%memoryPath)
474  !
475  ! -- Initialize
476  this%mshape(1) = this%nodes
477  !
478  ! -- Determine size of buff memory
479  if (this%nodes < this%nodesuser) then
480  isize = this%nodesuser
481  else
482  isize = this%nodes
483  end if
484  !
485  ! -- Allocate the arrays
486  call mem_allocate(this%dbuff, isize, 'DBUFF', this%name_model)
487  call mem_allocate(this%ibuff, isize, 'IBUFF', this%name_model)
488  end subroutine allocate_arrays
489 
490  !> @brief Convert a string to a user nodenumber.
491  !!
492  !! If DIS or DISV, read indices. If DISU, read user node number directly.
493  !! If flag_string is present and true, the first token may be
494  !! non-numeric (e.g. boundary name). In this case, return -2.
495  !<
496  function nodeu_from_string(this, lloc, istart, istop, in, iout, line, &
497  flag_string, allow_zero) result(nodeu)
498  ! -- dummy
499  class(disbasetype) :: this
500  integer(I4B), intent(inout) :: lloc
501  integer(I4B), intent(inout) :: istart
502  integer(I4B), intent(inout) :: istop
503  integer(I4B), intent(in) :: in
504  integer(I4B), intent(in) :: iout
505  character(len=*), intent(inout) :: line
506  logical, optional, intent(in) :: flag_string
507  logical, optional, intent(in) :: allow_zero
508  integer(I4B) :: nodeu
509 
510  nodeu = 0
511  call store_error('Programmer error: nodeu_from_string must be overridden', &
512  terminate=.true.)
513  end function nodeu_from_string
514 
515  !> @brief Convert a cellid string to a user nodenumber.
516  !!
517  !! If flag_string is present and true, the first token may be
518  !! non-numeric (e.g. boundary name). In this case, return -2.
519  !!
520  !! If allow_zero is present and true, and all indices are zero, the
521  !! result can be zero. If allow_zero is false, a zero in any index is an error.
522  !<
523  function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, &
524  allow_zero) result(nodeu)
525  ! -- dummy
526  class(disbasetype) :: this
527  character(len=*), intent(inout) :: cellid
528  integer(I4B), intent(in) :: inunit
529  integer(I4B), intent(in) :: iout
530  logical, optional, intent(in) :: flag_string
531  logical, optional, intent(in) :: allow_zero
532  integer(I4B) :: nodeu
533 
534  nodeu = 0
535  call store_error('Programmer error: nodeu_from_cellid must be overridden', &
536  terminate=.true.)
537  end function nodeu_from_cellid
538 
539  !> @brief Convert a string to a reduced nodenumber.
540  !!
541  !! If the model is unstructured; just read user nodenumber.
542  !! If flag_string argument is present and true, the first token in string
543  !! is allowed to be a string (e.g. boundary name). In this case, if a string
544  !! is encountered, return value as -2.
545  !<
546  function noder_from_string(this, lloc, istart, istop, in, iout, line, &
547  flag_string) result(noder)
548  ! -- dummy
549  class(disbasetype) :: this
550  integer(I4B), intent(inout) :: lloc
551  integer(I4B), intent(inout) :: istart
552  integer(I4B), intent(inout) :: istop
553  integer(I4B), intent(in) :: in
554  integer(I4B), intent(in) :: iout
555  character(len=*), intent(inout) :: line
556  logical, optional, intent(in) :: flag_string
557  integer(I4B) :: noder
558  ! -- local
559  integer(I4B) :: nodeu
560  character(len=LINELENGTH) :: nodestr
561  logical :: flag_string_local
562  !
563  if (present(flag_string)) then
564  flag_string_local = flag_string
565  else
566  flag_string_local = .false.
567  end if
568  nodeu = this%nodeu_from_string(lloc, istart, istop, in, iout, line, &
569  flag_string_local)
570  !
571  ! -- Convert user-based nodenumber to reduced node number
572  if (nodeu > 0) then
573  noder = this%get_nodenumber(nodeu, 0)
574  else
575  noder = nodeu
576  end if
577  if (noder <= 0 .and. .not. flag_string_local) then
578  call this%nodeu_to_string(nodeu, nodestr)
579  write (errmsg, *) &
580  ' Cell is outside active grid domain: '// &
581  trim(adjustl(nodestr))
582  call store_error(errmsg)
583  end if
584  end function noder_from_string
585 
586  !> @brief Convert cellid string to reduced nodenumber
587  !!
588  !! If flag_string argument is present and true, the first token in string
589  !! is allowed to be a string (e.g. boundary name). In this case, if a string
590  !! is encountered, return value as -2.
591  !! If allow_zero argument is present and true, if all indices equal zero, the
592  !! result can be zero. If allow_zero is false, a zero in any index is an error.
593  !<
594  function noder_from_cellid(this, cellid, inunit, iout, flag_string, &
595  allow_zero) result(noder)
596  ! -- return
597  integer(I4B) :: noder
598  ! -- dummy
599  class(disbasetype) :: this
600  character(len=*), intent(inout) :: cellid
601  integer(I4B), intent(in) :: inunit
602  integer(I4B), intent(in) :: iout
603  logical, optional, intent(in) :: flag_string
604  logical, optional, intent(in) :: allow_zero
605  ! -- local
606  integer(I4B) :: nodeu
607  logical :: allowzerolocal
608  character(len=LINELENGTH) :: nodestr
609  logical :: flag_string_local
610  !
611  if (present(flag_string)) then
612  flag_string_local = flag_string
613  else
614  flag_string_local = .false.
615  end if
616  if (present(allow_zero)) then
617  allowzerolocal = allow_zero
618  else
619  allowzerolocal = .false.
620  end if
621  !
622  nodeu = this%nodeu_from_cellid(cellid, inunit, iout, flag_string_local, &
623  allowzerolocal)
624  !
625  ! -- Convert user-based nodenumber to reduced node number
626  if (nodeu > 0) then
627  noder = this%get_nodenumber(nodeu, 0)
628  else
629  noder = nodeu
630  end if
631  if (noder <= 0 .and. .not. flag_string_local) then
632  call this%nodeu_to_string(nodeu, nodestr)
633  write (errmsg, *) &
634  ' Cell is outside active grid domain: '// &
635  trim(adjustl(nodestr))
636  call store_error(errmsg)
637  end if
638  end function noder_from_cellid
639 
640  !> @brief Indicates whether the grid discretization supports layers.
641  logical function supports_layers(this)
642  class(disbasetype) :: this
643  supports_layers = .false.
644  call store_error('Programmer error: supports_layers must be overridden', &
645  terminate=.true.)
646  end function supports_layers
647 
648  !> @brief Return number of cells per layer.
649  !! This is nodes for a DISU grid, as there are no layers.
650  function get_ncpl(this)
651  integer(I4B) :: get_ncpl
652  class(disbasetype) :: this
653  get_ncpl = 0
654  call store_error('Programmer error: get_ncpl must be overridden', &
655  terminate=.true.)
656  end function get_ncpl
657 
658  !> @brief Return volume of cell n based on x value passed.
659  function get_cell_volume(this, n, x)
660  ! -- return
661  real(dp) :: get_cell_volume
662  ! -- dummy
663  class(disbasetype) :: this
664  integer(I4B), intent(in) :: n
665  real(dp), intent(in) :: x
666  ! -- local
667  real(dp) :: tp
668  real(dp) :: bt
669  real(dp) :: sat
670  real(dp) :: thick
671 
673  tp = this%top(n)
674  bt = this%bot(n)
675  sat = squadraticsaturation(tp, bt, x)
676  thick = (tp - bt) * sat
677  get_cell_volume = this%area(n) * thick
678  end function get_cell_volume
679 
680  !> @brief Get a 2D array of cell polygon vertices, in
681  !! clockwise order starting with the lower left corner.
682  subroutine get_polyverts(this, ic, polyverts, closed)
683  class(disbasetype), intent(inout) :: this
684  integer(I4B), intent(in) :: ic !< cell number (reduced)
685  real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices (column-major indexing)
686  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
687 
688  errmsg = 'Programmer error: get_polyverts must be overridden'
689  call store_error(errmsg, terminate=.true.)
690  end subroutine get_polyverts
691 
692  !> @brief Get the number of cell polygon vertices.
693  function get_npolyverts(this, ic, closed) result(npolyverts)
694  class(disbasetype), intent(inout) :: this
695  integer(I4B), intent(in) :: ic !< cell number (reduced)
696  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
697  integer(I4B) :: npolyverts
698  npolyverts = 0 ! suppress compiler warning
699  errmsg = 'Programmer error: get_npolyverts must be overridden'
700  call store_error(errmsg, terminate=.true.)
701  end function get_npolyverts
702 
703  !> @brief Get the maximum number of cell polygon vertices.
704  function get_max_npolyverts(this, closed) result(max_npolyverts)
705  class(disbasetype), intent(inout) :: this
706  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
707  integer(I4B) :: max_npolyverts
708  max_npolyverts = 0 ! suppress compiler warning
709  errmsg = 'Programmer error: get_max_npolyverts must be overridden'
710  call store_error(errmsg, terminate=.true.)
711  end function get_max_npolyverts
712 
713  !> @brief Read an integer array
714  subroutine read_int_array(this, line, lloc, istart, istop, iout, in, &
715  iarray, aname)
716  ! -- dummy
717  class(disbasetype), intent(inout) :: this
718  character(len=*), intent(inout) :: line
719  integer(I4B), intent(inout) :: lloc
720  integer(I4B), intent(inout) :: istart
721  integer(I4B), intent(inout) :: istop
722  integer(I4B), intent(in) :: in
723  integer(I4B), intent(in) :: iout
724  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: iarray
725  character(len=*), intent(in) :: aname
726 
727  errmsg = 'Programmer error: read_int_array must be overridden'
728  call store_error(errmsg, terminate=.true.)
729  end subroutine read_int_array
730 
731  !> @brief Read a double precision array
732  subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, &
733  darray, aname)
734  ! -- dummy
735  class(disbasetype), intent(inout) :: this
736  character(len=*), intent(inout) :: line
737  integer(I4B), intent(inout) :: lloc
738  integer(I4B), intent(inout) :: istart
739  integer(I4B), intent(inout) :: istop
740  integer(I4B), intent(in) :: in
741  integer(I4B), intent(in) :: iout
742  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray
743  character(len=*), intent(in) :: aname
744 
745  errmsg = 'Programmer error: read_dbl_array must be overridden'
746  call store_error(errmsg, terminate=.true.)
747  end subroutine read_dbl_array
748 
749  !> @brief Fill an integer array
750  subroutine fill_int_array(this, ibuff1, ibuff2)
751  ! -- dummy
752  class(disbasetype), intent(inout) :: this
753  integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibuff1
754  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: ibuff2
755  ! -- local
756  integer(I4B) :: nodeu
757  integer(I4B) :: noder
758 
759  do nodeu = 1, this%nodesuser
760  noder = this%get_nodenumber(nodeu, 0)
761  if (noder <= 0) cycle
762  ibuff2(noder) = ibuff1(nodeu)
763  end do
764  end subroutine fill_int_array
765 
766  !> @brief Fill a double precision array
767  subroutine fill_dbl_array(this, buff1, buff2)
768  ! -- dummy
769  class(disbasetype), intent(inout) :: this
770  real(DP), dimension(:), pointer, contiguous, intent(in) :: buff1
771  real(DP), dimension(:), pointer, contiguous, intent(inout) :: buff2
772  ! -- local
773  integer(I4B) :: nodeu
774  integer(I4B) :: noder
775 
776  do nodeu = 1, this%nodesuser
777  noder = this%get_nodenumber(nodeu, 0)
778  if (noder <= 0) cycle
779  buff2(noder) = buff1(nodeu)
780  end do
781  end subroutine fill_dbl_array
782 
783  !> @brief Read a list using the list reader.
784  !!
785  !! Convert user node numbers to reduced numbers.
786  !! Terminate if any nodenumbers are within an inactive domain.
787  !! Set up time series and multiply by iauxmultcol if it exists.
788  !! Write the list to iout if iprpak is set.
789  !<
790  subroutine read_list(this, line_reader, in, iout, iprpak, nlist, &
791  inamedbound, iauxmultcol, nodelist, rlist, auxvar, &
792  auxname, boundname, label, pkgname, tsManager, iscloc, &
793  indxconvertflux)
794  ! -- modules
799  use inputoutputmodule, only: urword
802  ! -- dummy
803  class(disbasetype) :: this
804  type(longlinereadertype), intent(inout) :: line_reader
805  integer(I4B), intent(in) :: in
806  integer(I4B), intent(in) :: iout
807  integer(I4B), intent(in) :: iprpak
808  integer(I4B), intent(inout) :: nlist
809  integer(I4B), intent(in) :: inamedbound
810  integer(I4B), intent(in) :: iauxmultcol
811  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: nodelist
812  real(DP), dimension(:, :), pointer, contiguous, intent(inout) :: rlist
813  real(DP), dimension(:, :), pointer, contiguous, intent(inout) :: auxvar
814  character(len=LENAUXNAME), dimension(:), intent(inout) :: auxname
815  character(len=LENBOUNDNAME), dimension(:), pointer, contiguous, &
816  intent(inout) :: boundname
817  character(len=*), intent(in) :: label
818  character(len=*), intent(in) :: pkgName
819  type(timeseriesmanagertype) :: tsManager
820  integer(I4B), intent(in) :: iscloc
821  integer(I4B), intent(in), optional :: indxconvertflux
822  ! -- local
823  integer(I4B) :: l
824  integer(I4B) :: nodeu, noder
825  character(len=LINELENGTH) :: nodestr
826  integer(I4B) :: ii, jj
827  real(DP), pointer :: bndElem => null()
828  type(listreadertype) :: lstrdobj
829  type(timeserieslinktype), pointer :: tsLinkBnd => null()
830  type(timeserieslinktype), pointer :: tsLinkAux => null()
831  !
832  ! -- Read the list
833  call lstrdobj%read_list(line_reader, in, iout, nlist, inamedbound, &
834  this%mshape, nodelist, rlist, auxvar, auxname, &
835  boundname, label)
836  !
837  ! -- Go through all locations where a text string was found instead of
838  ! a double precision value and make time-series links to rlist
839  if (lstrdobj%ntxtrlist > 0) then
840  do l = 1, lstrdobj%ntxtrlist
841  ii = lstrdobj%idxtxtrow(l)
842  jj = lstrdobj%idxtxtcol(l)
843  tslinkbnd => null()
844  bndelem => rlist(jj, ii)
845  call read_value_or_time_series(lstrdobj%txtrlist(l), ii, jj, bndelem, &
846  pkgname, 'BND', tsmanager, iprpak, &
847  tslinkbnd)
848  if (associated(tslinkbnd)) then
849  !
850  ! -- If iauxmultcol is active and this column is the column
851  ! to be scaled, then assign tsLinkBnd%RMultiplier to auxvar
852  ! multiplier
853  if (iauxmultcol > 0 .and. jj == iscloc) then
854  tslinkbnd%RMultiplier => auxvar(iauxmultcol, ii)
855  end if
856  !
857  ! -- If boundaries are named, save the name in the link
858  if (lstrdobj%inamedbound == 1) then
859  tslinkbnd%BndName = lstrdobj%boundname(tslinkbnd%IRow)
860  end if
861  !
862  ! -- if the value is a flux and needs to be converted to a flow
863  ! then set the tsLinkBnd appropriately
864  if (present(indxconvertflux)) then
865  if (indxconvertflux == jj) then
866  tslinkbnd%convertflux = .true.
867  nodeu = nodelist(ii)
868  noder = this%get_nodenumber(nodeu, 0)
869  tslinkbnd%CellArea = this%get_area(noder)
870  end if
871  end if
872  !
873  end if
874  end do
875  end if
876  !
877  ! -- Make time-series substitutions for auxvar
878  if (lstrdobj%ntxtauxvar > 0) then
879  do l = 1, lstrdobj%ntxtauxvar
880  ii = lstrdobj%idxtxtauxrow(l)
881  jj = lstrdobj%idxtxtauxcol(l)
882  tslinkaux => null()
883  bndelem => auxvar(jj, ii)
884  call read_value_or_time_series(lstrdobj%txtauxvar(l), ii, jj, bndelem, &
885  pkgname, 'AUX', tsmanager, iprpak, &
886  tslinkaux)
887  if (lstrdobj%inamedbound == 1) then
888  if (associated(tslinkaux)) then
889  tslinkaux%BndName = lstrdobj%boundname(tslinkaux%IRow)
890  end if
891  end if
892  end do
893  end if
894  !
895  ! -- Multiply rlist by the multiplier column in auxvar
896  if (iauxmultcol > 0) then
897  do l = 1, nlist
898  rlist(iscloc, l) = rlist(iscloc, l) * auxvar(iauxmultcol, l)
899  end do
900  end if
901  !
902  ! -- Write the list to iout if requested
903  if (iprpak /= 0) then
904  call lstrdobj%write_list()
905  end if
906  !
907  ! -- Convert user nodenumbers to reduced nodenumbers, if necessary.
908  ! Conversion to reduced nodenumbers must be done last, after the
909  ! list is written so that correct indices are written to the list.
910  if (this%nodes < this%nodesuser) then
911  do l = 1, nlist
912  nodeu = nodelist(l)
913  noder = this%get_nodenumber(nodeu, 0)
914  if (noder <= 0) then
915  call this%nodeu_to_string(nodeu, nodestr)
916  write (errmsg, *) &
917  ' Cell is outside active grid domain: '// &
918  trim(adjustl(nodestr))
919  call store_error(errmsg)
920  end if
921  nodelist(l) = noder
922  end do
923  !
924  ! -- Check for errors and terminate if encountered
925  if (count_errors() > 0) then
926  write (errmsg, *) count_errors(), ' errors encountered.'
927  call store_error(errmsg)
928  call store_error_unit(in)
929  end if
930  end if
931  end subroutine read_list
932 
933  !> @brief Read a 2d double array into col icolbnd of darray.
934  !!
935  !! For cells that are outside of the active domain,
936  !! do not copy the array value into darray.
937  !<
938  subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, &
939  icolbnd, aname, inunit, iout)
940  ! -- dummy
941  class(disbasetype) :: this
942  integer(I4B), intent(in) :: ncolbnd
943  integer(I4B), intent(in) :: maxbnd
944  integer(I4B), dimension(maxbnd) :: nodelist
945  real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray
946  integer(I4B), intent(in) :: icolbnd
947  character(len=*), intent(in) :: aname
948  integer(I4B), intent(in) :: inunit
949  integer(I4B), intent(in) :: iout
950 
951  errmsg = 'Programmer error: read_layer_array must be overridden'
952  call store_error(errmsg, terminate=.true.)
953  end subroutine read_layer_array
954 
955  !> @brief Record a double precision array.
956  !!
957  !! The array is written to a formatted or unformatted external file
958  !! depending on the arguments.
959  subroutine record_array(this, darray, iout, iprint, idataun, aname, &
960  cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
961  ! -- dummy
962  class(disbasetype), intent(inout) :: this
963  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record
964  integer(I4B), intent(in) :: iout !< ascii output unit number
965  integer(I4B), intent(in) :: iprint !< whether to print the array
966  integer(I4B), intent(in) :: idataun !< binary output unit number
967  character(len=*), intent(in) :: aname !< text descriptor
968  character(len=*), intent(in) :: cdatafmp !< write format
969  integer(I4B), intent(in) :: nvaluesp !< values per line
970  integer(I4B), intent(in) :: nwidthp !< number width
971  character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E)
972  real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain
973 
974  errmsg = 'Programmer error: record_array must be overridden'
975  call store_error(errmsg, terminate=.true.)
976  end subroutine record_array
977 
978  !> @brief Record a connection-based double precision array
979  subroutine record_connection_array(this, flowja, ibinun, iout)
980  ! -- dummy
981  class(disbasetype) :: this
982  real(DP), dimension(:), intent(in) :: flowja
983  integer(I4B), intent(in) :: ibinun
984  integer(I4B), intent(in) :: iout
985  ! -- local
986  character(len=16), dimension(1) :: text
987  ! -- data
988  data text(1)/' FLOW-JA-FACE'/
989 
990  ! -- write full ja array
991  call ubdsv1(kstp, kper, text(1), ibinun, flowja, size(flowja), 1, 1, &
992  iout, delt, pertim, totim)
993  end subroutine record_connection_array
994 
995  !> @brief Convert reduced node number to string (nodenumber), (k,j) or (k,i,j)
996  subroutine noder_to_string(this, noder, str)
997  ! -- dummy
998  class(disbasetype) :: this
999  integer(I4B), intent(in) :: noder
1000  character(len=*), intent(inout) :: str
1001  ! -- local
1002  integer(I4B) :: nodeu
1003 
1004  nodeu = this%get_nodeuser(noder)
1005  call this%nodeu_to_string(nodeu, str)
1006  end subroutine noder_to_string
1007 
1008  !> @brief Convert reduced node number to array (nodenumber), (k,j) or (k,i,j)
1009  subroutine noder_to_array(this, noder, arr)
1010  ! -- dummy
1011  class(disbasetype) :: this
1012  integer(I4B), intent(in) :: noder
1013  integer(I4B), dimension(:), intent(inout) :: arr
1014  ! -- local
1015  integer(I4B) :: nodeu
1016 
1017  nodeu = this%get_nodeuser(noder)
1018  call this%nodeu_to_array(nodeu, arr)
1019  end subroutine noder_to_array
1020 
1021  !> @brief Record list header for imeth=6
1022  subroutine record_srcdst_list_header(this, text, textmodel, textpackage, &
1023  dstmodel, dstpackage, naux, auxtxt, &
1024  ibdchn, nlist, iout)
1025  class(disbasetype) :: this
1026  character(len=16), intent(in) :: text
1027  character(len=16), intent(in) :: textmodel
1028  character(len=16), intent(in) :: textpackage
1029  character(len=16), intent(in) :: dstmodel
1030  character(len=16), intent(in) :: dstpackage
1031  integer(I4B), intent(in) :: naux
1032  character(len=16), dimension(:), intent(in) :: auxtxt
1033  integer(I4B), intent(in) :: ibdchn
1034  integer(I4B), intent(in) :: nlist
1035  integer(I4B), intent(in) :: iout
1036 
1037  errmsg = 'Programmer error: record_srcdst_list_header must be overridden'
1038  call store_error(errmsg, terminate=.true.)
1039  end subroutine record_srcdst_list_header
1040 
1041  !> @brief Record list header
1042  subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, &
1043  naux, aux, olconv, olconv2)
1044  ! -- dummy
1045  class(disbasetype) :: this
1046  integer(I4B), intent(in) :: ibdchn
1047  integer(I4B), intent(in) :: noder
1048  integer(I4B), intent(in) :: noder2
1049  real(DP), intent(in) :: q
1050  integer(I4B), intent(in) :: naux
1051  real(DP), dimension(naux), intent(in) :: aux
1052  logical, optional, intent(in) :: olconv
1053  logical, optional, intent(in) :: olconv2
1054  ! -- local
1055  logical :: lconv
1056  logical :: lconv2
1057  integer(I4B) :: nodeu
1058  integer(I4B) :: nodeu2
1059  !
1060  ! -- Use ubdsvb to write list header
1061  if (present(olconv)) then
1062  lconv = olconv
1063  else
1064  lconv = .true.
1065  end if
1066  if (lconv) then
1067  nodeu = this%get_nodeuser(noder)
1068  else
1069  nodeu = noder
1070  end if
1071  if (present(olconv2)) then
1072  lconv2 = olconv2
1073  else
1074  lconv2 = .true.
1075  end if
1076  if (lconv2) then
1077  nodeu2 = this%get_nodeuser(noder2)
1078  else
1079  nodeu2 = noder2
1080  end if
1081  call ubdsvd(ibdchn, nodeu, nodeu2, q, naux, aux)
1082  end subroutine record_srcdst_list_entry
1083 
1084  !> @brief Convert an integer array to nodelist.
1085  !!
1086  !! For DIS/DISV, the array is layer number, for DISU it's node number.
1087  !<
1088  subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname)
1089  class(disbasetype) :: this
1090  integer(I4B), intent(in) :: maxbnd
1091  integer(I4B), dimension(:), pointer, contiguous :: darray
1092  integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
1093  integer(I4B), intent(inout) :: nbound
1094  character(len=*), intent(in) :: aname
1095 
1096  errmsg = 'Programmer error: nlarray_to_nodelist must be overridden'
1097  call store_error(errmsg, terminate=.true.)
1098  end subroutine nlarray_to_nodelist
1099 
1100  !> @brief Find the first highest active cell beneath cell n
1101  subroutine highest_active(this, n, ibound)
1102  ! -- dummy
1103  class(disbasetype) :: this
1104  integer(I4B), intent(inout) :: n
1105  integer(I4B), dimension(:), intent(in) :: ibound
1106  ! -- locals
1107  integer(I4B) :: m, ii, iis
1108  logical done, bottomcell
1109  !
1110  ! -- Loop through connected cells until the highest active one (including a
1111  ! constant head cell) is found. Return that cell as n.
1112  done = .false.
1113  do while (.not. done)
1114  bottomcell = .true.
1115  cloop: do ii = this%con%ia(n) + 1, this%con%ia(n + 1) - 1
1116  m = this%con%ja(ii)
1117  iis = this%con%jas(ii)
1118  if (this%con%ihc(iis) == 0 .and. m > n) then
1119  !
1120  ! -- this cannot be a bottom cell
1121  bottomcell = .false.
1122  !
1123  ! -- vertical down
1124  if (ibound(m) /= 0) then
1125  n = m
1126  done = .true.
1127  exit cloop
1128  else
1129  n = m
1130  exit cloop
1131  end if
1132  end if
1133  end do cloop
1134  if (bottomcell) done = .true.
1135  end do
1136  end subroutine highest_active
1137 
1138  !> @brief Find the first saturated cell beneath cell n
1139  subroutine highest_saturated(this, n, sat)
1140  ! -- dummy
1141  class(disbasetype) :: this
1142  integer(I4B), intent(inout) :: n
1143  real(DP), dimension(:), intent(in) :: sat
1144  ! -- locals
1145  integer(I4B) :: m, ii, iis
1146  logical(LGP) :: is_done, bottomcell
1147  !
1148  ! -- Loop through connected cells until the highest saturated one (including a
1149  ! constant head cell) is found. Return that cell as n.
1150  is_done = .false.
1151  do while (.not. is_done)
1152  bottomcell = .true.
1153  cloop: do ii = this%con%ia(n) + 1, this%con%ia(n + 1) - 1
1154  m = this%con%ja(ii)
1155  iis = this%con%jas(ii)
1156  if (this%con%ihc(iis) == 0 .and. m > n) then
1157  !
1158  ! -- this cannot be a bottom cell
1159  bottomcell = .false.
1160  !
1161  ! -- vertical down
1162  if (sat(m) > dzero) then
1163  n = m
1164  is_done = .true.
1165  exit cloop
1166  else
1167  n = m
1168  exit cloop
1169  end if
1170  end if
1171  end do cloop
1172  if (bottomcell) is_done = .true.
1173  end do
1174  end subroutine highest_saturated
1175 
1176  !> @brief Return the cell area for the given node
1177  function get_area(this, node) result(area)
1178  class(disbasetype) :: this
1179  integer(I4B), intent(in) :: node !< reduced node number
1180  real(dp) :: area
1181 
1182  area = this%area(node)
1183  end function get_area
1184 
1185  !> @ brief Calculate the area factor for the cell connection
1186  !!
1187  !! Function calculates the area factor for the cell connection. The sum of
1188  !! all area factors for all cell connections to overlying or underlying
1189  !! cells cells will be 1.
1190  !!
1191  !! TODO: confirm that this works for cells that are only partially covered
1192  !! by overlying or underlying cells.
1193  !<
1194  function get_area_factor(this, node, idx_conn) result(area_factor)
1195  ! -- return
1196  real(dp) :: area_factor !< connection cell area factor
1197  ! -- dummy
1198  class(disbasetype) :: this
1199  integer(I4B), intent(in) :: node !< cell node number
1200  integer(I4B), intent(in) :: idx_conn !< connection index
1201  ! -- local
1202  real(dp) :: area_node
1203  real(dp) :: area_conn
1204  !
1205  ! -- calculate the cell area fraction
1206  area_node = this%area(node)
1207  area_conn = this%con%hwva(idx_conn)
1208  !
1209  ! -- return the cell area factor
1210  area_factor = area_conn / area_node
1211  end function get_area_factor
1212 
1213  !> @ brief Calculate the flow width between two cells
1214  !!
1215  !! This should only be called for connections with IHC > 0.
1216  !! Routine is needed, so it can be overridden by the linear
1217  !! network discretization, which allows for a separate flow
1218  !< width for each cell.
1219  !<
1220  subroutine get_flow_width(this, n, m, idx_conn, width_n, width_m)
1221  ! dummy
1222  class(disbasetype) :: this
1223  integer(I4B), intent(in) :: n !< cell node number
1224  integer(I4B), intent(in) :: m !< cell node number
1225  integer(I4B), intent(in) :: idx_conn !< connection index
1226  real(DP), intent(out) :: width_n !< flow width for cell n
1227  real(DP), intent(out) :: width_m !< flow width for cell m
1228  ! local
1229  integer(I4B) :: isympos
1230 
1231  ! For general case, width_n = width_m
1232  isympos = this%con%jas(idx_conn)
1233  width_n = this%con%hwva(isympos)
1234  width_m = width_n
1235 
1236  end subroutine get_flow_width
1237 
1238  !> @Brief return true if grid is three dimensional
1239  function is_3d(this) result(r)
1240  ! dummy
1241  class(disbasetype) :: this
1242  ! return
1243  logical(LGP) :: r
1244  r = .false.
1245  select case (this%get_dis_enum())
1246  case (dis, disv, disu)
1247  r = .true.
1248  end select
1249  end function is_3d
1250 
1251  !> @Brief return true if grid is two dimensional
1252  function is_2d(this) result(r)
1253  ! dummy
1254  class(disbasetype) :: this
1255  ! return
1256  logical(LGP) :: r
1257  r = .false.
1258  select case (this%get_dis_enum())
1259  case (dis2d, disv2d, disu2d)
1260  r = .true.
1261  end select
1262  end function is_2d
1263 
1264  !> @Brief return true if grid is one dimensional
1265  function is_1d(this) result(r)
1266  ! dummy
1267  class(disbasetype) :: this
1268  ! return
1269  logical(LGP) :: r
1270  r = .false.
1271  select case (this%get_dis_enum())
1272  case (dis1d, disv1d, disu1d)
1273  r = .true.
1274  end select
1275  end function is_1d
1276 
1277 end module basedismodule
real(dp) function get_area(this, node)
Return the cell area for the given node.
subroutine dis_df(this)
Define the discretization.
subroutine noder_to_string(this, noder, str)
Convert reduced node number to string (nodenumber), (k,j) or (k,i,j)
integer(i4b) function get_nodenumber_idx2(this, k, j, icheck)
real(dp) function get_area_factor(this, node, idx_conn)
@ brief Calculate the area factor for the cell connection
subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, icolbnd, aname, inunit, iout)
Read a 2d double array into col icolbnd of darray.
integer(i4b) function nodeu_from_string(this, lloc, istart, istop, in, iout, line, flag_string, allow_zero)
Convert a string to a user nodenumber.
subroutine allocate_scalars(this, name_model, input_mempath)
Allocate and initialize scalar variables.
subroutine dis_mc(this, moffset, idxglo, matrix_sln)
Map cell connections in the numerical solution coefficient matrix.
logical(lgp) function is_1d(this)
@Brief return true if grid is one dimensional
integer(i4b) function get_max_npolyverts(this, closed)
Get the maximum number of cell polygon vertices.
subroutine highest_active(this, n, ibound)
Find the first highest active cell beneath cell n.
subroutine record_array(this, darray, iout, iprint, idataun, aname, cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
Record a double precision array.
subroutine record_srcdst_list_header(this, text, textmodel, textpackage, dstmodel, dstpackage, naux, auxtxt, ibdchn, nlist, iout)
Record list header for imeth=6.
integer(i4b) function get_nodenumber_idx3(this, k, i, j, icheck)
subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, naux, aux, olconv, olconv2)
Record list header.
subroutine write_grb(this, icelltype)
Write a binary grid file.
subroutine, public dis_transform_xy(x, y, xorigin, yorigin, angrot, xglo, yglo)
Get global (x, y) coordinates from cell-local coordinates.
subroutine read_int_array(this, line, lloc, istart, istop, iout, in, iarray, aname)
Read an integer array.
logical(lgp) function is_3d(this)
@Brief return true if grid is three dimensional
logical function supports_layers(this)
Indicates whether the grid discretization supports layers.
subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, ipos)
Get normal vector components between the cell and a given neighbor. The normal points outward from th...
subroutine nodeu_to_string(this, nodeu, str)
Convert a user nodenumber to a string (nodenumber), (k,j), or (k,i,j)
integer(i4b) function noder_from_string(this, lloc, istart, istop, in, iout, line, flag_string)
Convert a string to a reduced nodenumber.
subroutine fill_int_array(this, ibuff1, ibuff2)
Fill an integer array.
subroutine get_polyverts(this, ic, polyverts, closed)
Get a 2D array of cell polygon vertices, in clockwise order starting with the lower left corner.
subroutine dis_da(this)
@brier Deallocate variables
subroutine fill_dbl_array(this, buff1, buff2)
Fill a double precision array.
subroutine dis_ac(this, moffset, sparse)
Add connections to sparse cell connectivity matrix.
subroutine record_connection_array(this, flowja, ibinun, iout)
Record a connection-based double precision array.
subroutine nodeu_to_array(this, nodeu, arr)
Convert a user nodenumber to an array (nodenumber), (k,j), or (k,i,j)
integer(i4b) function get_nodenumber_idx1(this, nodeu, icheck)
integer(i4b) function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, allow_zero)
Convert a cellid string to a user nodenumber.
subroutine highest_saturated(this, n, sat)
Find the first saturated cell beneath cell n.
logical(lgp) function is_2d(this)
@Brief return true if grid is two dimensional
subroutine noder_to_array(this, noder, arr)
Convert reduced node number to array (nodenumber), (k,j) or (k,i,j)
subroutine get_dis_type(this, dis_type)
Get the discretization type (DIS, DISV, or DISU)
integer(i4b) function get_nodeuser(this, noder)
Convert a reduced nodenumber to a user node number.
integer(i4b) function get_npolyverts(this, ic, closed)
Get the number of cell polygon vertices.
subroutine get_flow_width(this, n, m, idx_conn, width_n, width_m)
@ brief Calculate the flow width between two cells
subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, darray, aname)
Read a double precision array.
subroutine dis_ar(this, icelltype)
Allocate and setup variables, and write binary grid file.
subroutine allocate_arrays(this)
Allocate and initialize arrays.
integer(i4b) function get_dis_enum(this)
Get the discretization type enumeration.
integer(i4b) function noder_from_cellid(this, cellid, inunit, iout, flag_string, allow_zero)
Convert cellid string to reduced nodenumber.
subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname)
Convert an integer array to nodelist.
integer(i4b) function get_ncpl(this)
Return number of cells per layer. This is nodes for a DISU grid, as there are no layers.
real(dp) function get_cell_volume(this, n, x)
Return volume of cell n based on x value passed.
subroutine read_list(this, line_reader, in, iout, iprpak, nlist, inamedbound, iauxmultcol, nodelist, rlist, auxvar, auxname, boundname, label, pkgname, tsManager, iscloc, indxconvertflux)
Read a list using the list reader.
subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, xcomp, ycomp, zcomp, conlen)
Get unit vector components between the cell and a given neighbor. Saturation must be provided to comp...
This module contains block parser methods.
Definition: BlockParser.f90:7
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 lenmodelname
maximum length of the model name
Definition: Constants.f90:22
@ disu
DISV6 discretization.
Definition: Constants.f90:157
@ dis
DIS6 discretization.
Definition: Constants.f90:155
@ dis1d
DIS1D6 discretization.
Definition: Constants.f90:159
@ disv2d
DISV2D6 discretization.
Definition: Constants.f90:164
@ disv1d
DISV1D6 discretization.
Definition: Constants.f90:160
@ dis2d
DIS2D6 discretization.
Definition: Constants.f90:163
@ disv
DISU6 discretization.
Definition: Constants.f90:156
@ disu1d
DISU1D6 discretization.
Definition: Constants.f90:161
@ disundef
undefined discretization
Definition: Constants.f90:153
@ disu2d
DISU2D6 discretization.
Definition: Constants.f90:165
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 dpio180
real constant
Definition: Constants.f90:130
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
subroutine, public ubdsv1(kstp, kper, text, ibdchn, buff, ncol, nrow, nlay, iout, delt, pertim, totim)
Record cell-by-cell flow terms for one component of flow as a 3-D array with extra record to indicate...
subroutine, public ubdsvd(ibdchn, n, n2, q, naux, aux)
Write one value of cell-by-cell flow using a list structure.
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
Generic List Reader Module.
Definition: ListReader.f90:3
This module contains the LongLineReaderType.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
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
real(dp) function squadraticsaturation(top, bot, x, eps)
@ brief sQuadraticSaturation
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
subroutine, public read_value_or_time_series(textInput, ii, jj, bndElem, pkgName, auxOrBnd, tsManager, iprpak, tsLink)
Call this subroutine if the time-series link is available or needed.