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