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