MODFLOW 6  version 6.5.0.dev2
MODFLOW 6 Code Documentation
Dis2d.f90
Go to the documentation of this file.
1 module dis2dmodule
2 
4  use kindmodule, only: dp, i4b, lgp
5  use constantsmodule, only: linelength, dhalf, done, dzero, &
7  use basedismodule, only: disbasetype
16  use tdismodule, only: kstp, kper, pertim, totim, delt
17 
18  implicit none
19  private
20  public dis2d_cr, dis2dtype
21 
22  !> @brief Structured grid discretization
23  type, extends(disbasetype) :: dis2dtype
24  integer(I4B), pointer :: nrow => null() !< number of rows
25  integer(I4B), pointer :: ncol => null() !< number of columns
26  real(dp), dimension(:), pointer, contiguous :: delr => null() !< spacing along a row
27  real(dp), dimension(:), pointer, contiguous :: delc => null() !< spacing along a column
28  real(dp), dimension(:, :), pointer, contiguous :: botm => null() !< bottom elevations for each cell (ncol, nrow)
29  integer(I4B), dimension(:, :), pointer, contiguous :: idomain => null() !< idomain (ncol, nrow)
30  real(dp), dimension(:), pointer, contiguous :: cellx => null() !< cell center x coordinate for column j
31  real(dp), dimension(:), pointer, contiguous :: celly => null() !< cell center y coordinate for row i
32 
33  contains
34 
35  procedure :: dis_df => dis3d_df
36  procedure :: dis_da => dis3d_da
37  procedure :: get_dis_type => get_dis_type
38  procedure :: get_dis_enum => get_dis_enum
39  procedure, public :: record_array
40  procedure, public :: read_layer_array
41  procedure, public :: record_srcdst_list_header
42  procedure, public :: nlarray_to_nodelist
43  ! -- helper functions
44  procedure :: get_nodenumber_idx1
45  procedure :: get_nodenumber_idx2
46  procedure :: nodeu_to_string
47  procedure :: nodeu_to_array
48  procedure :: nodeu_from_string
49  procedure :: nodeu_from_cellid
50  procedure :: supports_layers
51  procedure :: get_ncpl
52  procedure :: get_polyverts
53  procedure :: connection_vector
54  procedure :: connection_normal
55  ! -- private
56  procedure :: source_options
57  procedure :: source_dimensions
58  procedure :: source_griddata
59  procedure :: log_options
60  procedure :: log_dimensions
61  procedure :: log_griddata
62  procedure :: grid_finalize
63  procedure :: write_grb
64  procedure :: allocate_scalars
65  procedure :: allocate_arrays
66  !
67  ! -- Read a node-sized model array (reduced or not)
68  procedure :: read_int_array
69  procedure :: read_dbl_array
70  end type dis2dtype
71 
72  !> @brief Simplifies tracking parameters sourced from the input context.
74  logical :: length_units = .false.
75  logical :: nogrb = .false.
76  logical :: xorigin = .false.
77  logical :: yorigin = .false.
78  logical :: angrot = .false.
79  logical :: nrow = .false.
80  logical :: ncol = .false.
81  logical :: delr = .false.
82  logical :: delc = .false.
83  logical :: top = .false.
84  logical :: botm = .false.
85  logical :: idomain = .false.
86  end type disfoundtype
87 
88 contains
89 
90  !> @brief Create a new structured discretization object
91  !<
92  subroutine dis2d_cr(dis, name_model, input_mempath, inunit, iout)
93  ! -- dummy
94  class(disbasetype), pointer :: dis
95  character(len=*), intent(in) :: name_model
96  character(len=*), intent(in) :: input_mempath
97  integer(I4B), intent(in) :: inunit
98  integer(I4B), intent(in) :: iout
99  ! -- locals
100  type(dis2dtype), pointer :: disnew
101  ! -- formats
102  character(len=*), parameter :: fmtheader = &
103  "(1X, /1X, 'DIS -- STRUCTURED GRID DISCRETIZATION PACKAGE,', &
104  &' VERSION 2 : 3/27/2014 - INPUT READ FROM MEMPATH: ', A, /)"
105  !
106  allocate (disnew)
107  dis => disnew
108  call disnew%allocate_scalars(name_model, input_mempath)
109  dis%inunit = inunit
110  dis%iout = iout
111  !
112  ! -- If dis enabled
113  if (inunit > 0) then
114  !
115  ! -- Identify package
116  if (iout > 0) then
117  write (iout, fmtheader) dis%input_mempath
118  end if
119  end if
120  !
121  end subroutine dis2d_cr
122 
123  !> @brief Define the discretization
124  !<
125  subroutine dis3d_df(this)
126  ! -- dummy
127  class(dis2dtype) :: this
128  !
129  ! -- Transfer the data from the memory manager into this package object
130  if (this%inunit /= 0) then
131  !
132  ! -- source input options
133  call this%source_options()
134  !
135  ! -- source input dimensions
136  call this%source_dimensions()
137  !
138  ! -- source input griddata
139  call this%source_griddata()
140  end if
141  !
142  ! -- Final grid initialization
143  call this%grid_finalize()
144  !
145  end subroutine dis3d_df
146 
147  !> @brief Deallocate variables
148  !<
149  subroutine dis3d_da(this)
150  ! -- dummy
151  class(dis2dtype) :: this
152  !
153  ! -- Deallocate idm memory
154  call memorylist_remove(this%name_model, 'DIS2D', idm_context)
155  !
156  ! -- DisBaseType deallocate
157  call this%DisBaseType%dis_da()
158  !
159  ! -- Deallocate scalars
160  call mem_deallocate(this%nrow)
161  call mem_deallocate(this%ncol)
162  call mem_deallocate(this%delr)
163  call mem_deallocate(this%delc)
164  call mem_deallocate(this%cellx)
165  call mem_deallocate(this%celly)
166  !
167  ! -- Deallocate Arrays
168  call mem_deallocate(this%nodereduced)
169  call mem_deallocate(this%nodeuser)
170  call mem_deallocate(this%botm)
171  call mem_deallocate(this%idomain)
172  !
173  end subroutine dis3d_da
174 
175  !> @brief Copy options from IDM into package
176  !<
177  subroutine source_options(this)
178  ! -- dummy
179  class(dis2dtype) :: this
180  ! -- locals
181  character(len=LENVARNAME), dimension(3) :: lenunits = &
182  &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS']
183  type(disfoundtype) :: found
184  !
185  ! -- update defaults with idm sourced values
186  call mem_set_value(this%lenuni, 'LENGTH_UNITS', this%input_mempath, &
187  lenunits, found%length_units)
188  call mem_set_value(this%nogrb, 'NOGRB', this%input_mempath, found%nogrb)
189  call mem_set_value(this%xorigin, 'XORIGIN', this%input_mempath, found%xorigin)
190  call mem_set_value(this%yorigin, 'YORIGIN', this%input_mempath, found%yorigin)
191  call mem_set_value(this%angrot, 'ANGROT', this%input_mempath, found%angrot)
192  !
193  ! -- log values to list file
194  if (this%iout > 0) then
195  call this%log_options(found)
196  end if
197  !
198  end subroutine source_options
199 
200  !> @brief Write user options to list file
201  !<
202  subroutine log_options(this, found)
203  ! -- dummy
204  class(dis2dtype) :: this
205  type(disfoundtype), intent(in) :: found
206  !
207  write (this%iout, '(1x,a)') 'Setting Discretization Options'
208  !
209  if (found%length_units) then
210  write (this%iout, '(4x,a,i0)') 'Model length unit [0=UND, 1=FEET, &
211  &2=METERS, 3=CENTIMETERS] set as ', this%lenuni
212  end if
213  !
214  if (found%nogrb) then
215  write (this%iout, '(4x,a,i0)') 'Binary grid file [0=GRB, 1=NOGRB] &
216  &set as ', this%nogrb
217  end if
218  !
219  if (found%xorigin) then
220  write (this%iout, '(4x,a,G0)') 'XORIGIN = ', this%xorigin
221  end if
222  !
223  if (found%yorigin) then
224  write (this%iout, '(4x,a,G0)') 'YORIGIN = ', this%yorigin
225  end if
226  !
227  if (found%angrot) then
228  write (this%iout, '(4x,a,G0)') 'ANGROT = ', this%angrot
229  end if
230  !
231  write (this%iout, '(1x,a,/)') 'End Setting Discretization Options'
232  !
233  end subroutine log_options
234 
235  !> @brief Copy dimensions from IDM into package
236  !<
237  subroutine source_dimensions(this)
238  ! -- dummy
239  class(dis2dtype) :: this
240  ! -- locals
241  integer(I4B) :: i, j
242  type(disfoundtype) :: found
243  !
244  ! -- update defaults with idm sourced values
245  call mem_set_value(this%nrow, 'NROW', this%input_mempath, found%nrow)
246  call mem_set_value(this%ncol, 'NCOL', this%input_mempath, found%ncol)
247  !
248  ! -- log simulation values
249  if (this%iout > 0) then
250  call this%log_dimensions(found)
251  end if
252  !
253  ! -- verify dimensions were set
254  if (this%nrow < 1) then
255  call store_error( &
256  'NROW was not specified or was specified incorrectly.')
257  call store_error_filename(this%input_fname)
258  end if
259  if (this%ncol < 1) then
260  call store_error( &
261  'NCOL was not specified or was specified incorrectly.')
262  call store_error_filename(this%input_fname)
263  end if
264  !
265  ! -- calculate nodesuser
266  this%nodesuser = this%nrow * this%ncol
267  !
268  ! -- Allocate delr, delc, and non-reduced vectors for dis
269  call mem_allocate(this%delr, this%ncol, 'DELR', this%memoryPath)
270  call mem_allocate(this%delc, this%nrow, 'DELC', this%memoryPath)
271  call mem_allocate(this%idomain, this%ncol, this%nrow, 'IDOMAIN', &
272  this%memoryPath)
273  call mem_allocate(this%botm, this%ncol, this%nrow, 'BOTM', &
274  this%memoryPath)
275  call mem_allocate(this%cellx, this%ncol, 'CELLX', this%memoryPath)
276  call mem_allocate(this%celly, this%nrow, 'CELLY', this%memoryPath)
277  !
278  ! -- initialize all cells to be active (idomain = 1)
279  do i = 1, this%nrow
280  do j = 1, this%ncol
281  this%idomain(j, i) = 1
282  end do
283  end do
284  !
285  end subroutine source_dimensions
286 
287  !> @brief Write dimensions to list file
288  !<
289  subroutine log_dimensions(this, found)
290  ! -- dummy
291  class(dis2dtype) :: this
292  type(disfoundtype), intent(in) :: found
293  !
294  write (this%iout, '(1x,a)') 'Setting Discretization Dimensions'
295  !
296  if (found%nrow) then
297  write (this%iout, '(4x,a,i0)') 'NROW = ', this%nrow
298  end if
299  !
300  if (found%ncol) then
301  write (this%iout, '(4x,a,i0)') 'NCOL = ', this%ncol
302  end if
303  !
304  write (this%iout, '(1x,a,/)') 'End Setting Discretization Dimensions'
305  !
306  end subroutine log_dimensions
307 
308  !> @brief Copy grid data from IDM into package
309  !<
310  subroutine source_griddata(this)
311  ! -- dummy
312  class(dis2dtype) :: this
313  type(disfoundtype) :: found
314  !
315  ! -- update defaults with idm sourced values
316  call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr)
317  call mem_set_value(this%delc, 'DELC', this%input_mempath, found%delc)
318  call mem_set_value(this%botm, 'BOTM', this%input_mempath, found%botm)
319  call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
320  !
321  ! -- log simulation values
322  if (this%iout > 0) then
323  call this%log_griddata(found)
324  end if
325  !
326  end subroutine source_griddata
327 
328  !> @brief Write dimensions to list file
329  !<
330  subroutine log_griddata(this, found)
331  ! -- dummy
332  class(dis2dtype) :: this
333  type(disfoundtype), intent(in) :: found
334  !
335  write (this%iout, '(1x,a)') 'Setting Discretization Griddata'
336  !
337  if (found%delr) then
338  write (this%iout, '(4x,a)') 'DELR set from input file'
339  end if
340  !
341  if (found%delc) then
342  write (this%iout, '(4x,a)') 'DELC set from input file'
343  end if
344  !
345  if (found%top) then
346  write (this%iout, '(4x,a)') 'TOP set from input file'
347  end if
348  !
349  if (found%botm) then
350  write (this%iout, '(4x,a)') 'BOTM set from input file'
351  end if
352  !
353  if (found%idomain) then
354  write (this%iout, '(4x,a)') 'IDOMAIN set from input file'
355  end if
356  !
357  write (this%iout, '(1x,a,/)') 'End Setting Discretization Griddata'
358  !
359  end subroutine log_griddata
360 
361  !> @brief Finalize grid (check properties, allocate arrays, compute connections)
362  !<
363  subroutine grid_finalize(this)
364  ! -- modules
366  ! -- dummy
367  class(dis2dtype) :: this
368  ! -- locals
369  integer(I4B) :: i, j
370  integer(I4B) :: node
371  integer(I4B) :: noder
372  integer(I4B) :: nrsize
373  ! -- formats
374  character(len=*), parameter :: fmtdz = &
375  "('CELL (',i0,',',i0,',',i0,') THICKNESS <= 0. ', &
376  &'TOP, BOT: ',2(1pg24.15))"
377  character(len=*), parameter :: fmtnr = &
378  "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',&
379  &/1x, 'Number of user nodes: ',I0,&
380  &/1X, 'Number of nodes in solution: ', I0, //)"
381  !
382  ! -- count active cells
383  this%nodes = 0
384  do i = 1, this%nrow
385  do j = 1, this%ncol
386  if (this%idomain(j, i) > 0) this%nodes = this%nodes + 1
387  end do
388  end do
389  !
390  ! -- Check to make sure nodes is a valid number
391  if (this%nodes == 0) then
392  call store_error('Model does not have any active nodes. &
393  &Ensure IDOMAIN array has some values greater &
394  &than zero.')
395  call store_error_filename(this%input_fname)
396  end if
397  !
398  ! -- Write message if reduced grid
399  if (this%nodes < this%nodesuser) then
400  write (this%iout, fmtnr) this%nodesuser, this%nodes
401  end if
402  !
403  ! -- Array size is now known, so allocate
404  call this%allocate_arrays()
405  !
406  ! -- Fill the nodereduced array with the reduced nodenumber, or
407  ! a negative number to indicate it is a pass-through cell, or
408  ! a zero to indicate that the cell is excluded from the
409  ! solution.
410  if (this%nodes < this%nodesuser) then
411  node = 1
412  noder = 1
413  do i = 1, this%nrow
414  do j = 1, this%ncol
415  if (this%idomain(j, i) > 0) then
416  this%nodereduced(node) = noder
417  noder = noder + 1
418  elseif (this%idomain(j, i) < 0) then
419  this%nodereduced(node) = -1
420  else
421  this%nodereduced(node) = 0
422  end if
423  node = node + 1
424  end do
425  end do
426  end if
427  !
428  ! -- allocate and fill nodeuser if a reduced grid
429  if (this%nodes < this%nodesuser) then
430  node = 1
431  noder = 1
432  do i = 1, this%nrow
433  do j = 1, this%ncol
434  if (this%idomain(j, i) > 0) then
435  this%nodeuser(noder) = node
436  noder = noder + 1
437  end if
438  node = node + 1
439  end do
440  end do
441  end if
442  !
443  ! -- fill x,y coordinate arrays
444  this%cellx(1) = dhalf * this%delr(1)
445  this%celly(this%nrow) = dhalf * this%delc(this%nrow)
446  do j = 2, this%ncol
447  this%cellx(j) = this%cellx(j - 1) + dhalf * this%delr(j - 1) + &
448  dhalf * this%delr(j)
449  end do
450  ! -- row number increases in negative y direction:
451  do i = this%nrow - 1, 1, -1
452  this%celly(i) = this%celly(i + 1) + dhalf * this%delc(i + 1) + &
453  dhalf * this%delc(i)
454  end do
455  !
456  ! -- Move botm into bot, and calculate area
457  node = 0
458  do i = 1, this%nrow
459  do j = 1, this%ncol
460  node = node + 1
461  noder = node
462  if (this%nodes < this%nodesuser) noder = this%nodereduced(node)
463  if (noder <= 0) cycle
464  this%bot(noder) = this%botm(j, i)
465  this%area(noder) = this%delr(j) * this%delc(i)
466  this%xc(noder) = this%cellx(j)
467  this%yc(noder) = this%celly(i)
468  end do
469  end do
470  !
471  ! -- create and fill the connections object
472  nrsize = 0
473  if (this%nodes < this%nodesuser) nrsize = this%nodes
474  allocate (this%con)
475  call this%con%disconnections(this%name_model, this%nodes, &
476  this%ncol, this%nrow, 1, &
477  nrsize, this%delr, this%delc, &
478  this%top, this%bot, this%nodereduced, &
479  this%nodeuser)
480  this%nja = this%con%nja
481  this%njas = this%con%njas
482  !
483  end subroutine grid_finalize
484 
485  !> @brief Write a binary grid file
486  !<
487  subroutine write_grb(this, icelltype)
488  ! -- modules
489  use openspecmodule, only: access, form
490  ! -- dummy
491  class(dis2dtype) :: this
492  integer(I4B), dimension(:), intent(in) :: icelltype
493  ! -- local
494  integer(I4B) :: iunit, ntxt
495  integer(I4B), parameter :: lentxt = 100
496  character(len=50) :: txthdr
497  character(len=lentxt) :: txt
498  character(len=LINELENGTH) :: fname
499  character(len=*), parameter :: fmtgrdsave = &
500  "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', &
501  &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)"
502  !
503  ! -- Initialize
504  ntxt = 14
505  !
506  ! -- Open the file
507  fname = trim(this%input_fname)//'.grb'
508  iunit = getunit()
509  write (this%iout, fmtgrdsave) iunit, trim(adjustl(fname))
510  call openfile(iunit, this%iout, trim(adjustl(fname)), 'DATA(BINARY)', &
511  form, access, 'REPLACE')
512  !
513  ! -- write header information
514  write (txthdr, '(a)') 'GRID DIS2D'
515  txthdr(50:50) = new_line('a')
516  write (iunit) txthdr
517  write (txthdr, '(a)') 'VERSION 1'
518  txthdr(50:50) = new_line('a')
519  write (iunit) txthdr
520  write (txthdr, '(a, i0)') 'NTXT ', ntxt
521  txthdr(50:50) = new_line('a')
522  write (iunit) txthdr
523  write (txthdr, '(a, i0)') 'LENTXT ', lentxt
524  txthdr(50:50) = new_line('a')
525  write (iunit) txthdr
526  !
527  ! -- write variable definitions
528  write (txt, '(3a, i0)') 'NCELLS ', 'INTEGER ', 'NDIM 0 # ', this%nodesuser
529  txt(lentxt:lentxt) = new_line('a')
530  write (iunit) txt
531  write (txt, '(3a, i0)') 'NROW ', 'INTEGER ', 'NDIM 0 # ', this%nrow
532  txt(lentxt:lentxt) = new_line('a')
533  write (iunit) txt
534  write (txt, '(3a, i0)') 'NCOL ', 'INTEGER ', 'NDIM 0 # ', this%ncol
535  txt(lentxt:lentxt) = new_line('a')
536  write (iunit) txt
537  write (txt, '(3a, i0)') 'NJA ', 'INTEGER ', 'NDIM 0 # ', this%nja
538  txt(lentxt:lentxt) = new_line('a')
539  write (iunit) txt
540  write (txt, '(3a, 1pg24.15)') 'XORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%xorigin
541  txt(lentxt:lentxt) = new_line('a')
542  write (iunit) txt
543  write (txt, '(3a, 1pg24.15)') 'YORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%yorigin
544  txt(lentxt:lentxt) = new_line('a')
545  write (iunit) txt
546  write (txt, '(3a, 1pg24.15)') 'ANGROT ', 'DOUBLE ', 'NDIM 0 # ', this%angrot
547  txt(lentxt:lentxt) = new_line('a')
548  write (iunit) txt
549  write (txt, '(3a, i0)') 'DELR ', 'DOUBLE ', 'NDIM 1 ', this%ncol
550  txt(lentxt:lentxt) = new_line('a')
551  write (iunit) txt
552  write (txt, '(3a, i0)') 'DELC ', 'DOUBLE ', 'NDIM 1 ', this%nrow
553  txt(lentxt:lentxt) = new_line('a')
554  write (iunit) txt
555  write (txt, '(3a, i0)') 'BOTM ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser
556  txt(lentxt:lentxt) = new_line('a')
557  write (iunit) txt
558  write (txt, '(3a, i0)') 'IA ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1
559  txt(lentxt:lentxt) = new_line('a')
560  write (iunit) txt
561  write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', size(this%con%jausr)
562  txt(lentxt:lentxt) = new_line('a')
563  write (iunit) txt
564  write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
565  txt(lentxt:lentxt) = new_line('a')
566  write (iunit) txt
567  write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
568  txt(lentxt:lentxt) = new_line('a')
569  write (iunit) txt
570  !
571  ! -- write data
572  write (iunit) this%nodesuser ! ncells
573  write (iunit) this%nrow ! nrow
574  write (iunit) this%ncol ! ncol
575  write (iunit) this%nja ! nja
576  write (iunit) this%xorigin ! xorigin
577  write (iunit) this%yorigin ! yorigin
578  write (iunit) this%angrot ! angrot
579  write (iunit) this%delr ! delr
580  write (iunit) this%delc ! delc
581  write (iunit) this%botm ! botm
582  write (iunit) this%con%iausr ! iausr
583  write (iunit) this%con%jausr ! jausr
584  write (iunit) this%idomain ! idomain
585  write (iunit) icelltype ! icelltype
586  !
587  ! -- Close the file
588  close (iunit)
589  !
590  end subroutine write_grb
591 
592  !> @brief Convert a user nodenumber to a string (nodenumber) or (k,i,j)
593  !<
594  subroutine nodeu_to_string(this, nodeu, str)
595  ! -- dummy
596  class(dis2dtype) :: this
597  integer(I4B), intent(in) :: nodeu
598  character(len=*), intent(inout) :: str
599  ! -- local
600  integer(I4B) :: i, j, k
601  character(len=10) :: istr, jstr
602  !
603  call get_ijk(nodeu, this%nrow, this%ncol, 1, i, j, k)
604  write (istr, '(i10)') i
605  write (jstr, '(i10)') j
606  str = '('//trim(adjustl(istr))//','// &
607  trim(adjustl(jstr))//')'
608  !
609  end subroutine nodeu_to_string
610 
611  !> @brief Convert a user nodenumber to an array (nodenumber) or (i,j)
612  !<
613  subroutine nodeu_to_array(this, nodeu, arr)
614  ! -- dummy
615  class(dis2dtype) :: this
616  integer(I4B), intent(in) :: nodeu
617  integer(I4B), dimension(:), intent(inout) :: arr
618  ! -- local
619  integer(I4B) :: isize
620  integer(I4B) :: i, j, k
621  !
622  ! -- check the size of arr
623  isize = size(arr)
624  if (isize /= this%ndim) then
625  write (errmsg, '(a,i0,a,i0,a)') &
626  'Program error: nodeu_to_array size of array (', isize, &
627  ') is not equal to the discretization dimension (', this%ndim, ')'
628  call store_error(errmsg, terminate=.true.)
629  end if
630  !
631  ! -- get k, i, j
632  call get_ijk(nodeu, this%nrow, this%ncol, 1, i, j, k)
633  !
634  ! -- fill array
635  arr(1) = i
636  arr(2) = j
637  !
638  end subroutine nodeu_to_array
639 
640  !> @brief Get reduced node number from user node number
641  !<
642  function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
643  ! -- return
644  integer(I4B) :: nodenumber
645  ! -- dummy
646  class(dis2dtype), intent(in) :: this
647  integer(I4B), intent(in) :: nodeu
648  integer(I4B), intent(in) :: icheck
649  !
650  ! -- check the node number if requested
651  if (icheck /= 0) then
652  !
653  ! -- If within valid range, convert to reduced nodenumber
654  if (nodeu < 1 .or. nodeu > this%nodesuser) then
655  write (errmsg, '(a,i0,a)') &
656  'Node number (', nodeu, &
657  ') less than 1 or greater than the number of nodes.'
658  call store_error(errmsg)
659  nodenumber = 0
660  else
661  nodenumber = nodeu
662  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
663  end if
664  else
665  nodenumber = nodeu
666  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
667  end if
668  !
669  end function get_nodenumber_idx1
670 
671  !> @brief Get reduced node number from layer, row and column indices
672  !<
673  function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber)
674  ! -- return
675  integer(I4B) :: nodenumber
676  ! -- dummy
677  class(dis2dtype), intent(in) :: this
678  integer(I4B), intent(in) :: k, j
679  integer(I4B), intent(in) :: icheck
680  ! -- local
681  integer(I4B) :: nodeu, i
682  ! formats
683  character(len=*), parameter :: fmterr = &
684  "('Error in structured-grid cell indices: row = ',i0,&
685  &', column = ',i0)"
686  !
687  i = k
688  nodeu = get_node(1, i, j, 1, this%nrow, this%ncol)
689  if (nodeu < 1) then
690  write (errmsg, fmterr) i, j
691  call store_error(errmsg, terminate=.true.)
692  end if
693  nodenumber = nodeu
694  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
695  !
696  ! -- check the node number if requested
697  if (icheck /= 0) then
698  !
699  if (i < 1 .or. i > this%nrow) &
700  call store_error('Row less than one or greater than nrow')
701  if (j < 1 .or. j > this%ncol) &
702  call store_error('Column less than one or greater than ncol')
703  !
704  ! -- Error if outside of range
705  if (nodeu < 1 .or. nodeu > this%nodesuser) then
706  write (errmsg, '(a,i0,a)') &
707  'Node number (', nodeu, ')less than 1 or greater than nodes.'
708  call store_error(errmsg)
709  end if
710  end if
711  !
712  end function get_nodenumber_idx2
713 
714  !> @brief Allocate and initialize scalar variables
715  !<
716  subroutine allocate_scalars(this, name_model, input_mempath)
717  ! -- dummy
718  class(dis2dtype) :: this
719  character(len=*), intent(in) :: name_model
720  character(len=*), intent(in) :: input_mempath
721  !
722  ! -- Allocate parent scalars
723  call this%DisBaseType%allocate_scalars(name_model, input_mempath)
724  !
725  ! -- Allocate
726  call mem_allocate(this%nrow, 'NROW', this%memoryPath)
727  call mem_allocate(this%ncol, 'NCOL', this%memoryPath)
728  !
729  ! -- Initialize
730  this%nrow = 0
731  this%ncol = 0
732  this%ndim = 2
733  !
734  end subroutine allocate_scalars
735 
736  !> @brief Allocate and initialize arrays
737  !<
738  subroutine allocate_arrays(this)
739  ! -- dummy
740  class(dis2dtype) :: this
741  !
742  ! -- Allocate arrays in DisBaseType (mshape, top, bot, area)
743  call this%DisBaseType%allocate_arrays()
744  !
745  ! -- Allocate arrays for DisType
746  if (this%nodes < this%nodesuser) then
747  call mem_allocate(this%nodeuser, this%nodes, 'NODEUSER', this%memoryPath)
748  call mem_allocate(this%nodereduced, this%nodesuser, 'NODEREDUCED', &
749  this%memoryPath)
750  else
751  call mem_allocate(this%nodeuser, 1, 'NODEUSER', this%memoryPath)
752  call mem_allocate(this%nodereduced, 1, 'NODEREDUCED', this%memoryPath)
753  end if
754  !
755  ! -- Initialize
756  this%mshape(1) = this%nrow
757  this%mshape(2) = this%ncol
758  !
759  end subroutine allocate_arrays
760 
761  !> @brief Convert a string to a user nodenumber
762  !!
763  !! Parse layer, row and column and return user nodenumber.
764  !! If flag_string is present and true, the first token may be
765  !! non-numeric (e.g. boundary name). In this case, return -2.
766  !<
767  function nodeu_from_string(this, lloc, istart, istop, in, iout, line, &
768  flag_string, allow_zero) result(nodeu)
769  ! -- dummy
770  class(dis2dtype) :: this
771  integer(I4B), intent(inout) :: lloc
772  integer(I4B), intent(inout) :: istart
773  integer(I4B), intent(inout) :: istop
774  integer(I4B), intent(in) :: in
775  integer(I4B), intent(in) :: iout
776  character(len=*), intent(inout) :: line
777  logical, optional, intent(in) :: flag_string
778  logical, optional, intent(in) :: allow_zero
779  integer(I4B) :: nodeu
780  ! -- local
781  integer(I4B) :: i, j, nrow, ncol
782  integer(I4B) :: lloclocal, ndum, istat, n
783  real(dp) :: r
784  !
785  if (present(flag_string)) then
786  if (flag_string) then
787  ! Check to see if first token in line can be read as an integer.
788  lloclocal = lloc
789  call urword(line, lloclocal, istart, istop, 1, ndum, r, iout, in)
790  read (line(istart:istop), *, iostat=istat) n
791  if (istat /= 0) then
792  ! First token in line is not an integer; return flag to this effect.
793  nodeu = -2
794  return
795  end if
796  end if
797  end if
798  !
799  nrow = this%mshape(1)
800  ncol = this%mshape(2)
801  !
802  call urword(line, lloc, istart, istop, 2, i, r, iout, in)
803  call urword(line, lloc, istart, istop, 2, j, r, iout, in)
804  !
805  if (i == 0 .and. j == 0) then
806  if (present(allow_zero)) then
807  if (allow_zero) then
808  nodeu = 0
809  return
810  end if
811  end if
812  end if
813  !
814  errmsg = ""
815  !
816  if (i < 1 .or. i > nrow) then
817  write (errmsg, '(a,1x,a,i0,a)') &
818  trim(adjustl(errmsg)), 'Row number in list (', i, &
819  ') is outside of the grid.'
820  end if
821  if (j < 1 .or. j > ncol) then
822  write (errmsg, '(a,1x,a,i0,a)') &
823  trim(adjustl(errmsg)), 'Column number in list (', j, &
824  ') is outside of the grid.'
825  end if
826  !
827  nodeu = get_node(1, i, j, 1, nrow, ncol)
828  !
829  if (nodeu < 1 .or. nodeu > this%nodesuser) then
830  write (errmsg, '(a,1x,a,i0,a)') &
831  trim(adjustl(errmsg)), &
832  "Node number in list (", nodeu, ") is outside of the grid. "// &
833  "Cell number cannot be determined in line '"// &
834  trim(adjustl(line))//"'."
835  end if
836  !
837  if (len_trim(adjustl(errmsg)) > 0) then
838  call store_error(errmsg)
839  call store_error_unit(in)
840  end if
841  !
842  end function nodeu_from_string
843 
844  !> @brief Convert a cellid string to a user nodenumber
845  !!
846  !! If flag_string is present and true, the first token may be
847  !! non-numeric (e.g. boundary name). In this case, return -2.
848  !!
849  !! If allow_zero is present and true, and all indices are zero, the
850  !! result can be zero. If allow_zero is false, a zero in any index is an error.
851  !<
852  function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, &
853  allow_zero) result(nodeu)
854  ! -- return
855  integer(I4B) :: nodeu
856  ! -- dummy
857  class(dis2dtype) :: this
858  character(len=*), intent(inout) :: cellid
859  integer(I4B), intent(in) :: inunit
860  integer(I4B), intent(in) :: iout
861  logical, optional, intent(in) :: flag_string
862  logical, optional, intent(in) :: allow_zero
863  ! -- local
864  integer(I4B) :: lloclocal, istart, istop, ndum, n
865  integer(I4B) :: i, j, nrow, ncol
866  integer(I4B) :: istat
867  real(dp) :: r
868  !
869  if (present(flag_string)) then
870  if (flag_string) then
871  ! Check to see if first token in cellid can be read as an integer.
872  lloclocal = 1
873  call urword(cellid, lloclocal, istart, istop, 1, ndum, r, iout, inunit)
874  read (cellid(istart:istop), *, iostat=istat) n
875  if (istat /= 0) then
876  ! First token in cellid is not an integer; return flag to this effect.
877  nodeu = -2
878  return
879  end if
880  end if
881  end if
882  !
883  nrow = this%mshape(1)
884  ncol = this%mshape(2)
885  !
886  lloclocal = 1
887  call urword(cellid, lloclocal, istart, istop, 2, i, r, iout, inunit)
888  call urword(cellid, lloclocal, istart, istop, 2, j, r, iout, inunit)
889  !
890  if (i == 0 .and. j == 0) then
891  if (present(allow_zero)) then
892  if (allow_zero) then
893  nodeu = 0
894  return
895  end if
896  end if
897  end if
898  !
899  errmsg = ""
900  !
901  if (i < 1 .or. i > nrow) then
902  write (errmsg, '(a,1x,a,i0,a)') &
903  trim(adjustl(errmsg)), 'Row number in list (', i, &
904  ') is outside of the grid.'
905  end if
906  if (j < 1 .or. j > ncol) then
907  write (errmsg, '(a,1x,a,i0,a)') &
908  trim(adjustl(errmsg)), 'Column number in list (', j, &
909  ') is outside of the grid.'
910  end if
911  !
912  nodeu = get_node(1, i, j, 1, nrow, ncol)
913  !
914  if (nodeu < 1 .or. nodeu > this%nodesuser) then
915  write (errmsg, '(a,1x,a,i0,a)') &
916  trim(adjustl(errmsg)), &
917  "Cell number cannot be determined for cellid ("// &
918  trim(adjustl(cellid))//") and results in a user "// &
919  "node number (", nodeu, ") that is outside of the grid."
920  end if
921  !
922  if (len_trim(adjustl(errmsg)) > 0) then
923  call store_error(errmsg)
924  call store_error_unit(inunit)
925  end if
926  !
927  end function nodeu_from_cellid
928 
929  !> @brief Indicates whether the grid discretization supports layers
930  !<
931  logical function supports_layers(this)
932  ! -- dummy
933  class(dis2dtype) :: this
934  !
935  supports_layers = .true.
936  !
937  end function supports_layers
938 
939  !> @brief Return number of cells per layer (nrow * ncol)
940  !<
941  function get_ncpl(this)
942  integer(I4B) :: get_ncpl
943  class(dis2dtype) :: this
944  get_ncpl = this%nrow * this%ncol
945  end function get_ncpl
946 
947  !> @brief Get normal vector components between the cell and a given neighbor
948  !!
949  !! The normal points outward from the shared face between noden and nodem.
950  !<
951  subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, &
952  ipos)
953  ! -- dummy
954  class(dis2dtype) :: this
955  integer(I4B), intent(in) :: noden !< cell (reduced nn)
956  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
957  integer(I4B), intent(in) :: ihc !< horizontal connection flag
958  real(DP), intent(inout) :: xcomp
959  real(DP), intent(inout) :: ycomp
960  real(DP), intent(inout) :: zcomp
961  integer(I4B), intent(in) :: ipos
962  ! -- local
963  integer(I4B) :: nodeu1, i1, j1, k1
964  integer(I4B) :: nodeu2, i2, j2, k2
965  !
966  ! -- Set vector components based on ihc
967  if (ihc == 0) then
968  xcomp = dzero
969  ycomp = dzero
970  if (nodem < noden) then
971  !
972  ! -- nodem must be above noden, so upward connection
973  zcomp = done
974  else
975  !
976  ! -- nodem must be below noden, so downward connection
977  zcomp = -done
978  end if
979  else
980  xcomp = dzero
981  ycomp = dzero
982  zcomp = dzero
983  nodeu1 = this%get_nodeuser(noden)
984  nodeu2 = this%get_nodeuser(nodem)
985  call get_ijk(nodeu1, this%nrow, this%ncol, 1, i1, j1, k1)
986  call get_ijk(nodeu2, this%nrow, this%ncol, 1, i2, j2, k2)
987  if (i2 < i1) then ! back
988  ycomp = done
989  elseif (j2 < j1) then ! left
990  xcomp = -done
991  elseif (j2 > j1) then ! right
992  xcomp = done
993  else ! front
994  ycomp = -done
995  end if
996  !
997  end if
998  !
999  end subroutine connection_normal
1000 
1001  !> @brief Get unit vector components between the cell and a given neighbor
1002  !!
1003  !< For dis2d there is no z component
1004  subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, &
1005  xcomp, ycomp, zcomp, conlen)
1006  ! modules
1007  use disvgeom, only: line_unit_vector
1008  ! dummy
1009  class(dis2dtype) :: this
1010  integer(I4B), intent(in) :: noden !< cell (reduced nn)
1011  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
1012  logical, intent(in) :: nozee !< not used for dis2d
1013  real(DP), intent(in) :: satn !< not used for dis2d
1014  real(DP), intent(in) :: satm !< not used for dis2d
1015  integer(I4B), intent(in) :: ihc !< not used for dis2d (always horizontal)
1016  real(DP), intent(inout) :: xcomp !< x component of the connection vector
1017  real(DP), intent(inout) :: ycomp !< y component of the connection vector
1018  real(DP), intent(inout) :: zcomp !< z component, which is always zero
1019  real(DP), intent(inout) :: conlen !< calculated connection length
1020  ! local
1021  real(DP) :: z1, z2
1022  real(DP) :: x1, y1, x2, y2
1023  real(DP) :: ds
1024  integer(I4B) :: i1, i2, j1, j2, k1, k2
1025  integer(I4B) :: nodeu1, nodeu2, ipos
1026 
1027  ! Calculate vector components
1028  z1 = dzero
1029  z2 = dzero
1030  ipos = this%con%getjaindex(noden, nodem)
1031  ds = this%con%cl1(this%con%jas(ipos)) + this%con%cl2(this%con%jas(ipos))
1032  nodeu1 = this%get_nodeuser(noden)
1033  nodeu2 = this%get_nodeuser(nodem)
1034  call get_ijk(nodeu1, this%nrow, this%ncol, 1, i1, j1, k1)
1035  call get_ijk(nodeu2, this%nrow, this%ncol, 1, i2, j2, k2)
1036  x1 = dzero
1037  x2 = dzero
1038  y1 = dzero
1039  y2 = dzero
1040  if (i2 < i1) then ! back
1041  y2 = ds
1042  elseif (j2 < j1) then ! left
1043  x2 = -ds
1044  elseif (j2 > j1) then ! right
1045  x2 = ds
1046  else ! front
1047  y2 = -ds
1048  end if
1049  call line_unit_vector(x1, y1, z1, x2, y2, z2, xcomp, ycomp, zcomp, conlen)
1050  end subroutine connection_vector
1051 
1052  !> @brief Get the discretization type
1053  !<
1054  subroutine get_dis_type(this, dis_type)
1055  ! -- dummy
1056  class(dis2dtype), intent(in) :: this
1057  character(len=*), intent(out) :: dis_type
1058  !
1059  dis_type = "DIS2D"
1060  !
1061  end subroutine get_dis_type
1062 
1063  !> @brief Get the discretization type enumeration
1064  function get_dis_enum(this) result(dis_enum)
1065  use constantsmodule, only: dis2d
1066  class(dis2dtype), intent(in) :: this
1067  integer(I4B) :: dis_enum
1068  dis_enum = dis2d
1069  end function get_dis_enum
1070 
1071  !> @brief Get a 2D array of polygon vertices, listed in
1072  !!
1073  !! clockwise order beginning with the lower left corner
1074  !<
1075  subroutine get_polyverts(this, ic, polyverts, closed)
1076  ! -- dummy
1077  class(dis2dtype), intent(inout) :: this
1078  integer(I4B), intent(in) :: ic !< cell number (reduced)
1079  real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices (column-major indexing)
1080  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
1081  ! -- local
1082  integer(I4B) :: icu, nverts, irow, jcol, klay
1083  real(DP) :: cellx, celly, dxhalf, dyhalf
1084  logical(LGP) :: lclosed
1085  !
1086  nverts = 4
1087  !
1088  ! check closed option
1089  if (.not. (present(closed))) then
1090  lclosed = .false.
1091  else
1092  lclosed = closed
1093  end if
1094  !
1095  ! allocate vertices array
1096  if (lclosed) then
1097  allocate (polyverts(2, nverts + 1))
1098  else
1099  allocate (polyverts(2, nverts))
1100  end if
1101  !
1102  ! set vertices
1103  icu = this%get_nodeuser(ic)
1104  call get_ijk(icu, this%nrow, this%ncol, 1, irow, jcol, klay)
1105  cellx = this%cellx(jcol)
1106  celly = this%celly(irow)
1107  dxhalf = dhalf * this%delr(jcol)
1108  dyhalf = dhalf * this%delc(irow)
1109  polyverts(:, 1) = (/cellx - dxhalf, celly - dyhalf/) ! SW
1110  polyverts(:, 2) = (/cellx - dxhalf, celly + dyhalf/) ! NW
1111  polyverts(:, 3) = (/cellx + dxhalf, celly + dyhalf/) ! NE
1112  polyverts(:, 4) = (/cellx + dxhalf, celly - dyhalf/) ! SE
1113  !
1114  ! close if enabled
1115  if (lclosed) &
1116  polyverts(:, nverts + 1) = polyverts(:, 1)
1117  !
1118  end subroutine
1119 
1120  !> @brief Read an integer array
1121  !< TODO: REMOVE?
1122  subroutine read_int_array(this, line, lloc, istart, istop, iout, in, &
1123  iarray, aname)
1124  ! -- dummy
1125  class(dis2dtype), intent(inout) :: this
1126  character(len=*), intent(inout) :: line
1127  integer(I4B), intent(inout) :: lloc
1128  integer(I4B), intent(inout) :: istart
1129  integer(I4B), intent(inout) :: istop
1130  integer(I4B), intent(in) :: in
1131  integer(I4B), intent(in) :: iout
1132  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: iarray
1133  character(len=*), intent(in) :: aname
1134  ! -- local
1135 
1136  !
1137  end subroutine read_int_array
1138 
1139  !> @brief Read a double precision array
1140  !< TODO: REMOVE?
1141  subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, &
1142  darray, aname)
1143  ! -- dummy
1144  class(dis2dtype), intent(inout) :: this
1145  character(len=*), intent(inout) :: line
1146  integer(I4B), intent(inout) :: lloc
1147  integer(I4B), intent(inout) :: istart
1148  integer(I4B), intent(inout) :: istop
1149  integer(I4B), intent(in) :: in
1150  integer(I4B), intent(in) :: iout
1151  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray
1152  character(len=*), intent(in) :: aname
1153  !
1154  end subroutine read_dbl_array
1155 
1156  !> @brief Read a 2d double array into col icolbnd of darray
1157  !!
1158  !! For cells that are outside of the active domain,
1159  !! do not copy the array value into darray.
1160  !< TODO: REMOVE?
1161  subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, &
1162  icolbnd, aname, inunit, iout)
1163  ! -- dummy
1164  class(dis2dtype) :: this
1165  integer(I4B), intent(in) :: maxbnd
1166  integer(I4B), dimension(maxbnd) :: nodelist
1167  integer(I4B), intent(in) :: ncolbnd
1168  real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray
1169  integer(I4B), intent(in) :: icolbnd
1170  character(len=*), intent(in) :: aname
1171  integer(I4B), intent(in) :: inunit
1172  integer(I4B), intent(in) :: iout
1173  !
1174  end subroutine read_layer_array
1175 
1176  !> @brief Record a double precision array.
1177  !!
1178  !! The array is written to a formatted or unformatted external file
1179  !! depending on the arguments.
1180  !<
1181  subroutine record_array(this, darray, iout, iprint, idataun, aname, &
1182  cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
1183  ! -- dummy
1184  class(dis2dtype), intent(inout) :: this
1185  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record
1186  integer(I4B), intent(in) :: iout !< ascii output unit number
1187  integer(I4B), intent(in) :: iprint !< whether to print the array
1188  integer(I4B), intent(in) :: idataun !< binary output unit number, if negative don't write by layers, write entire array
1189  character(len=*), intent(in) :: aname !< text descriptor
1190  character(len=*), intent(in) :: cdatafmp !< write format
1191  integer(I4B), intent(in) :: nvaluesp !< values per line
1192  integer(I4B), intent(in) :: nwidthp !< number width
1193  character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E)
1194  real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain
1195  ! -- local
1196  integer(I4B) :: k, ifirst
1197  integer(I4B) :: nlay
1198  integer(I4B) :: nrow
1199  integer(I4B) :: ncol
1200  integer(I4B) :: nval
1201  integer(I4B) :: nodeu, noder
1202  integer(I4B) :: istart, istop
1203  real(DP), dimension(:), pointer, contiguous :: dtemp
1204  ! -- formats
1205  character(len=*), parameter :: fmthsv = &
1206  "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, &
1207  &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)"
1208  !
1209  ! -- set variables
1210  nlay = 1
1211  nrow = this%mshape(1)
1212  ncol = this%mshape(2)
1213  !
1214  ! -- If this is a reduced model, then copy the values from darray into
1215  ! dtemp.
1216  if (this%nodes < this%nodesuser) then
1217  nval = this%nodes
1218  dtemp => this%dbuff
1219  do nodeu = 1, this%nodesuser
1220  noder = this%get_nodenumber(nodeu, 0)
1221  if (noder <= 0) then
1222  dtemp(nodeu) = dinact
1223  cycle
1224  end if
1225  dtemp(nodeu) = darray(noder)
1226  end do
1227  else
1228  nval = this%nodes
1229  dtemp => darray
1230  end if
1231  !
1232  ! -- Print to iout if iprint /= 0
1233  if (iprint /= 0) then
1234  istart = 1
1235  do k = 1, nlay
1236  istop = istart + nrow * ncol - 1
1237  call ulaprufw(ncol, nrow, kstp, kper, k, iout, dtemp(istart:istop), &
1238  aname, cdatafmp, nvaluesp, nwidthp, editdesc)
1239  istart = istop + 1
1240  end do
1241  end if
1242  !
1243  ! -- Save array to an external file.
1244  if (idataun > 0) then
1245  ! -- write to binary file by layer
1246  ifirst = 1
1247  istart = 1
1248  do k = 1, nlay
1249  istop = istart + nrow * ncol - 1
1250  if (ifirst == 1) write (iout, fmthsv) &
1251  trim(adjustl(aname)), idataun, &
1252  kstp, kper
1253  ifirst = 0
1254  call ulasav(dtemp(istart:istop), aname, kstp, kper, &
1255  pertim, totim, ncol, nrow, k, idataun)
1256  istart = istop + 1
1257  end do
1258  elseif (idataun < 0) then
1259  !
1260  ! -- write entire array as one record
1261  call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, &
1262  iout, delt, pertim, totim)
1263  end if
1264  !
1265  end subroutine record_array
1266 
1267  !> @brief Record list header for imeth=6
1268  !<
1269  subroutine record_srcdst_list_header(this, text, textmodel, textpackage, &
1270  dstmodel, dstpackage, naux, auxtxt, &
1271  ibdchn, nlist, iout)
1272  ! -- dummy
1273  class(dis2dtype) :: this
1274  character(len=16), intent(in) :: text
1275  character(len=16), intent(in) :: textmodel
1276  character(len=16), intent(in) :: textpackage
1277  character(len=16), intent(in) :: dstmodel
1278  character(len=16), intent(in) :: dstpackage
1279  integer(I4B), intent(in) :: naux
1280  character(len=16), dimension(:), intent(in) :: auxtxt
1281  integer(I4B), intent(in) :: ibdchn
1282  integer(I4B), intent(in) :: nlist
1283  integer(I4B), intent(in) :: iout
1284  ! -- local
1285  integer(I4B) :: nlay, nrow, ncol
1286  !
1287  nlay = 1
1288  nrow = this%mshape(1)
1289  ncol = this%mshape(2)
1290  !
1291  ! -- Use ubdsv06 to write list header
1292  call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, &
1293  ibdchn, naux, auxtxt, ncol, nrow, nlay, &
1294  nlist, iout, delt, pertim, totim)
1295  !
1296  end subroutine record_srcdst_list_header
1297 
1298  !> @brief Convert an integer array (layer numbers) to nodelist
1299  !<
1300  subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname)
1301  ! -- dummy
1302  class(dis2dtype) :: this
1303  integer(I4B), intent(in) :: maxbnd
1304  integer(I4B), dimension(:), pointer, contiguous :: darray
1305  integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
1306  integer(I4B), intent(inout) :: nbound
1307  character(len=*), intent(in) :: aname
1308  ! -- local
1309  integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr
1310  !
1311  ! -- set variables
1312  nlay = 1
1313  nrow = this%mshape(1)
1314  ncol = this%mshape(2)
1315  !
1316  if (this%ndim > 1) then
1317  !
1318  nval = ncol * nrow
1319  !
1320  ! -- Copy array into nodelist
1321  ipos = 1
1322  ierr = 0
1323  do ir = 1, nrow
1324  do ic = 1, ncol
1325  nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
1326  il = darray(nodeu)
1327  if (il < 1 .or. il > nlay) then
1328  write (errmsg, '(a,1x,i0)') 'Invalid layer number:', il
1329  call store_error(errmsg, terminate=.true.)
1330  end if
1331  nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
1332  noder = this%get_nodenumber(nodeu, 0)
1333  if (ipos > maxbnd) then
1334  ierr = ipos
1335  else
1336  nodelist(ipos) = noder
1337  end if
1338  ipos = ipos + 1
1339  end do
1340  end do
1341  !
1342  ! -- Check for errors
1343  nbound = ipos - 1
1344  if (ierr > 0) then
1345  write (errmsg, '(a,1x,i0)') &
1346  'MAXBOUND dimension is too small.'// &
1347  'INCREASE MAXBOUND TO:', ierr
1348  call store_error(errmsg, terminate=.true.)
1349  end if
1350  !
1351  ! -- If nbound < maxbnd, then initialize nodelist to zero in this range
1352  if (nbound < maxbnd) then
1353  do ipos = nbound + 1, maxbnd
1354  nodelist(ipos) = 0
1355  end do
1356  end if
1357  !
1358  else
1359  !
1360  ! -- For unstructured, read nodelist directly, then check node numbers
1361  nodelist = darray
1362  do noder = 1, maxbnd
1363  if (noder < 1 .or. noder > this%nodes) then
1364  write (errmsg, '(a,1x,i0)') 'Invalid node number:', noder
1365  call store_error(errmsg, terminate=.true.)
1366  end if
1367  end do
1368  nbound = maxbnd
1369  !
1370  end if
1371  !
1372  end subroutine nlarray_to_nodelist
1373 
1374 end module dis2dmodule
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:44
@ dis2d
DIS2D6 discretization.
Definition: Constants.f90:162
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:67
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:64
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:26
real(dp), parameter done
real constant 1
Definition: Constants.f90:75
subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, icolbnd, aname, inunit, iout)
Read a 2d double array into col icolbnd of darray.
Definition: Dis2d.f90:1163
subroutine read_int_array(this, line, lloc, istart, istop, iout, in, iarray, aname)
Read an integer array.
Definition: Dis2d.f90:1124
integer(i4b) function get_nodenumber_idx2(this, k, j, icheck)
Get reduced node number from layer, row and column indices.
Definition: Dis2d.f90:674
subroutine log_options(this, found)
Write user options to list file.
Definition: Dis2d.f90:203
subroutine, public dis2d_cr(dis, name_model, input_mempath, inunit, iout)
Create a new structured discretization object.
Definition: Dis2d.f90:93
subroutine nodeu_to_string(this, nodeu, str)
Convert a user nodenumber to a string (nodenumber) or (k,i,j)
Definition: Dis2d.f90:595
subroutine record_srcdst_list_header(this, text, textmodel, textpackage, dstmodel, dstpackage, naux, auxtxt, ibdchn, nlist, iout)
Record list header for imeth=6.
Definition: Dis2d.f90:1272
subroutine get_dis_type(this, dis_type)
Get the discretization type.
Definition: Dis2d.f90:1055
integer(i4b) function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, allow_zero)
Convert a cellid string to a user nodenumber.
Definition: Dis2d.f90:854
integer(i4b) function get_dis_enum(this)
Get the discretization type enumeration.
Definition: Dis2d.f90:1065
subroutine dis3d_da(this)
Deallocate variables.
Definition: Dis2d.f90:150
subroutine nodeu_to_array(this, nodeu, arr)
Convert a user nodenumber to an array (nodenumber) or (i,j)
Definition: Dis2d.f90:614
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.
Definition: Dis2d.f90:1006
subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, darray, aname)
Read a double precision array.
Definition: Dis2d.f90:1143
subroutine source_options(this)
Copy options from IDM into package.
Definition: Dis2d.f90:178
subroutine log_griddata(this, found)
Write dimensions to list file.
Definition: Dis2d.f90:331
integer(i4b) function get_ncpl(this)
Return number of cells per layer (nrow * ncol)
Definition: Dis2d.f90:942
subroutine record_array(this, darray, iout, iprint, idataun, aname, cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
Record a double precision array.
Definition: Dis2d.f90:1183
subroutine allocate_scalars(this, name_model, input_mempath)
Allocate and initialize scalar variables.
Definition: Dis2d.f90:717
subroutine grid_finalize(this)
Finalize grid (check properties, allocate arrays, compute connections)
Definition: Dis2d.f90:364
logical function supports_layers(this)
Indicates whether the grid discretization supports layers.
Definition: Dis2d.f90:932
subroutine log_dimensions(this, found)
Write dimensions to list file.
Definition: Dis2d.f90:290
subroutine allocate_arrays(this)
Allocate and initialize arrays.
Definition: Dis2d.f90:739
integer(i4b) function get_nodenumber_idx1(this, nodeu, icheck)
Get reduced node number from user node number.
Definition: Dis2d.f90:643
subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, ipos)
Get normal vector components between the cell and a given neighbor.
Definition: Dis2d.f90:953
integer(i4b) function nodeu_from_string(this, lloc, istart, istop, in, iout, line, flag_string, allow_zero)
Convert a string to a user nodenumber.
Definition: Dis2d.f90:769
subroutine source_dimensions(this)
Copy dimensions from IDM into package.
Definition: Dis2d.f90:238
subroutine source_griddata(this)
Copy grid data from IDM into package.
Definition: Dis2d.f90:311
subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname)
Convert an integer array (layer numbers) to nodelist.
Definition: Dis2d.f90:1301
subroutine write_grb(this, icelltype)
Write a binary grid file.
Definition: Dis2d.f90:488
subroutine get_polyverts(this, ic, polyverts, closed)
Get a 2D array of polygon vertices, listed in.
Definition: Dis2d.f90:1076
subroutine dis3d_df(this)
Define the discretization.
Definition: Dis2d.f90:126
subroutine, public line_unit_vector(x0, y0, z0, x1, y1, z1, xcomp, ycomp, zcomp, vmag)
Calculate the vector components (xcomp, ycomp, and zcomp) for a line defined by two points,...
Definition: DisvGeom.f90:508
integer(i4b) function, public get_node(ilay, irow, icol, nlay, nrow, ncol)
Get node number, given layer, row, and column indices for a structured grid. If any argument is inval...
Definition: GeomUtil.f90:81
subroutine, public get_ijk(nodenumber, nrow, ncol, nlay, irow, icol, ilay)
Get row, column and layer indices from node number and grid dimensions. If nodenumber is invalid,...
Definition: GeomUtil.f90:98
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...
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public ulaprufw(ncol, nrow, kstp, kper, ilay, iout, buf, text, userfmt, nvalues, nwidth, editdesc)
Print 1 layer array with user formatting in wrap format.
subroutine, public ubdsv06(kstp, kper, text, modelnam1, paknam1, modelnam2, paknam2, ibdchn, naux, auxtxt, ncol, nrow, nlay, nlist, iout, delt, pertim, totim)
Write header records for cell-by-cell flow terms for one component of flow.
subroutine, public ulasav(buf, text, kstp, kper, pertim, totim, ncol, nrow, ilay, ichn)
Save 1 layer array on disk.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
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
subroutine, public memorylist_remove(component, subcomponent, context)
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
Definition: Sim.f90:168
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
character(len=linelength) idm_context
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
Structured grid discretization.
Definition: Dis2d.f90:23
Simplifies tracking parameters sourced from the input context.
Definition: Dis2d.f90:73