MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
dismodule Module Reference

Data Types

type  distype
 Structured grid discretization. More...
 
type  disfoundtype
 Simplifies tracking parameters sourced from the input context. More...
 

Functions/Subroutines

subroutine, public dis_cr (dis, name_model, input_mempath, inunit, iout)
 Create a new structured discretization object. More...
 
subroutine dis3d_df (this)
 Define the discretization. More...
 
subroutine dis3d_da (this)
 Deallocate variables. More...
 
subroutine source_options (this)
 Copy options from IDM into package. More...
 
subroutine log_options (this, found)
 Write user options to list file. More...
 
subroutine source_dimensions (this)
 Copy dimensions from IDM into package. More...
 
subroutine log_dimensions (this, found)
 Write dimensions to list file. More...
 
subroutine source_griddata (this)
 Copy grid data from IDM into package. More...
 
subroutine log_griddata (this, found)
 Write dimensions to list file. More...
 
subroutine grid_finalize (this)
 Finalize grid (check properties, allocate arrays, compute connections) More...
 
subroutine write_grb (this, icelltype)
 Write a binary grid file. More...
 
subroutine nodeu_to_string (this, nodeu, str)
 Convert a user nodenumber to a string (nodenumber) or (k,i,j) More...
 
subroutine nodeu_to_array (this, nodeu, arr)
 Convert a user nodenumber to an array (nodenumber) or (k,i,j) More...
 
integer(i4b) function get_nodenumber_idx1 (this, nodeu, icheck)
 Get reduced node number from user node number. More...
 
integer(i4b) function get_nodenumber_idx3 (this, k, i, j, icheck)
 Get reduced node number from layer, row and column indices. More...
 
subroutine allocate_scalars (this, name_model, input_mempath)
 Allocate and initialize scalar variables. More...
 
subroutine allocate_arrays (this)
 Allocate and initialize arrays. More...
 
integer(i4b) function nodeu_from_string (this, lloc, istart, istop, in, iout, line, flag_string, allow_zero)
 Convert a string to a user nodenumber. More...
 
integer(i4b) function nodeu_from_cellid (this, cellid, inunit, iout, flag_string, allow_zero)
 Convert a cellid string to a user nodenumber. More...
 
logical function supports_layers (this)
 Indicates whether the grid discretization supports layers. More...
 
integer(i4b) function get_ncpl (this)
 Return number of cells per layer (nrow * ncol) More...
 
subroutine connection_normal (this, noden, nodem, ihc, xcomp, ycomp, zcomp, ipos)
 Get normal vector components between the cell and a given neighbor. More...
 
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. More...
 
subroutine get_dis_type (this, dis_type)
 Get the discretization type. More...
 
integer(i4b) function get_dis_enum (this)
 Get the discretization type enumeration. More...
 
subroutine get_polyverts (this, ic, polyverts, closed)
 Get a 2D array of polygon vertices, listed in. More...
 
integer(i4b) function get_npolyverts (this, ic, closed)
 Get the number of cell polygon vertices. More...
 
integer(i4b) function get_max_npolyverts (this, closed)
 Get the maximum number of cell polygon vertices. More...
 
subroutine read_int_array (this, line, lloc, istart, istop, iout, in, iarray, aname)
 Read an integer array. More...
 
subroutine read_dbl_array (this, line, lloc, istart, istop, iout, in, darray, aname)
 Read a double precision array. More...
 
subroutine read_layer_array (this, nodelist, darray, ncolbnd, maxbnd, icolbnd, aname, inunit, iout)
 Read a 2d double array into col icolbnd of darray. More...
 
subroutine record_array (this, darray, iout, iprint, idataun, aname, cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
 Record a double precision array. More...
 
subroutine record_srcdst_list_header (this, text, textmodel, textpackage, dstmodel, dstpackage, naux, auxtxt, ibdchn, nlist, iout)
 Record list header for imeth=6. More...
 
subroutine nlarray_to_nodelist (this, darray, nodelist, maxbnd, nbound, aname)
 Convert an integer array (layer numbers) to nodelist. More...
 

Function/Subroutine Documentation

◆ allocate_arrays()

subroutine dismodule::allocate_arrays ( class(distype this)
private

Definition at line 842 of file Dis.f90.

843  ! -- dummy
844  class(DisType) :: this
845  !
846  ! -- Allocate arrays in DisBaseType (mshape, top, bot, area)
847  call this%DisBaseType%allocate_arrays()
848  !
849  ! -- Allocate arrays for DisType
850  if (this%nodes < this%nodesuser) then
851  call mem_allocate(this%nodeuser, this%nodes, 'NODEUSER', this%memoryPath)
852  call mem_allocate(this%nodereduced, this%nodesuser, 'NODEREDUCED', &
853  this%memoryPath)
854  else
855  call mem_allocate(this%nodeuser, 1, 'NODEUSER', this%memoryPath)
856  call mem_allocate(this%nodereduced, 1, 'NODEREDUCED', this%memoryPath)
857  end if
858  !
859  ! -- Initialize
860  this%mshape(1) = this%nlay
861  this%mshape(2) = this%nrow
862  this%mshape(3) = this%ncol
863  !

◆ allocate_scalars()

subroutine dismodule::allocate_scalars ( class(distype this,
character(len=*), intent(in)  name_model,
character(len=*), intent(in)  input_mempath 
)
private

Definition at line 818 of file Dis.f90.

819  ! -- dummy
820  class(DisType) :: this
821  character(len=*), intent(in) :: name_model
822  character(len=*), intent(in) :: input_mempath
823  !
824  ! -- Allocate parent scalars
825  call this%DisBaseType%allocate_scalars(name_model, input_mempath)
826  !
827  ! -- Allocate
828  call mem_allocate(this%nlay, 'NLAY', this%memoryPath)
829  call mem_allocate(this%nrow, 'NROW', this%memoryPath)
830  call mem_allocate(this%ncol, 'NCOL', this%memoryPath)
831  !
832  ! -- Initialize
833  this%nlay = 0
834  this%nrow = 0
835  this%ncol = 0
836  this%ndim = 3
837  !

◆ connection_normal()

subroutine dismodule::connection_normal ( class(distype this,
integer(i4b), intent(in)  noden,
integer(i4b), intent(in)  nodem,
integer(i4b), intent(in)  ihc,
real(dp), intent(inout)  xcomp,
real(dp), intent(inout)  ycomp,
real(dp), intent(inout)  zcomp,
integer(i4b), intent(in)  ipos 
)
private

The normal points outward from the shared face between noden and nodem.

Parameters
[in]nodencell (reduced nn)
[in]nodemneighbor (reduced nn)
[in]ihchorizontal connection flag

Definition at line 1068 of file Dis.f90.

1070  ! -- dummy
1071  class(DisType) :: this
1072  integer(I4B), intent(in) :: noden !< cell (reduced nn)
1073  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
1074  integer(I4B), intent(in) :: ihc !< horizontal connection flag
1075  real(DP), intent(inout) :: xcomp
1076  real(DP), intent(inout) :: ycomp
1077  real(DP), intent(inout) :: zcomp
1078  integer(I4B), intent(in) :: ipos
1079  ! -- local
1080  integer(I4B) :: nodeu1, i1, j1, k1
1081  integer(I4B) :: nodeu2, i2, j2, k2
1082  !
1083  ! -- Set vector components based on ihc
1084  if (ihc == 0) then
1085  xcomp = dzero
1086  ycomp = dzero
1087  if (nodem < noden) then
1088  !
1089  ! -- nodem must be above noden, so upward connection
1090  zcomp = done
1091  else
1092  !
1093  ! -- nodem must be below noden, so downward connection
1094  zcomp = -done
1095  end if
1096  else
1097  xcomp = dzero
1098  ycomp = dzero
1099  zcomp = dzero
1100  nodeu1 = this%get_nodeuser(noden)
1101  nodeu2 = this%get_nodeuser(nodem)
1102  call get_ijk(nodeu1, this%nrow, this%ncol, this%nlay, i1, j1, k1)
1103  call get_ijk(nodeu2, this%nrow, this%ncol, this%nlay, i2, j2, k2)
1104  if (i2 < i1) then ! back
1105  ycomp = done
1106  elseif (j2 < j1) then ! left
1107  xcomp = -done
1108  elseif (j2 > j1) then ! right
1109  xcomp = done
1110  else ! front
1111  ycomp = -done
1112  end if
1113  !
1114  end if
1115  !
Here is the call graph for this function:

◆ connection_vector()

subroutine dismodule::connection_vector ( class(distype this,
integer(i4b), intent(in)  noden,
integer(i4b), intent(in)  nodem,
logical, intent(in)  nozee,
real(dp), intent(in)  satn,
real(dp), intent(in)  satm,
integer(i4b), intent(in)  ihc,
real(dp), intent(inout)  xcomp,
real(dp), intent(inout)  ycomp,
real(dp), intent(inout)  zcomp,
real(dp), intent(inout)  conlen 
)
private

Saturation must be provided to compute cell center vertical coordinates. Also return the straight-line connection length.

Parameters
[in]nodencell (reduced nn)
[in]nodemneighbor (reduced nn)
[in]ihchorizontal connection flag

Definition at line 1123 of file Dis.f90.

1125  ! -- modules
1126  use disvgeom, only: line_unit_vector
1127  ! -- dummy
1128  class(DisType) :: this
1129  integer(I4B), intent(in) :: noden !< cell (reduced nn)
1130  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
1131  logical, intent(in) :: nozee
1132  real(DP), intent(in) :: satn
1133  real(DP), intent(in) :: satm
1134  integer(I4B), intent(in) :: ihc !< horizontal connection flag
1135  real(DP), intent(inout) :: xcomp
1136  real(DP), intent(inout) :: ycomp
1137  real(DP), intent(inout) :: zcomp
1138  real(DP), intent(inout) :: conlen
1139  ! -- local
1140  real(DP) :: z1, z2
1141  real(DP) :: x1, y1, x2, y2
1142  real(DP) :: ds
1143  integer(I4B) :: i1, i2, j1, j2, k1, k2
1144  integer(I4B) :: nodeu1, nodeu2, ipos
1145  !
1146  ! -- Set vector components based on ihc
1147  if (ihc == 0) then
1148  !
1149  ! -- vertical connection; set zcomp positive upward
1150  xcomp = dzero
1151  ycomp = dzero
1152  if (nodem < noden) then
1153  zcomp = done
1154  else
1155  zcomp = -done
1156  end if
1157  z1 = this%bot(noden) + dhalf * (this%top(noden) - this%bot(noden))
1158  z2 = this%bot(nodem) + dhalf * (this%top(nodem) - this%bot(nodem))
1159  conlen = abs(z2 - z1)
1160  else
1161  !
1162  if (nozee) then
1163  z1 = dzero
1164  z2 = dzero
1165  else
1166  z1 = this%bot(noden) + dhalf * satn * (this%top(noden) - this%bot(noden))
1167  z2 = this%bot(nodem) + dhalf * satm * (this%top(nodem) - this%bot(nodem))
1168  end if
1169  ipos = this%con%getjaindex(noden, nodem)
1170  ds = this%con%cl1(this%con%jas(ipos)) + this%con%cl2(this%con%jas(ipos))
1171  nodeu1 = this%get_nodeuser(noden)
1172  nodeu2 = this%get_nodeuser(nodem)
1173  call get_ijk(nodeu1, this%nrow, this%ncol, this%nlay, i1, j1, k1)
1174  call get_ijk(nodeu2, this%nrow, this%ncol, this%nlay, i2, j2, k2)
1175  x1 = dzero
1176  x2 = dzero
1177  y1 = dzero
1178  y2 = dzero
1179  if (i2 < i1) then ! back
1180  y2 = ds
1181  elseif (j2 < j1) then ! left
1182  x2 = -ds
1183  elseif (j2 > j1) then ! right
1184  x2 = ds
1185  else ! front
1186  y2 = -ds
1187  end if
1188  call line_unit_vector(x1, y1, z1, x2, y2, z2, xcomp, ycomp, zcomp, conlen)
1189  end if
1190  !
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:475
Here is the call graph for this function:

◆ dis3d_da()

subroutine dismodule::dis3d_da ( class(distype this)
private

Definition at line 155 of file Dis.f90.

156  ! -- dummy
157  class(DisType) :: this
158  !
159  ! -- Deallocate idm memory
160  call memorystore_remove(this%name_model, 'DIS', idm_context)
161  !
162  ! -- DisBaseType deallocate
163  call this%DisBaseType%dis_da()
164  !
165  ! -- Deallocate scalars
166  call mem_deallocate(this%nlay)
167  call mem_deallocate(this%nrow)
168  call mem_deallocate(this%ncol)
169  call mem_deallocate(this%delr)
170  call mem_deallocate(this%delc)
171  call mem_deallocate(this%cellx)
172  call mem_deallocate(this%celly)
173  !
174  ! -- Deallocate Arrays
175  call mem_deallocate(this%nodereduced)
176  call mem_deallocate(this%nodeuser)
177  call mem_deallocate(this%top2d)
178  call mem_deallocate(this%bot3d)
179  call mem_deallocate(this%idomain)
180  !
Here is the call graph for this function:

◆ dis3d_df()

subroutine dismodule::dis3d_df ( class(distype this)
private

Definition at line 131 of file Dis.f90.

132  ! -- dummy
133  class(DisType) :: this
134  !
135  ! -- Transfer the data from the memory manager into this package object
136  if (this%inunit /= 0) then
137  !
138  ! -- source input options
139  call this%source_options()
140  !
141  ! -- source input dimensions
142  call this%source_dimensions()
143  !
144  ! -- source input griddata
145  call this%source_griddata()
146  end if
147  !
148  ! -- Final grid initialization
149  call this%grid_finalize()
150  !

◆ dis_cr()

subroutine, public dismodule::dis_cr ( class(disbasetype), pointer  dis,
character(len=*), intent(in)  name_model,
character(len=*), intent(in)  input_mempath,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout 
)

Definition at line 98 of file Dis.f90.

99  ! -- dummy
100  class(DisBaseType), pointer :: dis
101  character(len=*), intent(in) :: name_model
102  character(len=*), intent(in) :: input_mempath
103  integer(I4B), intent(in) :: inunit
104  integer(I4B), intent(in) :: iout
105  ! -- locals
106  type(DisType), pointer :: disnew
107  ! -- formats
108  character(len=*), parameter :: fmtheader = &
109  "(1X, /1X, 'DIS -- STRUCTURED GRID DISCRETIZATION PACKAGE,', &
110  &' VERSION 2 : 3/27/2014 - INPUT READ FROM MEMPATH: ', A, /)"
111  !
112  allocate (disnew)
113  dis => disnew
114  call disnew%allocate_scalars(name_model, input_mempath)
115  dis%inunit = inunit
116  dis%iout = iout
117  !
118  ! -- If dis enabled
119  if (inunit > 0) then
120  !
121  ! -- Identify package
122  if (iout > 0) then
123  write (iout, fmtheader) dis%input_mempath
124  end if
125  end if
126  !
Here is the caller graph for this function:

◆ get_dis_enum()

integer(i4b) function dismodule::get_dis_enum ( class(distype), intent(in)  this)
private

Definition at line 1205 of file Dis.f90.

1206  use constantsmodule, only: dis
1207  class(DisType), intent(in) :: this
1208  integer(I4B) :: dis_enum
1209  dis_enum = dis
This module contains simulation constants.
Definition: Constants.f90:9
@ dis
DIS6 discretization.
Definition: Constants.f90:155

◆ get_dis_type()

subroutine dismodule::get_dis_type ( class(distype), intent(in)  this,
character(len=*), intent(out)  dis_type 
)

Definition at line 1195 of file Dis.f90.

1196  ! -- dummy
1197  class(DisType), intent(in) :: this
1198  character(len=*), intent(out) :: dis_type
1199  !
1200  dis_type = "DIS"
1201  !

◆ get_max_npolyverts()

integer(i4b) function dismodule::get_max_npolyverts ( class(distype), intent(inout)  this,
logical(lgp), intent(in), optional  closed 
)
private
Parameters
[in]closedwhether to close the polygon, duplicating a vertex

Definition at line 1274 of file Dis.f90.

1275  class(DisType), intent(inout) :: this
1276  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
1277  integer(I4B) :: max_npolyverts
1278  max_npolyverts = 4
1279  if (present(closed)) then
1280  if (closed) max_npolyverts = 5
1281  end if

◆ get_ncpl()

integer(i4b) function dismodule::get_ncpl ( class(distype this)
private

Definition at line 1058 of file Dis.f90.

1059  integer(I4B) :: get_ncpl
1060  class(DisType) :: this
1061  get_ncpl = this%nrow * this%ncol

◆ get_nodenumber_idx1()

integer(i4b) function dismodule::get_nodenumber_idx1 ( class(distype), intent(in)  this,
integer(i4b), intent(in)  nodeu,
integer(i4b), intent(in)  icheck 
)
private

Definition at line 743 of file Dis.f90.

744  ! -- return
745  integer(I4B) :: nodenumber
746  ! -- dummy
747  class(DisType), intent(in) :: this
748  integer(I4B), intent(in) :: nodeu
749  integer(I4B), intent(in) :: icheck
750  !
751  ! -- check the node number if requested
752  if (icheck /= 0) then
753  !
754  ! -- If within valid range, convert to reduced nodenumber
755  if (nodeu < 1 .or. nodeu > this%nodesuser) then
756  write (errmsg, '(a,i0,a)') &
757  'Node number (', nodeu, &
758  ') less than 1 or greater than the number of nodes.'
759  call store_error(errmsg)
760  nodenumber = 0
761  else
762  nodenumber = nodeu
763  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
764  end if
765  else
766  nodenumber = nodeu
767  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
768  end if
769  !
Here is the call graph for this function:

◆ get_nodenumber_idx3()

integer(i4b) function dismodule::get_nodenumber_idx3 ( class(distype), intent(in)  this,
integer(i4b), intent(in)  k,
integer(i4b), intent(in)  i,
integer(i4b), intent(in)  j,
integer(i4b), intent(in)  icheck 
)
private

Definition at line 774 of file Dis.f90.

775  ! -- return
776  integer(I4B) :: nodenumber
777  ! -- dummy
778  class(DisType), intent(in) :: this
779  integer(I4B), intent(in) :: k, i, j
780  integer(I4B), intent(in) :: icheck
781  ! -- local
782  integer(I4B) :: nodeu
783  ! formats
784  character(len=*), parameter :: fmterr = &
785  "('Error in structured-grid cell indices: layer = ',i0,', &
786  &row = ',i0,', column = ',i0)"
787  !
788  nodeu = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
789  if (nodeu < 1) then
790  write (errmsg, fmterr) k, i, j
791  call store_error(errmsg, terminate=.true.)
792  end if
793  nodenumber = nodeu
794  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
795  !
796  ! -- check the node number if requested
797  if (icheck /= 0) then
798  !
799  if (k < 1 .or. k > this%nlay) &
800  call store_error('Layer less than one or greater than nlay')
801  if (i < 1 .or. i > this%nrow) &
802  call store_error('Row less than one or greater than nrow')
803  if (j < 1 .or. j > this%ncol) &
804  call store_error('Column less than one or greater than ncol')
805  !
806  ! -- Error if outside of range
807  if (nodeu < 1 .or. nodeu > this%nodesuser) then
808  write (errmsg, '(a,i0,a)') &
809  'Node number (', nodeu, ')less than 1 or greater than nodes.'
810  call store_error(errmsg)
811  end if
812  end if
813  !
Here is the call graph for this function:

◆ get_npolyverts()

integer(i4b) function dismodule::get_npolyverts ( class(distype), intent(inout)  this,
integer(i4b), intent(in)  ic,
logical(lgp), intent(in), optional  closed 
)
private
Parameters
[in]iccell number (reduced)
[in]closedwhether to close the polygon, duplicating a vertex

Definition at line 1262 of file Dis.f90.

1263  class(DisType), intent(inout) :: this
1264  integer(I4B), intent(in) :: ic !< cell number (reduced)
1265  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
1266  integer(I4B) :: npolyverts
1267  npolyverts = 4
1268  if (present(closed)) then
1269  if (closed) npolyverts = 5
1270  end if

◆ get_polyverts()

subroutine dismodule::get_polyverts ( class(distype), intent(inout)  this,
integer(i4b), intent(in)  ic,
real(dp), dimension(:, :), intent(out), allocatable  polyverts,
logical(lgp), intent(in), optional  closed 
)

clockwise order beginning with the lower left corner

Parameters
[in]iccell number (reduced)
[out]polyvertspolygon vertices (column-major indexing)
[in]closedwhether to close the polygon, duplicating a vertex

Definition at line 1216 of file Dis.f90.

1217  ! -- dummy
1218  class(DisType), intent(inout) :: this
1219  integer(I4B), intent(in) :: ic !< cell number (reduced)
1220  real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices (column-major indexing)
1221  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
1222  ! -- local
1223  integer(I4B) :: icu, nverts, irow, jcol, klay
1224  real(DP) :: cellx, celly, dxhalf, dyhalf
1225  logical(LGP) :: lclosed
1226  !
1227  nverts = 4
1228  !
1229  ! check closed option
1230  if (.not. (present(closed))) then
1231  lclosed = .false.
1232  else
1233  lclosed = closed
1234  end if
1235  !
1236  ! allocate vertices array
1237  if (lclosed) then
1238  allocate (polyverts(2, nverts + 1))
1239  else
1240  allocate (polyverts(2, nverts))
1241  end if
1242  !
1243  ! set vertices
1244  icu = this%get_nodeuser(ic)
1245  call get_ijk(icu, this%nrow, this%ncol, this%nlay, irow, jcol, klay)
1246  cellx = this%cellx(jcol)
1247  celly = this%celly(irow)
1248  dxhalf = dhalf * this%delr(jcol)
1249  dyhalf = dhalf * this%delc(irow)
1250  polyverts(:, 1) = (/cellx - dxhalf, celly - dyhalf/) ! SW
1251  polyverts(:, 2) = (/cellx - dxhalf, celly + dyhalf/) ! NW
1252  polyverts(:, 3) = (/cellx + dxhalf, celly + dyhalf/) ! NE
1253  polyverts(:, 4) = (/cellx + dxhalf, celly - dyhalf/) ! SE
1254  !
1255  ! close if enabled
1256  if (lclosed) &
1257  polyverts(:, nverts + 1) = polyverts(:, 1)
1258  !
Here is the call graph for this function:

◆ grid_finalize()

subroutine dismodule::grid_finalize ( class(distype this)
private

Definition at line 385 of file Dis.f90.

386  ! -- modules
388  ! -- dummy
389  class(DisType) :: this
390  ! -- locals
391  integer(I4B) :: n, i, j, k
392  integer(I4B) :: node
393  integer(I4B) :: noder
394  integer(I4B) :: nrsize
395  real(DP) :: top
396  real(DP) :: dz
397  ! -- formats
398  character(len=*), parameter :: fmtdz = &
399  "('CELL (',i0,',',i0,',',i0,') THICKNESS <= 0. ', &
400  &'TOP, BOT: ',2(1pg24.15))"
401  character(len=*), parameter :: fmtnr = &
402  "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',&
403  &/1x, 'Number of user nodes: ',I0,&
404  &/1X, 'Number of nodes in solution: ', I0, //)"
405  !
406  ! -- count active cells
407  this%nodes = 0
408  do k = 1, this%nlay
409  do i = 1, this%nrow
410  do j = 1, this%ncol
411  if (this%idomain(j, i, k) > 0) this%nodes = this%nodes + 1
412  end do
413  end do
414  end do
415  !
416  ! -- Check to make sure nodes is a valid number
417  if (this%nodes == 0) then
418  call store_error('Model does not have any active nodes. &
419  &Ensure IDOMAIN array has some values greater &
420  &than zero.')
421  call store_error_filename(this%input_fname)
422  end if
423  !
424  ! -- Check cell thicknesses
425  n = 0
426  do k = 1, this%nlay
427  do i = 1, this%nrow
428  do j = 1, this%ncol
429  if (this%idomain(j, i, k) < 1) cycle
430  if (k > 1) then
431  top = this%bot3d(j, i, k - 1)
432  else
433  top = this%top2d(j, i)
434  end if
435  dz = top - this%bot3d(j, i, k)
436  if (dz <= dzero) then
437  n = n + 1
438  write (errmsg, fmt=fmtdz) k, i, j, top, this%bot3d(j, i, k)
439  call store_error(errmsg)
440  end if
441  end do
442  end do
443  end do
444  if (count_errors() > 0) then
445  call store_error_filename(this%input_fname)
446  end if
447  !
448  ! -- Write message if reduced grid
449  if (this%nodes < this%nodesuser) then
450  write (this%iout, fmtnr) this%nodesuser, this%nodes
451  end if
452  !
453  ! -- Array size is now known, so allocate
454  call this%allocate_arrays()
455  !
456  ! -- Fill the nodereduced array with the reduced nodenumber, or
457  ! a negative number to indicate it is a pass-through cell, or
458  ! a zero to indicate that the cell is excluded from the
459  ! solution.
460  if (this%nodes < this%nodesuser) then
461  node = 1
462  noder = 1
463  do k = 1, this%nlay
464  do i = 1, this%nrow
465  do j = 1, this%ncol
466  if (this%idomain(j, i, k) > 0) then
467  this%nodereduced(node) = noder
468  noder = noder + 1
469  elseif (this%idomain(j, i, k) < 0) then
470  this%nodereduced(node) = -1
471  else
472  this%nodereduced(node) = 0
473  end if
474  node = node + 1
475  end do
476  end do
477  end do
478  end if
479  !
480  ! -- allocate and fill nodeuser if a reduced grid
481  if (this%nodes < this%nodesuser) then
482  node = 1
483  noder = 1
484  do k = 1, this%nlay
485  do i = 1, this%nrow
486  do j = 1, this%ncol
487  if (this%idomain(j, i, k) > 0) then
488  this%nodeuser(noder) = node
489  noder = noder + 1
490  end if
491  node = node + 1
492  end do
493  end do
494  end do
495  end if
496  !
497  ! -- fill x,y coordinate arrays
498  this%cellx(1) = dhalf * this%delr(1)
499  this%celly(this%nrow) = dhalf * this%delc(this%nrow)
500  do j = 2, this%ncol
501  this%cellx(j) = this%cellx(j - 1) + dhalf * this%delr(j - 1) + &
502  dhalf * this%delr(j)
503  end do
504  ! -- row number increases in negative y direction:
505  do i = this%nrow - 1, 1, -1
506  this%celly(i) = this%celly(i + 1) + dhalf * this%delc(i + 1) + &
507  dhalf * this%delc(i)
508  end do
509  !
510  ! -- Move top2d and botm3d into top and bot, and calculate area
511  node = 0
512  do k = 1, this%nlay
513  do i = 1, this%nrow
514  do j = 1, this%ncol
515  node = node + 1
516  noder = node
517  if (this%nodes < this%nodesuser) noder = this%nodereduced(node)
518  if (noder <= 0) cycle
519  if (k > 1) then
520  top = this%bot3d(j, i, k - 1)
521  else
522  top = this%top2d(j, i)
523  end if
524  this%top(noder) = top
525  this%bot(noder) = this%bot3d(j, i, k)
526  this%area(noder) = this%delr(j) * this%delc(i)
527  this%xc(noder) = this%cellx(j)
528  this%yc(noder) = this%celly(i)
529  end do
530  end do
531  end do
532  !
533  ! -- create and fill the connections object
534  nrsize = 0
535  if (this%nodes < this%nodesuser) nrsize = this%nodes
536  allocate (this%con)
537  call this%con%disconnections(this%name_model, this%nodes, &
538  this%ncol, this%nrow, this%nlay, &
539  nrsize, this%delr, this%delc, &
540  this%top, this%bot, this%nodereduced, &
541  this%nodeuser)
542  this%nja = this%con%nja
543  this%njas = this%con%njas
544  !
Here is the call graph for this function:

◆ log_dimensions()

subroutine dismodule::log_dimensions ( class(distype this,
type(disfoundtype), intent(in)  found 
)
private

Definition at line 306 of file Dis.f90.

307  ! -- dummy
308  class(DisType) :: this
309  type(DisFoundType), intent(in) :: found
310  !
311  write (this%iout, '(1x,a)') 'Setting Discretization Dimensions'
312  !
313  if (found%nlay) then
314  write (this%iout, '(4x,a,i0)') 'NLAY = ', this%nlay
315  end if
316  !
317  if (found%nrow) then
318  write (this%iout, '(4x,a,i0)') 'NROW = ', this%nrow
319  end if
320  !
321  if (found%ncol) then
322  write (this%iout, '(4x,a,i0)') 'NCOL = ', this%ncol
323  end if
324  !
325  write (this%iout, '(1x,a,/)') 'End Setting Discretization Dimensions'
326  !

◆ log_griddata()

subroutine dismodule::log_griddata ( class(distype this,
type(disfoundtype), intent(in)  found 
)
private

Definition at line 352 of file Dis.f90.

353  ! -- dummy
354  class(DisType) :: this
355  type(DisFoundType), intent(in) :: found
356  !
357  write (this%iout, '(1x,a)') 'Setting Discretization Griddata'
358  !
359  if (found%delr) then
360  write (this%iout, '(4x,a)') 'DELR set from input file'
361  end if
362  !
363  if (found%delc) then
364  write (this%iout, '(4x,a)') 'DELC set from input file'
365  end if
366  !
367  if (found%top) then
368  write (this%iout, '(4x,a)') 'TOP set from input file'
369  end if
370  !
371  if (found%botm) then
372  write (this%iout, '(4x,a)') 'BOTM set from input file'
373  end if
374  !
375  if (found%idomain) then
376  write (this%iout, '(4x,a)') 'IDOMAIN set from input file'
377  end if
378  !
379  write (this%iout, '(1x,a,/)') 'End Setting Discretization Griddata'
380  !

◆ log_options()

subroutine dismodule::log_options ( class(distype this,
type(disfoundtype), intent(in)  found 
)
private

Definition at line 210 of file Dis.f90.

211  ! -- dummy
212  class(DisType) :: this
213  type(DisFoundType), intent(in) :: found
214  !
215  write (this%iout, '(1x,a)') 'Setting Discretization Options'
216  !
217  if (found%length_units) then
218  write (this%iout, '(4x,a,i0)') 'Model length unit [0=UND, 1=FEET, &
219  &2=METERS, 3=CENTIMETERS] set as ', this%lenuni
220  end if
221  !
222  if (found%nogrb) then
223  write (this%iout, '(4x,a,i0)') 'Binary grid file [0=GRB, 1=NOGRB] &
224  &set as ', this%nogrb
225  end if
226  !
227  if (found%xorigin) then
228  write (this%iout, '(4x,a,G0)') 'XORIGIN = ', this%xorigin
229  end if
230  !
231  if (found%yorigin) then
232  write (this%iout, '(4x,a,G0)') 'YORIGIN = ', this%yorigin
233  end if
234  !
235  if (found%angrot) then
236  write (this%iout, '(4x,a,G0)') 'ANGROT = ', this%angrot
237  end if
238  !
239  write (this%iout, '(1x,a,/)') 'End Setting Discretization Options'
240  !

◆ nlarray_to_nodelist()

subroutine dismodule::nlarray_to_nodelist ( class(distype this,
integer(i4b), dimension(:), pointer, contiguous  darray,
integer(i4b), dimension(maxbnd), intent(inout)  nodelist,
integer(i4b), intent(in)  maxbnd,
integer(i4b), intent(inout)  nbound,
character(len=*), intent(in)  aname 
)
private

Definition at line 1572 of file Dis.f90.

1573  ! -- dummy
1574  class(DisType) :: this
1575  integer(I4B), intent(in) :: maxbnd
1576  integer(I4B), dimension(:), pointer, contiguous :: darray
1577  integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
1578  integer(I4B), intent(inout) :: nbound
1579  character(len=*), intent(in) :: aname
1580  ! -- local
1581  integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr
1582  !
1583  ! -- set variables
1584  nlay = this%mshape(1)
1585  nrow = this%mshape(2)
1586  ncol = this%mshape(3)
1587  !
1588  if (this%ndim > 1) then
1589  !
1590  nval = ncol * nrow
1591  !
1592  ! -- Copy array into nodelist
1593  ipos = 1
1594  ierr = 0
1595  do ir = 1, nrow
1596  do ic = 1, ncol
1597  nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
1598  il = darray(nodeu)
1599  if (il < 1 .or. il > nlay) then
1600  write (errmsg, '(a,1x,i0)') 'Invalid layer number:', il
1601  call store_error(errmsg, terminate=.true.)
1602  end if
1603  nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
1604  noder = this%get_nodenumber(nodeu, 0)
1605  if (ipos > maxbnd) then
1606  ierr = ipos
1607  else
1608  nodelist(ipos) = noder
1609  end if
1610  ipos = ipos + 1
1611  end do
1612  end do
1613  !
1614  ! -- Check for errors
1615  nbound = ipos - 1
1616  if (ierr > 0) then
1617  write (errmsg, '(a,1x,i0)') &
1618  'MAXBOUND dimension is too small.'// &
1619  'INCREASE MAXBOUND TO:', ierr
1620  call store_error(errmsg, terminate=.true.)
1621  end if
1622  !
1623  ! -- If nbound < maxbnd, then initialize nodelist to zero in this range
1624  if (nbound < maxbnd) then
1625  do ipos = nbound + 1, maxbnd
1626  nodelist(ipos) = 0
1627  end do
1628  end if
1629  !
1630  else
1631  !
1632  ! -- For unstructured, read nodelist directly, then check node numbers
1633  nodelist = darray
1634  do noder = 1, maxbnd
1635  if (noder < 1 .or. noder > this%nodes) then
1636  write (errmsg, '(a,1x,i0)') 'Invalid node number:', noder
1637  call store_error(errmsg, terminate=.true.)
1638  end if
1639  end do
1640  nbound = maxbnd
1641  !
1642  end if
1643  !
Here is the call graph for this function:

◆ nodeu_from_cellid()

integer(i4b) function dismodule::nodeu_from_cellid ( class(distype this,
character(len=*), intent(inout)  cellid,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout,
logical, intent(in), optional  flag_string,
logical, intent(in), optional  allow_zero 
)
private

If flag_string is present and true, the first token may be non-numeric (e.g. boundary name). In this case, return -2.

If allow_zero is present and true, and all indices are zero, the result can be zero. If allow_zero is false, a zero in any index is an error.

Definition at line 963 of file Dis.f90.

965  ! -- return
966  integer(I4B) :: nodeu
967  ! -- dummy
968  class(DisType) :: this
969  character(len=*), intent(inout) :: cellid
970  integer(I4B), intent(in) :: inunit
971  integer(I4B), intent(in) :: iout
972  logical, optional, intent(in) :: flag_string
973  logical, optional, intent(in) :: allow_zero
974  ! -- local
975  integer(I4B) :: lloclocal, istart, istop, ndum, n
976  integer(I4B) :: k, i, j, nlay, nrow, ncol
977  integer(I4B) :: istat
978  real(DP) :: r
979  !
980  if (present(flag_string)) then
981  if (flag_string) then
982  ! Check to see if first token in cellid can be read as an integer.
983  lloclocal = 1
984  call urword(cellid, lloclocal, istart, istop, 1, ndum, r, iout, inunit)
985  read (cellid(istart:istop), *, iostat=istat) n
986  if (istat /= 0) then
987  ! First token in cellid is not an integer; return flag to this effect.
988  nodeu = -2
989  return
990  end if
991  end if
992  end if
993  !
994  nlay = this%mshape(1)
995  nrow = this%mshape(2)
996  ncol = this%mshape(3)
997  !
998  lloclocal = 1
999  call urword(cellid, lloclocal, istart, istop, 2, k, r, iout, inunit)
1000  call urword(cellid, lloclocal, istart, istop, 2, i, r, iout, inunit)
1001  call urword(cellid, lloclocal, istart, istop, 2, j, r, iout, inunit)
1002  !
1003  if (k == 0 .and. i == 0 .and. j == 0) then
1004  if (present(allow_zero)) then
1005  if (allow_zero) then
1006  nodeu = 0
1007  return
1008  end if
1009  end if
1010  end if
1011  !
1012  errmsg = ""
1013  !
1014  if (k < 1 .or. k > nlay) then
1015  write (errmsg, '(a,i0,a)') &
1016  'Layer number in list (', k, ') is outside of the grid.'
1017  end if
1018  if (i < 1 .or. i > nrow) then
1019  write (errmsg, '(a,1x,a,i0,a)') &
1020  trim(adjustl(errmsg)), 'Row number in list (', i, &
1021  ') is outside of the grid.'
1022  end if
1023  if (j < 1 .or. j > ncol) then
1024  write (errmsg, '(a,1x,a,i0,a)') &
1025  trim(adjustl(errmsg)), 'Column number in list (', j, &
1026  ') is outside of the grid.'
1027  end if
1028  !
1029  nodeu = get_node(k, i, j, nlay, nrow, ncol)
1030  !
1031  if (nodeu < 1 .or. nodeu > this%nodesuser) then
1032  write (errmsg, '(a,1x,a,i0,a)') &
1033  trim(adjustl(errmsg)), &
1034  "Cell number cannot be determined for cellid ("// &
1035  trim(adjustl(cellid))//") and results in a user "// &
1036  "node number (", nodeu, ") that is outside of the grid."
1037  end if
1038  !
1039  if (len_trim(adjustl(errmsg)) > 0) then
1040  call store_error(errmsg)
1041  call store_error_unit(inunit)
1042  end if
1043  !
Here is the call graph for this function:

◆ nodeu_from_string()

integer(i4b) function dismodule::nodeu_from_string ( class(distype this,
integer(i4b), intent(inout)  lloc,
integer(i4b), intent(inout)  istart,
integer(i4b), intent(inout)  istop,
integer(i4b), intent(in)  in,
integer(i4b), intent(in)  iout,
character(len=*), intent(inout)  line,
logical, intent(in), optional  flag_string,
logical, intent(in), optional  allow_zero 
)
private

Parse layer, row and column and return user nodenumber. If flag_string is present and true, the first token may be non-numeric (e.g. boundary name). In this case, return -2.

Definition at line 872 of file Dis.f90.

874  ! -- dummy
875  class(DisType) :: this
876  integer(I4B), intent(inout) :: lloc
877  integer(I4B), intent(inout) :: istart
878  integer(I4B), intent(inout) :: istop
879  integer(I4B), intent(in) :: in
880  integer(I4B), intent(in) :: iout
881  character(len=*), intent(inout) :: line
882  logical, optional, intent(in) :: flag_string
883  logical, optional, intent(in) :: allow_zero
884  integer(I4B) :: nodeu
885  ! -- local
886  integer(I4B) :: k, i, j, nlay, nrow, ncol
887  integer(I4B) :: lloclocal, ndum, istat, n
888  real(DP) :: r
889  !
890  if (present(flag_string)) then
891  if (flag_string) then
892  ! Check to see if first token in line can be read as an integer.
893  lloclocal = lloc
894  call urword(line, lloclocal, istart, istop, 1, ndum, r, iout, in)
895  read (line(istart:istop), *, iostat=istat) n
896  if (istat /= 0) then
897  ! First token in line is not an integer; return flag to this effect.
898  nodeu = -2
899  return
900  end if
901  end if
902  end if
903  !
904  nlay = this%mshape(1)
905  nrow = this%mshape(2)
906  ncol = this%mshape(3)
907  !
908  call urword(line, lloc, istart, istop, 2, k, r, iout, in)
909  call urword(line, lloc, istart, istop, 2, i, r, iout, in)
910  call urword(line, lloc, istart, istop, 2, j, r, iout, in)
911  !
912  if (k == 0 .and. i == 0 .and. j == 0) then
913  if (present(allow_zero)) then
914  if (allow_zero) then
915  nodeu = 0
916  return
917  end if
918  end if
919  end if
920  !
921  errmsg = ""
922  !
923  if (k < 1 .or. k > nlay) then
924  write (errmsg, '(a,i0,a)') &
925  'Layer number in list (', k, ') is outside of the grid.'
926  end if
927  if (i < 1 .or. i > nrow) then
928  write (errmsg, '(a,1x,a,i0,a)') &
929  trim(adjustl(errmsg)), 'Row number in list (', i, &
930  ') is outside of the grid.'
931  end if
932  if (j < 1 .or. j > ncol) then
933  write (errmsg, '(a,1x,a,i0,a)') &
934  trim(adjustl(errmsg)), 'Column number in list (', j, &
935  ') is outside of the grid.'
936  end if
937  !
938  nodeu = get_node(k, i, j, nlay, nrow, ncol)
939  !
940  if (nodeu < 1 .or. nodeu > this%nodesuser) then
941  write (errmsg, '(a,1x,a,i0,a)') &
942  trim(adjustl(errmsg)), &
943  "Node number in list (", nodeu, ") is outside of the grid. "// &
944  "Cell number cannot be determined in line '"// &
945  trim(adjustl(line))//"'."
946  end if
947  !
948  if (len_trim(adjustl(errmsg)) > 0) then
949  call store_error(errmsg)
950  call store_error_unit(in)
951  end if
952  !
Here is the call graph for this function:

◆ nodeu_to_array()

subroutine dismodule::nodeu_to_array ( class(distype this,
integer(i4b), intent(in)  nodeu,
integer(i4b), dimension(:), intent(inout)  arr 
)
private

Definition at line 713 of file Dis.f90.

714  ! -- dummy
715  class(DisType) :: this
716  integer(I4B), intent(in) :: nodeu
717  integer(I4B), dimension(:), intent(inout) :: arr
718  ! -- local
719  integer(I4B) :: isize
720  integer(I4B) :: i, j, k
721  !
722  ! -- check the size of arr
723  isize = size(arr)
724  if (isize /= this%ndim) then
725  write (errmsg, '(a,i0,a,i0,a)') &
726  'Program error: nodeu_to_array size of array (', isize, &
727  ') is not equal to the discretization dimension (', this%ndim, ')'
728  call store_error(errmsg, terminate=.true.)
729  end if
730  !
731  ! -- get k, i, j
732  call get_ijk(nodeu, this%nrow, this%ncol, this%nlay, i, j, k)
733  !
734  ! -- fill array
735  arr(1) = k
736  arr(2) = i
737  arr(3) = j
738  !
Here is the call graph for this function:

◆ nodeu_to_string()

subroutine dismodule::nodeu_to_string ( class(distype this,
integer(i4b), intent(in)  nodeu,
character(len=*), intent(inout)  str 
)

Definition at line 692 of file Dis.f90.

693  ! -- dummy
694  class(DisType) :: this
695  integer(I4B), intent(in) :: nodeu
696  character(len=*), intent(inout) :: str
697  ! -- local
698  integer(I4B) :: i, j, k
699  character(len=10) :: kstr, istr, jstr
700  !
701  call get_ijk(nodeu, this%nrow, this%ncol, this%nlay, i, j, k)
702  write (kstr, '(i10)') k
703  write (istr, '(i10)') i
704  write (jstr, '(i10)') j
705  str = '('//trim(adjustl(kstr))//','// &
706  trim(adjustl(istr))//','// &
707  trim(adjustl(jstr))//')'
708  !
Here is the call graph for this function:

◆ read_dbl_array()

subroutine dismodule::read_dbl_array ( class(distype), intent(inout)  this,
character(len=*), intent(inout)  line,
integer(i4b), intent(inout)  lloc,
integer(i4b), intent(inout)  istart,
integer(i4b), intent(inout)  istop,
integer(i4b), intent(in)  iout,
integer(i4b), intent(in)  in,
real(dp), dimension(:), intent(inout), pointer, contiguous  darray,
character(len=*), intent(in)  aname 
)
private

Definition at line 1345 of file Dis.f90.

1347  ! -- dummy
1348  class(DisType), intent(inout) :: this
1349  character(len=*), intent(inout) :: line
1350  integer(I4B), intent(inout) :: lloc
1351  integer(I4B), intent(inout) :: istart
1352  integer(I4B), intent(inout) :: istop
1353  integer(I4B), intent(in) :: in
1354  integer(I4B), intent(in) :: iout
1355  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray
1356  character(len=*), intent(in) :: aname
1357  ! -- local
1358  integer(I4B) :: ival
1359  real(DP) :: rval
1360  integer(I4B) :: nlay
1361  integer(I4B) :: nrow
1362  integer(I4B) :: ncol
1363  integer(I4B) :: nval
1364  real(DP), dimension(:), pointer, contiguous :: dtemp
1365  !
1366  ! -- Point the temporary pointer array, which is passed to the reading
1367  ! subroutine. The temporary array will point to dbuff if it is a
1368  ! reduced structured system, or to darray if it is an unstructured
1369  ! model.
1370  nlay = this%mshape(1)
1371  nrow = this%mshape(2)
1372  ncol = this%mshape(3)
1373  !
1374  if (this%nodes < this%nodesuser) then
1375  nval = this%nodesuser
1376  dtemp => this%dbuff
1377  else
1378  nval = this%nodes
1379  dtemp => darray
1380  end if
1381  !
1382  ! -- Read the array
1383  call urword(line, lloc, istart, istop, 1, ival, rval, iout, in)
1384  if (line(istart:istop) .EQ. 'LAYERED') then
1385  !
1386  ! -- Read structured input
1387  call readarray(in, dtemp, aname, this%ndim, ncol, nrow, nlay, nval, &
1388  iout, 1, nlay)
1389  else
1390  !
1391  ! -- Read unstructured input
1392  call readarray(in, dtemp, aname, this%ndim, nval, iout, 0)
1393  end if
1394  !
1395  ! -- If reduced model, then need to copy from dtemp(=>dbuff) to darray
1396  if (this%nodes < this%nodesuser) then
1397  call this%fill_grid_array(dtemp, darray)
1398  end if
1399  !
Here is the call graph for this function:

◆ read_int_array()

subroutine dismodule::read_int_array ( class(distype), intent(inout)  this,
character(len=*), intent(inout)  line,
integer(i4b), intent(inout)  lloc,
integer(i4b), intent(inout)  istart,
integer(i4b), intent(inout)  istop,
integer(i4b), intent(in)  iout,
integer(i4b), intent(in)  in,
integer(i4b), dimension(:), intent(inout), pointer, contiguous  iarray,
character(len=*), intent(in)  aname 
)
private

Definition at line 1286 of file Dis.f90.

1288  ! -- dummy
1289  class(DisType), intent(inout) :: this
1290  character(len=*), intent(inout) :: line
1291  integer(I4B), intent(inout) :: lloc
1292  integer(I4B), intent(inout) :: istart
1293  integer(I4B), intent(inout) :: istop
1294  integer(I4B), intent(in) :: in
1295  integer(I4B), intent(in) :: iout
1296  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: iarray
1297  character(len=*), intent(in) :: aname
1298  ! -- local
1299  integer(I4B) :: ival
1300  real(DP) :: rval
1301  integer(I4B) :: nlay
1302  integer(I4B) :: nrow
1303  integer(I4B) :: ncol
1304  integer(I4B) :: nval
1305  integer(I4B), dimension(:), pointer, contiguous :: itemp
1306  !
1307  ! -- Point the temporary pointer array, which is passed to the reading
1308  ! subroutine. The temporary array will point to ibuff if it is a
1309  ! reduced structured system, or to iarray if it is an unstructured
1310  ! model.
1311  nlay = this%mshape(1)
1312  nrow = this%mshape(2)
1313  ncol = this%mshape(3)
1314  !
1315  if (this%nodes < this%nodesuser) then
1316  nval = this%nodesuser
1317  itemp => this%ibuff
1318  else
1319  nval = this%nodes
1320  itemp => iarray
1321  end if
1322  !
1323  ! -- Read the array
1324  call urword(line, lloc, istart, istop, 1, ival, rval, iout, in)
1325  if (line(istart:istop) .EQ. 'LAYERED') then
1326  !
1327  ! -- Read layered input
1328  call readarray(in, itemp, aname, this%ndim, ncol, nrow, nlay, nval, &
1329  iout, 1, nlay)
1330  else
1331  !
1332  ! -- Read unstructured input
1333  call readarray(in, itemp, aname, this%ndim, nval, iout, 0)
1334  end if
1335  !
1336  ! -- If reduced model, then need to copy from itemp(=>ibuff) to iarray
1337  if (this%nodes < this%nodesuser) then
1338  call this%fill_grid_array(itemp, iarray)
1339  end if
1340  !
Here is the call graph for this function:

◆ read_layer_array()

subroutine dismodule::read_layer_array ( class(distype this,
integer(i4b), dimension(maxbnd)  nodelist,
real(dp), dimension(ncolbnd, maxbnd), intent(inout)  darray,
integer(i4b), intent(in)  ncolbnd,
integer(i4b), intent(in)  maxbnd,
integer(i4b), intent(in)  icolbnd,
character(len=*), intent(in)  aname,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout 
)
private

For cells that are outside of the active domain, do not copy the array value into darray.

Definition at line 1407 of file Dis.f90.

1409  ! -- dummy
1410  class(DisType) :: this
1411  integer(I4B), intent(in) :: maxbnd
1412  integer(I4B), dimension(maxbnd) :: nodelist
1413  integer(I4B), intent(in) :: ncolbnd
1414  real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray
1415  integer(I4B), intent(in) :: icolbnd
1416  character(len=*), intent(in) :: aname
1417  integer(I4B), intent(in) :: inunit
1418  integer(I4B), intent(in) :: iout
1419  ! -- local
1420  integer(I4B) :: ir, ic, ncol, nrow, nlay, nval, ipos, nodeu
1421  !
1422  ! -- set variables
1423  nlay = this%mshape(1)
1424  nrow = this%mshape(2)
1425  ncol = this%mshape(3)
1426  !
1427  ! -- Read the array
1428  nval = ncol * nrow
1429  call readarray(inunit, this%dbuff, aname, this%ndim, ncol, nrow, nlay, &
1430  nval, iout, 0, 0)
1431  !
1432  ! -- Copy array into bound. Note that this routine was substantially
1433  ! changed on 9/21/2021 to support changes to READASARRAYS input
1434  ! for recharge and evapotranspiration. nodelist and bound are of
1435  ! size nrow * ncol and correspond directly to dbuff.
1436  ipos = 1
1437  do ir = 1, nrow
1438  do ic = 1, ncol
1439  nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
1440  darray(icolbnd, ipos) = this%dbuff(nodeu)
1441  ipos = ipos + 1
1442  !
1443  end do
1444  end do
1445  !
Here is the call graph for this function:

◆ record_array()

subroutine dismodule::record_array ( class(distype), intent(inout)  this,
real(dp), dimension(:), intent(inout), pointer, contiguous  darray,
integer(i4b), intent(in)  iout,
integer(i4b), intent(in)  iprint,
integer(i4b), intent(in)  idataun,
character(len=*), intent(in)  aname,
character(len=*), intent(in)  cdatafmp,
integer(i4b), intent(in)  nvaluesp,
integer(i4b), intent(in)  nwidthp,
character(len=*), intent(in)  editdesc,
real(dp), intent(in)  dinact 
)
private

The array is written to a formatted or unformatted external file depending on the arguments.

Parameters
[in,out]darraydouble precision array to record
[in]ioutascii output unit number
[in]iprintwhether to print the array
[in]idataunbinary output unit number, if negative don't write by layers, write entire array
[in]anametext descriptor
[in]cdatafmpwrite format
[in]nvaluespvalues per line
[in]nwidthpnumber width
[in]editdescformat type (I, G, F, S, E)
[in]dinactdouble precision value for cells excluded from model domain

Definition at line 1453 of file Dis.f90.

1455  ! -- dummy
1456  class(DisType), intent(inout) :: this
1457  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record
1458  integer(I4B), intent(in) :: iout !< ascii output unit number
1459  integer(I4B), intent(in) :: iprint !< whether to print the array
1460  integer(I4B), intent(in) :: idataun !< binary output unit number, if negative don't write by layers, write entire array
1461  character(len=*), intent(in) :: aname !< text descriptor
1462  character(len=*), intent(in) :: cdatafmp !< write format
1463  integer(I4B), intent(in) :: nvaluesp !< values per line
1464  integer(I4B), intent(in) :: nwidthp !< number width
1465  character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E)
1466  real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain
1467  ! -- local
1468  integer(I4B) :: k, ifirst
1469  integer(I4B) :: nlay
1470  integer(I4B) :: nrow
1471  integer(I4B) :: ncol
1472  integer(I4B) :: nval
1473  integer(I4B) :: nodeu, noder
1474  integer(I4B) :: istart, istop
1475  real(DP), dimension(:), pointer, contiguous :: dtemp
1476  ! -- formats
1477  character(len=*), parameter :: fmthsv = &
1478  "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, &
1479  &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)"
1480  !
1481  ! -- set variables
1482  nlay = this%mshape(1)
1483  nrow = this%mshape(2)
1484  ncol = this%mshape(3)
1485  !
1486  ! -- If this is a reduced model, then copy the values from darray into
1487  ! dtemp.
1488  if (this%nodes < this%nodesuser) then
1489  nval = this%nodes
1490  dtemp => this%dbuff
1491  do nodeu = 1, this%nodesuser
1492  noder = this%get_nodenumber(nodeu, 0)
1493  if (noder <= 0) then
1494  dtemp(nodeu) = dinact
1495  cycle
1496  end if
1497  dtemp(nodeu) = darray(noder)
1498  end do
1499  else
1500  nval = this%nodes
1501  dtemp => darray
1502  end if
1503  !
1504  ! -- Print to iout if iprint /= 0
1505  if (iprint /= 0) then
1506  istart = 1
1507  do k = 1, nlay
1508  istop = istart + nrow * ncol - 1
1509  call ulaprufw(ncol, nrow, kstp, kper, k, iout, dtemp(istart:istop), &
1510  aname, cdatafmp, nvaluesp, nwidthp, editdesc)
1511  istart = istop + 1
1512  end do
1513  end if
1514  !
1515  ! -- Save array to an external file.
1516  if (idataun > 0) then
1517  ! -- write to binary file by layer
1518  ifirst = 1
1519  istart = 1
1520  do k = 1, nlay
1521  istop = istart + nrow * ncol - 1
1522  if (ifirst == 1) write (iout, fmthsv) &
1523  trim(adjustl(aname)), idataun, &
1524  kstp, kper
1525  ifirst = 0
1526  call ulasav(dtemp(istart:istop), aname, kstp, kper, &
1527  pertim, totim, ncol, nrow, k, idataun)
1528  istart = istop + 1
1529  end do
1530  elseif (idataun < 0) then
1531  !
1532  ! -- write entire array as one record
1533  call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, &
1534  iout, delt, pertim, totim)
1535  end if
1536  !
Here is the call graph for this function:

◆ record_srcdst_list_header()

subroutine dismodule::record_srcdst_list_header ( class(distype this,
character(len=16), intent(in)  text,
character(len=16), intent(in)  textmodel,
character(len=16), intent(in)  textpackage,
character(len=16), intent(in)  dstmodel,
character(len=16), intent(in)  dstpackage,
integer(i4b), intent(in)  naux,
character(len=16), dimension(:), intent(in)  auxtxt,
integer(i4b), intent(in)  ibdchn,
integer(i4b), intent(in)  nlist,
integer(i4b), intent(in)  iout 
)
private

Definition at line 1541 of file Dis.f90.

1544  ! -- dummy
1545  class(DisType) :: this
1546  character(len=16), intent(in) :: text
1547  character(len=16), intent(in) :: textmodel
1548  character(len=16), intent(in) :: textpackage
1549  character(len=16), intent(in) :: dstmodel
1550  character(len=16), intent(in) :: dstpackage
1551  integer(I4B), intent(in) :: naux
1552  character(len=16), dimension(:), intent(in) :: auxtxt
1553  integer(I4B), intent(in) :: ibdchn
1554  integer(I4B), intent(in) :: nlist
1555  integer(I4B), intent(in) :: iout
1556  ! -- local
1557  integer(I4B) :: nlay, nrow, ncol
1558  !
1559  nlay = this%mshape(1)
1560  nrow = this%mshape(2)
1561  ncol = this%mshape(3)
1562  !
1563  ! -- Use ubdsv06 to write list header
1564  call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, &
1565  ibdchn, naux, auxtxt, ncol, nrow, nlay, &
1566  nlist, iout, delt, pertim, totim)
1567  !
Here is the call graph for this function:

◆ source_dimensions()

subroutine dismodule::source_dimensions ( class(distype this)
private

Definition at line 245 of file Dis.f90.

246  ! -- dummy
247  class(DisType) :: this
248  ! -- locals
249  integer(I4B) :: i, j, k
250  type(DisFoundType) :: found
251  !
252  ! -- update defaults with idm sourced values
253  call mem_set_value(this%nlay, 'NLAY', this%input_mempath, found%nlay)
254  call mem_set_value(this%nrow, 'NROW', this%input_mempath, found%nrow)
255  call mem_set_value(this%ncol, 'NCOL', this%input_mempath, found%ncol)
256  !
257  ! -- log simulation values
258  if (this%iout > 0) then
259  call this%log_dimensions(found)
260  end if
261  !
262  ! -- verify dimensions were set
263  if (this%nlay < 1) then
264  call store_error( &
265  'NLAY was not specified or was specified incorrectly.')
266  call store_error_filename(this%input_fname)
267  end if
268  if (this%nrow < 1) then
269  call store_error( &
270  'NROW was not specified or was specified incorrectly.')
271  call store_error_filename(this%input_fname)
272  end if
273  if (this%ncol < 1) then
274  call store_error( &
275  'NCOL was not specified or was specified incorrectly.')
276  call store_error_filename(this%input_fname)
277  end if
278  !
279  ! -- calculate nodesuser
280  this%nodesuser = this%nlay * this%nrow * this%ncol
281  !
282  ! -- Allocate delr, delc, and non-reduced vectors for dis
283  call mem_allocate(this%delr, this%ncol, 'DELR', this%memoryPath)
284  call mem_allocate(this%delc, this%nrow, 'DELC', this%memoryPath)
285  call mem_allocate(this%idomain, this%ncol, this%nrow, this%nlay, 'IDOMAIN', &
286  this%memoryPath)
287  call mem_allocate(this%top2d, this%ncol, this%nrow, 'TOP2D', this%memoryPath)
288  call mem_allocate(this%bot3d, this%ncol, this%nrow, this%nlay, 'BOT3D', &
289  this%memoryPath)
290  call mem_allocate(this%cellx, this%ncol, 'CELLX', this%memoryPath)
291  call mem_allocate(this%celly, this%nrow, 'CELLY', this%memoryPath)
292  !
293  ! -- initialize all cells to be active (idomain = 1)
294  do k = 1, this%nlay
295  do i = 1, this%nrow
296  do j = 1, this%ncol
297  this%idomain(j, i, k) = 1
298  end do
299  end do
300  end do
301  !
Here is the call graph for this function:

◆ source_griddata()

subroutine dismodule::source_griddata ( class(distype this)
private

Definition at line 331 of file Dis.f90.

332  ! -- dummy
333  class(DisType) :: this
334  type(DisFoundType) :: found
335  !
336  ! -- update defaults with idm sourced values
337  call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr)
338  call mem_set_value(this%delc, 'DELC', this%input_mempath, found%delc)
339  call mem_set_value(this%top2d, 'TOP', this%input_mempath, found%top)
340  call mem_set_value(this%bot3d, 'BOTM', this%input_mempath, found%botm)
341  call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
342  !
343  ! -- log simulation values
344  if (this%iout > 0) then
345  call this%log_griddata(found)
346  end if
347  !

◆ source_options()

subroutine dismodule::source_options ( class(distype this)
private

Definition at line 185 of file Dis.f90.

186  ! -- dummy
187  class(DisType) :: this
188  ! -- locals
189  character(len=LENVARNAME), dimension(3) :: lenunits = &
190  &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS']
191  type(disfoundtype) :: found
192  !
193  ! -- update defaults with idm sourced values
194  call mem_set_value(this%lenuni, 'LENGTH_UNITS', this%input_mempath, &
195  lenunits, found%length_units)
196  call mem_set_value(this%nogrb, 'NOGRB', this%input_mempath, found%nogrb)
197  call mem_set_value(this%xorigin, 'XORIGIN', this%input_mempath, found%xorigin)
198  call mem_set_value(this%yorigin, 'YORIGIN', this%input_mempath, found%yorigin)
199  call mem_set_value(this%angrot, 'ANGROT', this%input_mempath, found%angrot)
200  !
201  ! -- log values to list file
202  if (this%iout > 0) then
203  call this%log_options(found)
204  end if
205  !

◆ supports_layers()

logical function dismodule::supports_layers ( class(distype this)
private

Definition at line 1048 of file Dis.f90.

1049  ! -- dummy
1050  class(DisType) :: this
1051  !
1052  supports_layers = .true.
1053  !

◆ write_grb()

subroutine dismodule::write_grb ( class(distype this,
integer(i4b), dimension(:), intent(in)  icelltype 
)

Definition at line 549 of file Dis.f90.

550  ! -- modules
551  use openspecmodule, only: access, form
552  use constantsmodule, only: lenbigline
553  ! -- dummy
554  class(DisType) :: this
555  integer(I4B), dimension(:), intent(in) :: icelltype
556  ! -- local
557  integer(I4B) :: iunit, ntxt, ncpl, version
558  integer(I4B), parameter :: lentxt = 100
559  character(len=50) :: txthdr
560  character(len=lentxt) :: txt
561  character(len=LINELENGTH) :: fname
562  character(len=LENBIGLINE) :: crs
563  logical(LGP) :: found_crs
564  character(len=*), parameter :: fmtgrdsave = &
565  "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', &
566  &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)"
567  !
568  ! -- Initialize
569  version = 1
570  ntxt = 16
571  ncpl = this%nrow * this%ncol
572  !
573  call mem_set_value(crs, 'CRS', this%input_mempath, found_crs)
574  !
575  ! -- set version
576  if (found_crs) then
577  ntxt = ntxt + 1
578  version = 2
579  end if
580  !
581  ! -- Open the file
582  fname = trim(this%output_fname)
583  iunit = getunit()
584  write (this%iout, fmtgrdsave) iunit, trim(adjustl(fname))
585  call openfile(iunit, this%iout, trim(adjustl(fname)), 'DATA(BINARY)', &
586  form, access, 'REPLACE')
587  !
588  ! -- write header information
589  write (txthdr, '(a)') 'GRID DIS'
590  txthdr(50:50) = new_line('a')
591  write (iunit) txthdr
592  write (txthdr, '(a, i0)') 'VERSION ', version
593  txthdr(50:50) = new_line('a')
594  write (iunit) txthdr
595  write (txthdr, '(a, i0)') 'NTXT ', ntxt
596  txthdr(50:50) = new_line('a')
597  write (iunit) txthdr
598  write (txthdr, '(a, i0)') 'LENTXT ', lentxt
599  txthdr(50:50) = new_line('a')
600  write (iunit) txthdr
601  !
602  ! -- write variable definitions
603  write (txt, '(3a, i0)') 'NCELLS ', 'INTEGER ', 'NDIM 0 # ', this%nodesuser
604  txt(lentxt:lentxt) = new_line('a')
605  write (iunit) txt
606  write (txt, '(3a, i0)') 'NLAY ', 'INTEGER ', 'NDIM 0 # ', this%nlay
607  txt(lentxt:lentxt) = new_line('a')
608  write (iunit) txt
609  write (txt, '(3a, i0)') 'NROW ', 'INTEGER ', 'NDIM 0 # ', this%nrow
610  txt(lentxt:lentxt) = new_line('a')
611  write (iunit) txt
612  write (txt, '(3a, i0)') 'NCOL ', 'INTEGER ', 'NDIM 0 # ', this%ncol
613  txt(lentxt:lentxt) = new_line('a')
614  write (iunit) txt
615  write (txt, '(3a, i0)') 'NJA ', 'INTEGER ', 'NDIM 0 # ', this%nja
616  txt(lentxt:lentxt) = new_line('a')
617  write (iunit) txt
618  write (txt, '(3a, 1pg24.15)') 'XORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%xorigin
619  txt(lentxt:lentxt) = new_line('a')
620  write (iunit) txt
621  write (txt, '(3a, 1pg24.15)') 'YORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%yorigin
622  txt(lentxt:lentxt) = new_line('a')
623  write (iunit) txt
624  write (txt, '(3a, 1pg24.15)') 'ANGROT ', 'DOUBLE ', 'NDIM 0 # ', this%angrot
625  txt(lentxt:lentxt) = new_line('a')
626  write (iunit) txt
627  write (txt, '(3a, i0)') 'DELR ', 'DOUBLE ', 'NDIM 1 ', this%ncol
628  txt(lentxt:lentxt) = new_line('a')
629  write (iunit) txt
630  write (txt, '(3a, i0)') 'DELC ', 'DOUBLE ', 'NDIM 1 ', this%nrow
631  txt(lentxt:lentxt) = new_line('a')
632  write (iunit) txt
633  write (txt, '(3a, i0)') 'TOP ', 'DOUBLE ', 'NDIM 1 ', ncpl
634  txt(lentxt:lentxt) = new_line('a')
635  write (iunit) txt
636  write (txt, '(3a, i0)') 'BOTM ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser
637  txt(lentxt:lentxt) = new_line('a')
638  write (iunit) txt
639  write (txt, '(3a, i0)') 'IA ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1
640  txt(lentxt:lentxt) = new_line('a')
641  write (iunit) txt
642  write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', size(this%con%jausr)
643  txt(lentxt:lentxt) = new_line('a')
644  write (iunit) txt
645  write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
646  txt(lentxt:lentxt) = new_line('a')
647  write (iunit) txt
648  write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
649  txt(lentxt:lentxt) = new_line('a')
650  write (iunit) txt
651  !
652  ! -- if version 2 write character array headers
653  if (version == 2) then
654  if (found_crs) then
655  write (txt, '(3a, i0)') 'CRS ', 'CHARACTER ', 'NDIM 1 ', &
656  len_trim(crs)
657  txt(lentxt:lentxt) = new_line('a')
658  write (iunit) txt
659  end if
660  end if
661  !
662  ! -- write data
663  write (iunit) this%nodesuser ! ncells
664  write (iunit) this%nlay ! nlay
665  write (iunit) this%nrow ! nrow
666  write (iunit) this%ncol ! ncol
667  write (iunit) this%nja ! nja
668  write (iunit) this%xorigin ! xorigin
669  write (iunit) this%yorigin ! yorigin
670  write (iunit) this%angrot ! angrot
671  write (iunit) this%delr ! delr
672  write (iunit) this%delc ! delc
673  write (iunit) this%top2d ! top2d
674  write (iunit) this%bot3d ! bot3d
675  write (iunit) this%con%iausr ! iausr
676  write (iunit) this%con%jausr ! jausr
677  write (iunit) this%idomain ! idomain
678  write (iunit) icelltype ! icelltype
679  !
680  ! -- if version 2 write character array data
681  if (version == 2) then
682  if (found_crs) write (iunit) trim(crs) ! crs user input
683  end if
684  !
685  ! -- Close the file
686  close (iunit)
687  !
integer(i4b), parameter lenbigline
maximum length of a big line
Definition: Constants.f90:15
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
Here is the call graph for this function: