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

Data Types

type  disvtype
 Vertex grid discretization. More...
 
type  disvfoundtype
 

Functions/Subroutines

subroutine, public disv_cr (dis, name_model, input_mempath, inunit, iout)
 Create a new discretization by vertices object. More...
 
subroutine disv_load (this)
 Transfer IDM data into this discretization object. More...
 
subroutine disv_df (this)
 Define the discretization. More...
 
subroutine disv_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 griddata found to list file. More...
 
subroutine grid_finalize (this)
 Finalize grid (check properties, allocate arrays, compute connections) More...
 
subroutine source_vertices (this)
 Load grid vertices from IDM into package. More...
 
subroutine define_cellverts (this, icell2d, ncvert, icvert)
 Build data structures to hold cell vertex info. More...
 
subroutine source_cell2d (this)
 Copy cell2d data from IDM into package. More...
 
subroutine connect (this)
 Build grid 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,j) More...
 
subroutine nodeu_to_array (this, nodeu, arr)
 Convert a user nodenumber to an array (nodenumber) or (k,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_idx2 (this, k, j, icheck)
 Get reduced node number from layer and within-layer node indices. 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 allocate_scalars (this, name_model, input_mempath)
 Allocate and initialize scalars. More...
 
subroutine allocate_arrays (this)
 Allocate and initialize arrays. More...
 
real(dp) function get_cell2d_area (this, icell2d)
 Get the signed area of the cell. 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)
 Get number of cells per layer (ncpl) More...
 
subroutine get_polyverts (this, ic, polyverts, closed)
 Get a 2D array of polygon vertices, listed in clockwise order beginning with the lower left corner. 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 disvmodule::allocate_arrays ( class(disvtype this)
private

Definition at line 1156 of file Disv.f90.

1157  ! -- dummy
1158  class(DisvType) :: this
1159  !
1160  ! -- Allocate arrays in DisBaseType (mshape, top, bot, area)
1161  call this%DisBaseType%allocate_arrays()
1162  !
1163  ! -- Allocate arrays for DisvType
1164  if (this%nodes < this%nodesuser) then
1165  call mem_allocate(this%nodeuser, this%nodes, 'NODEUSER', this%memoryPath)
1166  call mem_allocate(this%nodereduced, this%nodesuser, 'NODEREDUCED', &
1167  this%memoryPath)
1168  else
1169  call mem_allocate(this%nodeuser, 1, 'NODEUSER', this%memoryPath)
1170  call mem_allocate(this%nodereduced, 1, 'NODEREDUCED', this%memoryPath)
1171  end if
1172  ! -- Initialize
1173  this%mshape(1) = this%nlay
1174  this%mshape(2) = this%ncpl
1175  !

◆ allocate_scalars()

subroutine disvmodule::allocate_scalars ( class(disvtype this,
character(len=*), intent(in)  name_model,
character(len=*), intent(in)  input_mempath 
)

Definition at line 1132 of file Disv.f90.

1133  ! -- dummy
1134  class(DisvType) :: this
1135  character(len=*), intent(in) :: name_model
1136  character(len=*), intent(in) :: input_mempath
1137  !
1138  ! -- Allocate parent scalars
1139  call this%DisBaseType%allocate_scalars(name_model, input_mempath)
1140  !
1141  ! -- Allocate
1142  call mem_allocate(this%nlay, 'NLAY', this%memoryPath)
1143  call mem_allocate(this%ncpl, 'NCPL', this%memoryPath)
1144  call mem_allocate(this%nvert, 'NVERT', this%memoryPath)
1145  !
1146  ! -- Initialize
1147  this%nlay = 0
1148  this%ncpl = 0
1149  this%nvert = 0
1150  this%ndim = 2
1151  !

◆ connect()

subroutine disvmodule::connect ( class(disvtype this)
private

Definition at line 639 of file Disv.f90.

640  ! -- dummy
641  class(DisvType) :: this
642  ! -- local
643  integer(I4B) :: j, k
644  integer(I4B) :: noder, nrsize
645  integer(I4B) :: narea_eq_zero
646  integer(I4B) :: narea_lt_zero
647  real(DP) :: area
648  !
649  ! -- Initialize
650  narea_eq_zero = 0
651  narea_lt_zero = 0
652  !
653  ! -- Assign the cell area
654  do j = 1, this%ncpl
655  area = this%get_cell2d_area(j)
656  do k = 1, this%nlay
657  noder = this%get_nodenumber(k, j, 0)
658  if (noder > 0) this%area(noder) = area
659  end do
660  if (area < dzero) then
661  narea_lt_zero = narea_lt_zero + 1
662  write (errmsg, '(a,i0,a)') &
663  &'Calculated CELL2D area less than zero for cell ', j, '.'
664  call store_error(errmsg)
665  end if
666  if (area == dzero) then
667  narea_eq_zero = narea_eq_zero + 1
668  write (errmsg, '(a,i0,a)') &
669  'Calculated CELL2D area is zero for cell ', j, '.'
670  call store_error(errmsg)
671  end if
672  end do
673  !
674  ! -- check for errors
675  if (count_errors() > 0) then
676  if (narea_lt_zero > 0) then
677  write (errmsg, '(i0,a)') narea_lt_zero, &
678  ' cell(s) have an area less than zero. Calculated cell &
679  &areas must be greater than zero. Negative areas often &
680  &mean vertices are not listed in clockwise order.'
681  call store_error(errmsg)
682  end if
683  if (narea_eq_zero > 0) then
684  write (errmsg, '(i0,a)') narea_eq_zero, &
685  ' cell(s) have an area equal to zero. Calculated cell &
686  &areas must be greater than zero. Calculated cell &
687  &areas equal to zero indicate that the cell is not defined &
688  &by a valid polygon.'
689  call store_error(errmsg)
690  end if
691  call store_error_filename(this%input_fname)
692  end if
693  !
694  ! -- create and fill the connections object
695  nrsize = 0
696  if (this%nodes < this%nodesuser) nrsize = this%nodes
697  allocate (this%con)
698  call this%con%disvconnections(this%name_model, this%nodes, &
699  this%ncpl, this%nlay, nrsize, &
700  this%nvert, this%vertices, this%iavert, &
701  this%javert, this%cellxy, &
702  this%top, this%bot, &
703  this%nodereduced, this%nodeuser)
704  this%nja = this%con%nja
705  this%njas = this%con%njas
706  !
Here is the call graph for this function:

◆ connection_normal()

subroutine disvmodule::connection_normal ( class(disvtype 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 1008 of file Disv.f90.

1010  ! -- dummy
1011  class(DisvType) :: this
1012  integer(I4B), intent(in) :: noden !< cell (reduced nn)
1013  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
1014  integer(I4B), intent(in) :: ihc !< horizontal connection flag
1015  real(DP), intent(inout) :: xcomp
1016  real(DP), intent(inout) :: ycomp
1017  real(DP), intent(inout) :: zcomp
1018  integer(I4B), intent(in) :: ipos
1019  ! -- local
1020  real(DP) :: angle, dmult
1021  !
1022  ! -- Set vector components based on ihc
1023  if (ihc == 0) then
1024  xcomp = dzero
1025  ycomp = dzero
1026  if (nodem < noden) then
1027  !
1028  ! -- nodem must be above noden, so upward connection
1029  zcomp = done
1030  else
1031  !
1032  ! -- nodem must be below noden, so downward connection
1033  zcomp = -done
1034  end if
1035  else
1036  ! -- find from anglex, since anglex is symmetric, need to flip vector
1037  ! for lower triangle (nodem < noden)
1038  !ipos = this%con%getjaindex(noden, nodem)
1039  angle = this%con%anglex(this%con%jas(ipos))
1040  dmult = done
1041  if (nodem < noden) dmult = -done
1042  xcomp = cos(angle) * dmult
1043  ycomp = sin(angle) * dmult
1044  zcomp = dzero
1045  end if
1046  !

◆ connection_vector()

subroutine disvmodule::connection_vector ( class(disvtype 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 1054 of file Disv.f90.

1056  ! -- dummy
1057  class(DisvType) :: this
1058  integer(I4B), intent(in) :: noden !< cell (reduced nn)
1059  integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
1060  logical, intent(in) :: nozee
1061  real(DP), intent(in) :: satn
1062  real(DP), intent(in) :: satm
1063  integer(I4B), intent(in) :: ihc !< horizontal connection flag
1064  real(DP), intent(inout) :: xcomp
1065  real(DP), intent(inout) :: ycomp
1066  real(DP), intent(inout) :: zcomp
1067  real(DP), intent(inout) :: conlen
1068  ! -- local
1069  integer(I4B) :: nodeu, ncell2d, mcell2d, k
1070  real(DP) :: xn, xm, yn, ym, zn, zm
1071  !
1072  ! -- Set vector components based on ihc
1073  if (ihc == 0) then
1074  !
1075  ! -- vertical connection; set zcomp positive upward
1076  xcomp = dzero
1077  ycomp = dzero
1078  if (nodem < noden) then
1079  zcomp = done
1080  else
1081  zcomp = -done
1082  end if
1083  zn = this%bot(noden) + dhalf * (this%top(noden) - this%bot(noden))
1084  zm = this%bot(nodem) + dhalf * (this%top(nodem) - this%bot(nodem))
1085  conlen = abs(zm - zn)
1086  else
1087  !
1088  ! -- horizontal connection, with possible z component due to cell offsets
1089  ! and/or water table conditions
1090  if (nozee) then
1091  zn = dzero
1092  zm = dzero
1093  else
1094  zn = this%bot(noden) + dhalf * satn * (this%top(noden) - this%bot(noden))
1095  zm = this%bot(nodem) + dhalf * satm * (this%top(nodem) - this%bot(nodem))
1096  end if
1097  nodeu = this%get_nodeuser(noden)
1098  call get_jk(nodeu, this%ncpl, this%nlay, ncell2d, k)
1099  nodeu = this%get_nodeuser(nodem)
1100  call get_jk(nodeu, this%ncpl, this%nlay, mcell2d, k)
1101  xn = this%cellxy(1, ncell2d)
1102  yn = this%cellxy(2, ncell2d)
1103  xm = this%cellxy(1, mcell2d)
1104  ym = this%cellxy(2, mcell2d)
1105  call line_unit_vector(xn, yn, zn, xm, ym, zm, xcomp, ycomp, zcomp, &
1106  conlen)
1107  end if
1108  !
Here is the call graph for this function:

◆ define_cellverts()

subroutine disvmodule::define_cellverts ( class(disvtype this,
integer(i4b), dimension(:), intent(in), pointer, contiguous  icell2d,
integer(i4b), dimension(:), intent(in), pointer, contiguous  ncvert,
integer(i4b), dimension(:), intent(in), pointer, contiguous  icvert 
)
private

Definition at line 550 of file Disv.f90.

551  ! -- modules
552  use sparsemodule, only: sparsematrix
553  ! -- dummy
554  class(DisvType) :: this
555  integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icell2d
556  integer(I4B), dimension(:), contiguous, pointer, intent(in) :: ncvert
557  integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icvert
558  ! -- locals
559  type(sparsematrix) :: vert_spm
560  integer(I4B) :: i, j, ierr
561  integer(I4B) :: icv_idx, startvert, maxnnz = 5
562  !
563  ! -- initialize sparse matrix
564  call vert_spm%init(this%ncpl, this%nvert, maxnnz)
565  !
566  ! -- add sparse matrix connections from input memory paths
567  icv_idx = 1
568  do i = 1, this%ncpl
569  if (icell2d(i) /= i) call store_error('ICELL2D input sequence violation.')
570  do j = 1, ncvert(i)
571  call vert_spm%addconnection(i, icvert(icv_idx), 0)
572  if (j == 1) then
573  startvert = icvert(icv_idx)
574  elseif (j == ncvert(i) .and. (icvert(icv_idx) /= startvert)) then
575  call vert_spm%addconnection(i, startvert, 0)
576  end if
577  icv_idx = icv_idx + 1
578  end do
579  end do
580  !
581  ! -- allocate and fill iavert and javert
582  call mem_allocate(this%iavert, this%ncpl + 1, 'IAVERT', this%memoryPath)
583  call mem_allocate(this%javert, vert_spm%nnz, 'JAVERT', this%memoryPath)
584  call vert_spm%filliaja(this%iavert, this%javert, ierr)
585  call vert_spm%destroy()
586  !
Here is the call graph for this function:

◆ disv_cr()

subroutine, public disvmodule::disv_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 110 of file Disv.f90.

111  ! -- dummy
112  class(DisBaseType), pointer :: dis
113  character(len=*), intent(in) :: name_model
114  character(len=*), intent(in) :: input_mempath
115  integer(I4B), intent(in) :: inunit
116  integer(I4B), intent(in) :: iout
117  ! -- local
118  type(DisvType), pointer :: disnew
119  ! -- formats
120  character(len=*), parameter :: fmtheader = &
121  "(1X, /1X, 'DISV -- VERTEX GRID DISCRETIZATION PACKAGE,', &
122  &' VERSION 1 : 12/23/2015 - INPUT READ FROM MEMPATH: ', A, //)"
123  !
124  allocate (disnew)
125  dis => disnew
126  call disnew%allocate_scalars(name_model, input_mempath)
127  dis%inunit = inunit
128  dis%iout = iout
129  !
130  ! -- If disv enabled
131  if (inunit > 0) then
132  !
133  ! -- Identify package
134  if (iout > 0) then
135  write (iout, fmtheader) dis%input_mempath
136  end if
137  !
138  ! -- load disv
139  call disnew%disv_load()
140  end if
141  !
Here is the caller graph for this function:

◆ disv_da()

subroutine disvmodule::disv_da ( class(disvtype this)
private

Definition at line 171 of file Disv.f90.

172  ! -- dummy
173  class(DisvType) :: this
174  !
175  ! -- Deallocate idm memory
176  call memorystore_remove(this%name_model, 'DISV', idm_context)
177  call memorystore_remove(component=this%name_model, &
178  context=idm_context)
179  !
180  ! -- DisBaseType deallocate
181  call this%DisBaseType%dis_da()
182  !
183  ! -- Deallocate scalars
184  call mem_deallocate(this%nlay)
185  call mem_deallocate(this%ncpl)
186  call mem_deallocate(this%nvert)
187  !
188  ! -- Deallocate Arrays
189  call mem_deallocate(this%nodereduced)
190  call mem_deallocate(this%nodeuser)
191  call mem_deallocate(this%vertices)
192  call mem_deallocate(this%cellxy)
193  call mem_deallocate(this%iavert)
194  call mem_deallocate(this%javert)
195  call mem_deallocate(this%top1d)
196  call mem_deallocate(this%bot2d)
197  call mem_deallocate(this%idomain)
198  !
Here is the call graph for this function:

◆ disv_df()

subroutine disvmodule::disv_df ( class(disvtype this)
private

Definition at line 161 of file Disv.f90.

162  ! -- dummy
163  class(DisvType) :: this
164  !
165  call this%grid_finalize()
166  !

◆ disv_load()

subroutine disvmodule::disv_load ( class(disvtype this)
private

Definition at line 146 of file Disv.f90.

147  ! -- dummy
148  class(DisvType) :: this
149  !
150  ! -- source input data
151  call this%source_options()
152  call this%source_dimensions()
153  call this%source_griddata()
154  call this%source_vertices()
155  call this%source_cell2d()
156  !

◆ get_cell2d_area()

real(dp) function disvmodule::get_cell2d_area ( class(disvtype this,
integer(i4b), intent(in)  icell2d 
)
private

A negative result means points are in counter-clockwise orientation. Area is computed from the formula: a = 1/2 *[(x1*y2 + x2*y3 + x3*y4 + ... + xn*y1) - (x2*y1 + x3*y2 + x4*y3 + ... + x1*yn)]

Definition at line 1185 of file Disv.f90.

1186  ! -- dummy
1187  class(DisvType) :: this
1188  integer(I4B), intent(in) :: icell2d
1189  ! -- return
1190  real(DP) :: area
1191  ! -- local
1192  integer(I4B) :: ivert
1193  integer(I4B) :: nvert
1194  integer(I4B) :: icount
1195  integer(I4B) :: iv1
1196  real(DP) :: x
1197  real(DP) :: y
1198  real(DP) :: x1
1199  real(DP) :: y1
1200  !
1201  area = dzero
1202  nvert = this%iavert(icell2d + 1) - this%iavert(icell2d)
1203  icount = 1
1204  iv1 = this%javert(this%iavert(icell2d))
1205  x1 = this%vertices(1, iv1)
1206  y1 = this%vertices(2, iv1)
1207  do ivert = this%iavert(icell2d), this%iavert(icell2d + 1) - 1
1208  x = this%vertices(1, this%javert(ivert))
1209  if (icount < nvert) then
1210  y = this%vertices(2, this%javert(ivert + 1))
1211  else
1212  y = this%vertices(2, this%javert(this%iavert(icell2d)))
1213  end if
1214  area = area + (x - x1) * (y - y1)
1215  icount = icount + 1
1216  end do
1217  !
1218  icount = 1
1219  do ivert = this%iavert(icell2d), this%iavert(icell2d + 1) - 1
1220  y = this%vertices(2, this%javert(ivert))
1221  if (icount < nvert) then
1222  x = this%vertices(1, this%javert(ivert + 1))
1223  else
1224  x = this%vertices(1, this%javert(this%iavert(icell2d)))
1225  end if
1226  area = area - (x - x1) * (y - y1)
1227  icount = icount + 1
1228  end do
1229  !
1230  area = -done * area * dhalf
1231  !

◆ get_dis_enum()

integer(i4b) function disvmodule::get_dis_enum ( class(disvtype), intent(in)  this)
private

Definition at line 1123 of file Disv.f90.

1124  use constantsmodule, only: disv
1125  class(DisvType), intent(in) :: this
1126  integer(I4B) :: dis_enum
1127  dis_enum = disv
This module contains simulation constants.
Definition: Constants.f90:9
@ disv
DISU6 discretization.
Definition: Constants.f90:156

◆ get_dis_type()

subroutine disvmodule::get_dis_type ( class(disvtype), intent(in)  this,
character(len=*), intent(out)  dis_type 
)
private

Definition at line 1113 of file Disv.f90.

1114  ! -- dummy
1115  class(DisvType), intent(in) :: this
1116  character(len=*), intent(out) :: dis_type
1117  !
1118  dis_type = "DISV"
1119  !

◆ get_max_npolyverts()

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

Definition at line 1489 of file Disv.f90.

1490  class(DisvType), intent(inout) :: this
1491  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
1492  integer(I4B) :: max_npolyverts
1493  integer(I4B) :: ic
1494 
1495  max_npolyverts = 0
1496  do ic = 1, this%nodes
1497  max_npolyverts = max(max_npolyverts, this%get_npolyverts(ic, closed))
1498  end do

◆ get_ncpl()

integer(i4b) function disvmodule::get_ncpl ( class(disvtype this)
private

Definition at line 1414 of file Disv.f90.

1415  ! -- return
1416  integer(I4B) :: get_ncpl
1417  ! -- dummy
1418  class(DisvType) :: this
1419  !
1420  get_ncpl = this%ncpl
1421  !

◆ get_nodenumber_idx1()

integer(i4b) function disvmodule::get_nodenumber_idx1 ( class(disvtype), intent(in)  this,
integer(i4b), intent(in)  nodeu,
integer(i4b), intent(in)  icheck 
)
private

Definition at line 920 of file Disv.f90.

921  ! -- return
922  integer(I4B) :: nodenumber
923  ! -- dummy
924  class(DisvType), intent(in) :: this
925  integer(I4B), intent(in) :: nodeu
926  integer(I4B), intent(in) :: icheck
927  ! -- local
928  !
929  ! -- check the node number if requested
930  if (icheck /= 0) then
931  !
932  ! -- If within valid range, convert to reduced nodenumber
933  if (nodeu < 1 .or. nodeu > this%nodesuser) then
934  nodenumber = 0
935  write (errmsg, '(a,i0,a,i0,a)') &
936  'Node number (', nodeu, ') is less than 1 or greater than nodes (', &
937  this%nodesuser, ').'
938  call store_error(errmsg)
939  else
940  nodenumber = nodeu
941  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
942  end if
943  else
944  nodenumber = nodeu
945  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
946  end if
947  !
Here is the call graph for this function:

◆ get_nodenumber_idx2()

integer(i4b) function disvmodule::get_nodenumber_idx2 ( class(disvtype), intent(in)  this,
integer(i4b), intent(in)  k,
integer(i4b), intent(in)  j,
integer(i4b), intent(in)  icheck 
)
private

Definition at line 952 of file Disv.f90.

953  ! -- return
954  integer(I4B) :: nodenumber
955  ! -- dummy
956  class(DisvType), intent(in) :: this
957  integer(I4B), intent(in) :: k, j
958  integer(I4B), intent(in) :: icheck
959  ! -- local
960  integer(I4B) :: nodeu
961  ! -- formats
962  character(len=*), parameter :: fmterr = &
963  &"('Error in disv grid cell indices: layer = ',i0,', node = ',i0)"
964  !
965  nodeu = get_node(k, 1, j, this%nlay, 1, this%ncpl)
966  if (nodeu < 1) then
967  write (errmsg, fmterr) k, j
968  call store_error(errmsg, terminate=.true.)
969  end if
970  nodenumber = nodeu
971  if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu)
972  !
973  ! -- check the node number if requested
974  if (icheck /= 0) then
975  !
976  errmsg = ""
977  !
978  if (k < 1 .or. k > this%nlay) then
979  write (errmsg, '(a,i0,a)') &
980  'Layer number in list (', k, ') is outside of the grid.'
981  end if
982  if (j < 1 .or. j > this%ncpl) then
983  write (errmsg, '(a,1x,a,i0,a)') &
984  trim(adjustl(errmsg)), 'Node number in list (', j, &
985  ') is outside of the grid.'
986  end if
987  !
988  ! -- Error if outside of range
989  if (nodeu < 1 .or. nodeu > this%nodesuser) then
990  write (errmsg, '(a,1x,a,i0,a,i0,a)') &
991  trim(adjustl(errmsg)), &
992  'Node number (', nodeu, ') is less than 1 or greater '// &
993  'than nodes (', this%nodesuser, ').'
994  end if
995  !
996  if (len_trim(adjustl(errmsg)) > 0) then
997  call store_error(errmsg)
998  end if
999  !
1000  end if
1001  !
Here is the call graph for this function:

◆ get_npolyverts()

integer(i4b) function disvmodule::get_npolyverts ( class(disvtype), 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 1471 of file Disv.f90.

1472  class(DisvType), intent(inout) :: this
1473  integer(I4B), intent(in) :: ic !< cell number (reduced)
1474  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
1475  integer(I4B) :: npolyverts
1476  ! local
1477  integer(I4B) :: ncpl, icu, icu2d
1478 
1479  ncpl = this%get_ncpl()
1480  icu = this%get_nodeuser(ic)
1481  icu2d = icu - ((icu - 1) / ncpl) * ncpl
1482  npolyverts = this%iavert(icu2d + 1) - this%iavert(icu2d) - 1
1483  if (present(closed)) then
1484  if (closed) npolyverts = npolyverts + 1
1485  end if

◆ get_polyverts()

subroutine disvmodule::get_polyverts ( class(disvtype), intent(inout)  this,
integer(i4b), intent(in)  ic,
real(dp), dimension(:, :), intent(out), allocatable  polyverts,
logical(lgp), intent(in), optional  closed 
)
private
Parameters
[in]iccell number (reduced)
[out]polyvertspolygon vertices (column-major indexing)
[in]closedwhether to close the polygon, duplicating a vertex (default false)

Definition at line 1427 of file Disv.f90.

1428  ! -- dummy
1429  class(DisvType), intent(inout) :: this
1430  integer(I4B), intent(in) :: ic !< cell number (reduced)
1431  real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices (column-major indexing)
1432  logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex (default false)
1433  ! -- local
1434  integer(I4B) :: icu, icu2d, iavert, ncpl, nverts, m, j
1435  logical(LGP) :: lclosed
1436  !
1437  ! count vertices
1438  ncpl = this%get_ncpl()
1439  icu = this%get_nodeuser(ic)
1440  icu2d = icu - ((icu - 1) / ncpl) * ncpl
1441  nverts = this%iavert(icu2d + 1) - this%iavert(icu2d) - 1
1442  !
1443  ! check closed option
1444  if (.not. (present(closed))) then
1445  lclosed = .false.
1446  else
1447  lclosed = closed
1448  end if
1449  !
1450  ! allocate vertices array
1451  if (lclosed) then
1452  allocate (polyverts(2, nverts + 1))
1453  else
1454  allocate (polyverts(2, nverts))
1455  end if
1456  !
1457  ! set vertices
1458  iavert = this%iavert(icu2d)
1459  do m = 1, nverts
1460  j = this%javert(iavert - 1 + m)
1461  polyverts(:, m) = (/this%vertices(1, j), this%vertices(2, j)/)
1462  end do
1463  !
1464  ! close if enabled
1465  if (lclosed) &
1466  polyverts(:, nverts + 1) = polyverts(:, 1)
1467  !

◆ grid_finalize()

subroutine disvmodule::grid_finalize ( class(disvtype this)
private

Definition at line 392 of file Disv.f90.

393  ! -- dummy
394  class(DisvType) :: this
395  ! -- locals
396  integer(I4B) :: node, noder, j, k
397  real(DP) :: top
398  real(DP) :: dz
399  ! -- formats
400  character(len=*), parameter :: fmtdz = &
401  "('CELL (',i0,',',i0,') THICKNESS <= 0. ', &
402  &'TOP, BOT: ',2(1pg24.15))"
403  character(len=*), parameter :: fmtnr = &
404  "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',&
405  &/1x, 'Number of user nodes: ',I0,&
406  &/1X, 'Number of nodes in solution: ', I0, //)"
407  !
408  ! -- count active cells
409  this%nodes = 0
410  do k = 1, this%nlay
411  do j = 1, this%ncpl
412  if (this%idomain(j, k) > 0) this%nodes = this%nodes + 1
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  do k = 1, this%nlay
426  do j = 1, this%ncpl
427  if (this%idomain(j, k) == 0) cycle
428  if (this%idomain(j, k) > 0) then
429  if (k > 1) then
430  top = this%bot2d(j, k - 1)
431  else
432  top = this%top1d(j)
433  end if
434  dz = top - this%bot2d(j, k)
435  if (dz <= dzero) then
436  write (errmsg, fmt=fmtdz) k, j, top, this%bot2d(j, k)
437  call store_error(errmsg)
438  end if
439  end if
440  end do
441  end do
442  if (count_errors() > 0) then
443  call store_error_filename(this%input_fname)
444  end if
445  !
446  ! -- Write message if reduced grid
447  if (this%nodes < this%nodesuser) then
448  write (this%iout, fmtnr) this%nodesuser, this%nodes
449  end if
450  !
451  ! -- Array size is now known, so allocate
452  call this%allocate_arrays()
453  !
454  ! -- Fill the nodereduced array with the reduced nodenumber, or
455  ! a negative number to indicate it is a pass-through cell, or
456  ! a zero to indicate that the cell is excluded from the
457  ! solution.
458  if (this%nodes < this%nodesuser) then
459  node = 1
460  noder = 1
461  do k = 1, this%nlay
462  do j = 1, this%ncpl
463  if (this%idomain(j, k) > 0) then
464  this%nodereduced(node) = noder
465  noder = noder + 1
466  elseif (this%idomain(j, k) < 0) then
467  this%nodereduced(node) = -1
468  else
469  this%nodereduced(node) = 0
470  end if
471  node = node + 1
472  end do
473  end do
474  end if
475  !
476  ! -- allocate and fill nodeuser if a reduced grid
477  if (this%nodes < this%nodesuser) then
478  node = 1
479  noder = 1
480  do k = 1, this%nlay
481  do j = 1, this%ncpl
482  if (this%idomain(j, k) > 0) then
483  this%nodeuser(noder) = node
484  noder = noder + 1
485  end if
486  node = node + 1
487  end do
488  end do
489  end if
490  !
491  ! -- Move top1d and bot2d into top and bot
492  ! and set x and y center coordinates
493  node = 0
494  do k = 1, this%nlay
495  do j = 1, this%ncpl
496  node = node + 1
497  noder = node
498  if (this%nodes < this%nodesuser) noder = this%nodereduced(node)
499  if (noder <= 0) cycle
500  if (k > 1) then
501  top = this%bot2d(j, k - 1)
502  else
503  top = this%top1d(j)
504  end if
505  this%top(noder) = top
506  this%bot(noder) = this%bot2d(j, k)
507  this%xc(noder) = this%cellxy(1, j)
508  this%yc(noder) = this%cellxy(2, j)
509  end do
510  end do
511  !
512  ! -- Build connections
513  call this%connect()
514  !
Here is the call graph for this function:

◆ log_dimensions()

subroutine disvmodule::log_dimensions ( class(disvtype this,
type(disvfoundtype), intent(in)  found 
)
private

Definition at line 322 of file Disv.f90.

323  ! -- dummy
324  class(DisvType) :: this
325  type(DisvFoundType), intent(in) :: found
326  !
327  write (this%iout, '(1x,a)') 'Setting Discretization Dimensions'
328  !
329  if (found%nlay) then
330  write (this%iout, '(4x,a,i0)') 'NLAY = ', this%nlay
331  end if
332  !
333  if (found%ncpl) then
334  write (this%iout, '(4x,a,i0)') 'NCPL = ', this%ncpl
335  end if
336  !
337  if (found%nvert) then
338  write (this%iout, '(4x,a,i0)') 'NVERT = ', this%nvert
339  end if
340  !
341  write (this%iout, '(1x,a,/)') 'End Setting Discretization Dimensions'
342  !

◆ log_griddata()

subroutine disvmodule::log_griddata ( class(disvtype this,
type(disvfoundtype), intent(in)  found 
)
private

Definition at line 367 of file Disv.f90.

368  ! -- dummy
369  class(DisvType) :: this
370  type(DisvFoundType), intent(in) :: found
371  !
372  write (this%iout, '(1x,a)') 'Setting Discretization Griddata'
373  !
374  if (found%top) then
375  write (this%iout, '(4x,a)') 'TOP set from input file'
376  end if
377  !
378  if (found%botm) then
379  write (this%iout, '(4x,a)') 'BOTM set from input file'
380  end if
381  !
382  if (found%idomain) then
383  write (this%iout, '(4x,a)') 'IDOMAIN set from input file'
384  end if
385  !
386  write (this%iout, '(1x,a,/)') 'End Setting Discretization Griddata'
387  !

◆ log_options()

subroutine disvmodule::log_options ( class(disvtype this,
type(disvfoundtype), intent(in)  found 
)
private

Definition at line 228 of file Disv.f90.

229  ! -- dummy
230  class(DisvType) :: this
231  type(DisvFoundType), intent(in) :: found
232  !
233  write (this%iout, '(1x,a)') 'Setting Discretization Options'
234  !
235  if (found%length_units) then
236  write (this%iout, '(4x,a,i0)') 'Model length unit [0=UND, 1=FEET, &
237  &2=METERS, 3=CENTIMETERS] set as ', this%lenuni
238  end if
239  !
240  if (found%nogrb) then
241  write (this%iout, '(4x,a,i0)') 'Binary grid file [0=GRB, 1=NOGRB] &
242  &set as ', this%nogrb
243  end if
244  !
245  if (found%xorigin) then
246  write (this%iout, '(4x,a,G0)') 'XORIGIN = ', this%xorigin
247  end if
248  !
249  if (found%yorigin) then
250  write (this%iout, '(4x,a,G0)') 'YORIGIN = ', this%yorigin
251  end if
252  !
253  if (found%angrot) then
254  write (this%iout, '(4x,a,G0)') 'ANGROT = ', this%angrot
255  end if
256  !
257  write (this%iout, '(1x,a,/)') 'End Setting Discretization Options'
258  !

◆ nlarray_to_nodelist()

subroutine disvmodule::nlarray_to_nodelist ( class(disvtype 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 1787 of file Disv.f90.

1788  ! -- dummy
1789  class(DisvType) :: this
1790  integer(I4B), intent(in) :: maxbnd
1791  integer(I4B), dimension(:), pointer, contiguous :: darray
1792  integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
1793  integer(I4B), intent(inout) :: nbound
1794  character(len=*), intent(in) :: aname
1795  ! -- local
1796  integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr
1797  !
1798  ! -- set variables
1799  nlay = this%mshape(1)
1800  nrow = 1
1801  ncol = this%mshape(2)
1802  !
1803  nval = ncol * nrow
1804  !
1805  ! -- Copy array into nodelist
1806  ipos = 1
1807  ierr = 0
1808  do ir = 1, nrow
1809  do ic = 1, ncol
1810  nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
1811  il = darray(nodeu)
1812  if (il < 1 .or. il > nlay) then
1813  write (errmsg, '(a,i0,a)') &
1814  'Invalid layer number (', il, ').'
1815  call store_error(errmsg, terminate=.true.)
1816  end if
1817  nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
1818  noder = this%get_nodenumber(nodeu, 0)
1819  if (ipos > maxbnd) then
1820  ierr = ipos
1821  else
1822  nodelist(ipos) = noder
1823  end if
1824  ipos = ipos + 1
1825  end do
1826  end do
1827  !
1828  ! -- Check for errors
1829  nbound = ipos - 1
1830  if (ierr > 0) then
1831  write (errmsg, '(a,i0,a)') &
1832  'MAXBOUND dimension is too small. Increase MAXBOUND to ', ierr, '.'
1833  call store_error(errmsg, terminate=.true.)
1834  end if
1835  !
1836  ! -- If nbound < maxbnd, then initialize nodelist to zero in this range
1837  if (nbound < maxbnd) then
1838  do ipos = nbound + 1, maxbnd
1839  nodelist(ipos) = 0
1840  end do
1841  end if
1842  !
Here is the call graph for this function:

◆ nodeu_from_cellid()

integer(i4b) function disvmodule::nodeu_from_cellid ( class(disvtype 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 1325 of file Disv.f90.

1327  ! -- return
1328  integer(I4B) :: nodeu
1329  ! -- dummy
1330  class(DisvType) :: this
1331  character(len=*), intent(inout) :: cellid
1332  integer(I4B), intent(in) :: inunit
1333  integer(I4B), intent(in) :: iout
1334  logical, optional, intent(in) :: flag_string
1335  logical, optional, intent(in) :: allow_zero
1336  ! -- local
1337  integer(I4B) :: j, k, nlay, nrow, ncpl
1338  integer(I4B) :: lloclocal, ndum, istat, n
1339  integer(I4B) :: istart, istop
1340  real(DP) :: r
1341  !
1342  if (present(flag_string)) then
1343  if (flag_string) then
1344  ! Check to see if first token in cellid can be read as an integer.
1345  lloclocal = 1
1346  call urword(cellid, lloclocal, istart, istop, 1, ndum, r, iout, inunit)
1347  read (cellid(istart:istop), *, iostat=istat) n
1348  if (istat /= 0) then
1349  ! First token in cellid is not an integer; return flag to this effect.
1350  nodeu = -2
1351  return
1352  end if
1353  end if
1354  end if
1355  !
1356  nlay = this%mshape(1)
1357  nrow = 1
1358  ncpl = this%mshape(2)
1359  !
1360  lloclocal = 1
1361  call urword(cellid, lloclocal, istart, istop, 2, k, r, iout, inunit)
1362  call urword(cellid, lloclocal, istart, istop, 2, j, r, iout, inunit)
1363  !
1364  if (k == 0 .and. j == 0) then
1365  if (present(allow_zero)) then
1366  if (allow_zero) then
1367  nodeu = 0
1368  return
1369  end if
1370  end if
1371  end if
1372  !
1373  errmsg = ''
1374  !
1375  if (k < 1 .or. k > nlay) then
1376  write (errmsg, '(a,i0,a)') &
1377  'Layer number in list (', k, ') is outside of the grid.'
1378  end if
1379  if (j < 1 .or. j > ncpl) then
1380  write (errmsg, '(a,1x,a,i0,a)') &
1381  trim(adjustl(errmsg)), 'Cell2d number in list (', j, &
1382  ') is outside of the grid.'
1383  end if
1384  !
1385  nodeu = get_node(k, 1, j, nlay, nrow, ncpl)
1386  !
1387  if (nodeu < 1 .or. nodeu > this%nodesuser) then
1388  write (errmsg, '(a,1x,a,i0,a)') &
1389  trim(adjustl(errmsg)), &
1390  "Cell number cannot be determined for cellid ("// &
1391  trim(adjustl(cellid))//") and results in a user "// &
1392  "node number (", nodeu, ") that is outside of the grid."
1393  end if
1394  !
1395  if (len_trim(adjustl(errmsg)) > 0) then
1396  call store_error(errmsg)
1397  call store_error_unit(inunit)
1398  end if
1399  !
Here is the call graph for this function:

◆ nodeu_from_string()

integer(i4b) function disvmodule::nodeu_from_string ( class(disvtype 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 and within-layer cell number 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 1240 of file Disv.f90.

1242  ! -- dummy
1243  class(DisvType) :: this
1244  integer(I4B), intent(inout) :: lloc
1245  integer(I4B), intent(inout) :: istart
1246  integer(I4B), intent(inout) :: istop
1247  integer(I4B), intent(in) :: in
1248  integer(I4B), intent(in) :: iout
1249  character(len=*), intent(inout) :: line
1250  logical, optional, intent(in) :: flag_string
1251  logical, optional, intent(in) :: allow_zero
1252  integer(I4B) :: nodeu
1253  ! -- local
1254  integer(I4B) :: j, k, nlay, nrow, ncpl
1255  integer(I4B) :: lloclocal, ndum, istat, n
1256  real(DP) :: r
1257  !
1258  if (present(flag_string)) then
1259  if (flag_string) then
1260  ! Check to see if first token in line can be read as an integer.
1261  lloclocal = lloc
1262  call urword(line, lloclocal, istart, istop, 1, ndum, r, iout, in)
1263  read (line(istart:istop), *, iostat=istat) n
1264  if (istat /= 0) then
1265  ! First token in line is not an integer; return flag to this effect.
1266  nodeu = -2
1267  return
1268  end if
1269  end if
1270  end if
1271  !
1272  nlay = this%mshape(1)
1273  nrow = 1
1274  ncpl = this%mshape(2)
1275  !
1276  call urword(line, lloc, istart, istop, 2, k, r, iout, in)
1277  call urword(line, lloc, istart, istop, 2, j, r, iout, in)
1278  !
1279  if (k == 0 .and. j == 0) then
1280  if (present(allow_zero)) then
1281  if (allow_zero) then
1282  nodeu = 0
1283  return
1284  end if
1285  end if
1286  end if
1287  !
1288  errmsg = ''
1289  !
1290  if (k < 1 .or. k > nlay) then
1291  write (errmsg, '(a,i0,a)') &
1292  'Layer number in list (', k, ') is outside of the grid.'
1293  end if
1294  if (j < 1 .or. j > ncpl) then
1295  write (errmsg, '(a,1x,a,i0,a)') &
1296  trim(adjustl(errmsg)), 'Cell2d number in list (', j, &
1297  ') is outside of the grid.'
1298  end if
1299  !
1300  nodeu = get_node(k, 1, j, nlay, nrow, ncpl)
1301  !
1302  if (nodeu < 1 .or. nodeu > this%nodesuser) then
1303  write (errmsg, '(a,1x,a,i0,a)') &
1304  trim(adjustl(errmsg)), &
1305  "Node number in list (", nodeu, ") is outside of the grid. "// &
1306  "Cell number cannot be determined in line '"// &
1307  trim(adjustl(line))//"'."
1308  end if
1309  !
1310  if (len_trim(adjustl(errmsg)) > 0) then
1311  call store_error(errmsg)
1312  call store_error_unit(in)
1313  end if
1314  !
Here is the call graph for this function:

◆ nodeu_to_array()

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

Definition at line 891 of file Disv.f90.

892  ! -- dummy
893  class(DisvType) :: this
894  integer(I4B), intent(in) :: nodeu
895  integer(I4B), dimension(:), intent(inout) :: arr
896  ! -- local
897  integer(I4B) :: isize
898  integer(I4B) :: i, j, k
899  !
900  ! -- check the size of arr
901  isize = size(arr)
902  if (isize /= this%ndim) then
903  write (errmsg, '(a,i0,a,i0,a)') &
904  'Program error: nodeu_to_array size of array (', isize, &
905  ') is not equal to the discretization dimension (', this%ndim, ').'
906  call store_error(errmsg, terminate=.true.)
907  end if
908  !
909  ! -- get k, i, j
910  call get_ijk(nodeu, 1, this%ncpl, this%nlay, i, j, k)
911  !
912  ! -- fill array
913  arr(1) = k
914  arr(2) = j
915  !
Here is the call graph for this function:

◆ nodeu_to_string()

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

Definition at line 872 of file Disv.f90.

873  ! -- dummy
874  class(DisvType) :: this
875  integer(I4B), intent(in) :: nodeu
876  character(len=*), intent(inout) :: str
877  ! -- local
878  integer(I4B) :: i, j, k
879  character(len=10) :: kstr, jstr
880  !
881  call get_ijk(nodeu, 1, this%ncpl, this%nlay, i, j, k)
882  write (kstr, '(i10)') k
883  write (jstr, '(i10)') j
884  str = '('//trim(adjustl(kstr))//','// &
885  trim(adjustl(jstr))//')'
886  !
Here is the call graph for this function:

◆ read_dbl_array()

subroutine disvmodule::read_dbl_array ( class(disvtype), 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 1562 of file Disv.f90.

1564  ! -- dummy
1565  class(DisvType), intent(inout) :: this
1566  character(len=*), intent(inout) :: line
1567  integer(I4B), intent(inout) :: lloc
1568  integer(I4B), intent(inout) :: istart
1569  integer(I4B), intent(inout) :: istop
1570  integer(I4B), intent(in) :: in
1571  integer(I4B), intent(in) :: iout
1572  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray
1573  character(len=*), intent(in) :: aname
1574  ! -- local
1575  integer(I4B) :: ival
1576  real(DP) :: rval
1577  integer(I4B) :: nlay
1578  integer(I4B) :: nrow
1579  integer(I4B) :: ncol
1580  integer(I4B) :: nval
1581  real(DP), dimension(:), pointer, contiguous :: dtemp
1582  !
1583  ! -- Point the temporary pointer array, which is passed to the reading
1584  ! subroutine. The temporary array will point to dbuff if it is a
1585  ! reduced structured system, or to darray if it is an unstructured
1586  ! model.
1587  nlay = this%mshape(1)
1588  nrow = 1
1589  ncol = this%mshape(2)
1590  !
1591  if (this%nodes < this%nodesuser) then
1592  nval = this%nodesuser
1593  dtemp => this%dbuff
1594  else
1595  nval = this%nodes
1596  dtemp => darray
1597  end if
1598  !
1599  ! -- Read the array
1600  call urword(line, lloc, istart, istop, 1, ival, rval, iout, in)
1601  if (line(istart:istop) .EQ. 'LAYERED') then
1602  !
1603  ! -- Read structured input
1604  call readarray(in, dtemp, aname, this%ndim, ncol, nrow, nlay, nval, &
1605  iout, 1, nlay)
1606  else
1607  !
1608  ! -- Read unstructured input
1609  call readarray(in, dtemp, aname, this%ndim, nval, iout, 0)
1610  end if
1611  !
1612  ! -- If reduced model, then need to copy from dtemp(=>dbuff) to darray
1613  if (this%nodes < this%nodesuser) then
1614  call this%fill_grid_array(dtemp, darray)
1615  end if
1616  !
Here is the call graph for this function:

◆ read_int_array()

subroutine disvmodule::read_int_array ( class(disvtype), 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 1503 of file Disv.f90.

1505  ! -- dummy
1506  class(DisvType), intent(inout) :: this
1507  character(len=*), intent(inout) :: line
1508  integer(I4B), intent(inout) :: lloc
1509  integer(I4B), intent(inout) :: istart
1510  integer(I4B), intent(inout) :: istop
1511  integer(I4B), intent(in) :: in
1512  integer(I4B), intent(in) :: iout
1513  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: iarray
1514  character(len=*), intent(in) :: aname
1515  ! -- local
1516  integer(I4B) :: ival
1517  real(DP) :: rval
1518  integer(I4B) :: nlay
1519  integer(I4B) :: nrow
1520  integer(I4B) :: ncol
1521  integer(I4B) :: nval
1522  integer(I4B), dimension(:), pointer, contiguous :: itemp
1523  !
1524  ! -- Point the temporary pointer array, which is passed to the reading
1525  ! subroutine. The temporary array will point to ibuff if it is a
1526  ! reduced structured system, or to iarray if it is an unstructured
1527  ! model.
1528  nlay = this%mshape(1)
1529  nrow = 1
1530  ncol = this%mshape(2)
1531  !
1532  if (this%nodes < this%nodesuser) then
1533  nval = this%nodesuser
1534  itemp => this%ibuff
1535  else
1536  nval = this%nodes
1537  itemp => iarray
1538  end if
1539  !
1540  ! -- Read the array
1541  call urword(line, lloc, istart, istop, 1, ival, rval, iout, in)
1542  if (line(istart:istop) .EQ. 'LAYERED') then
1543  !
1544  ! -- Read layered input
1545  call readarray(in, itemp, aname, this%ndim, ncol, nrow, nlay, nval, &
1546  iout, 1, nlay)
1547  else
1548  !
1549  ! -- Read unstructured input
1550  call readarray(in, itemp, aname, this%ndim, nval, iout, 0)
1551  end if
1552  !
1553  ! -- If reduced model, then need to copy from itemp(=>ibuff) to iarray
1554  if (this%nodes < this%nodesuser) then
1555  call this%fill_grid_array(itemp, iarray)
1556  end if
1557  !
Here is the call graph for this function:

◆ read_layer_array()

subroutine disvmodule::read_layer_array ( class(disvtype 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 1624 of file Disv.f90.

1626  ! -- dummy
1627  class(DisvType) :: this
1628  integer(I4B), intent(in) :: ncolbnd
1629  integer(I4B), intent(in) :: maxbnd
1630  integer(I4B), dimension(maxbnd) :: nodelist
1631  real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray
1632  integer(I4B), intent(in) :: icolbnd
1633  character(len=*), intent(in) :: aname
1634  integer(I4B), intent(in) :: inunit
1635  integer(I4B), intent(in) :: iout
1636  ! -- local
1637  integer(I4B) :: ir, ic, ncol, nrow, nlay, nval, ipos, nodeu
1638  !
1639  ! -- set variables
1640  nlay = this%mshape(1)
1641  nrow = 1
1642  ncol = this%mshape(2)
1643  !
1644  ! -- Read the array
1645  nval = ncol * nrow
1646  call readarray(inunit, this%dbuff, aname, this%ndim, nval, iout, 0)
1647  !
1648  ! -- Copy array into bound. Note that this routine was substantially
1649  ! changed on 9/21/2021 to support changes to READASARRAYS input
1650  ! for recharge and evapotranspiration. nodelist and bound are of
1651  ! size nrow * ncol and correspond directly to dbuff.
1652  ipos = 1
1653  do ir = 1, nrow
1654  do ic = 1, ncol
1655  nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
1656  darray(icolbnd, ipos) = this%dbuff(nodeu)
1657  ipos = ipos + 1
1658  end do
1659  end do
1660  !
Here is the call graph for this function:

◆ record_array()

subroutine disvmodule::record_array ( class(disvtype), 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 1668 of file Disv.f90.

1670  ! -- dummy
1671  class(DisvType), intent(inout) :: this
1672  real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record
1673  integer(I4B), intent(in) :: iout !< ascii output unit number
1674  integer(I4B), intent(in) :: iprint !< whether to print the array
1675  integer(I4B), intent(in) :: idataun !< binary output unit number, if negative don't write by layers, write entire array
1676  character(len=*), intent(in) :: aname !< text descriptor
1677  character(len=*), intent(in) :: cdatafmp !< write format
1678  integer(I4B), intent(in) :: nvaluesp !< values per line
1679  integer(I4B), intent(in) :: nwidthp !< number width
1680  character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E)
1681  real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain
1682  ! -- local
1683  integer(I4B) :: k, ifirst
1684  integer(I4B) :: nlay
1685  integer(I4B) :: nrow
1686  integer(I4B) :: ncol
1687  integer(I4B) :: nval
1688  integer(I4B) :: nodeu, noder
1689  integer(I4B) :: istart, istop
1690  real(DP), dimension(:), pointer, contiguous :: dtemp
1691  ! -- formats
1692  character(len=*), parameter :: fmthsv = &
1693  "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, &
1694  &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)"
1695  !
1696  ! -- set variables
1697  nlay = this%mshape(1)
1698  nrow = 1
1699  ncol = this%mshape(2)
1700  !
1701  ! -- If this is a reduced model, then copy the values from darray into
1702  ! dtemp.
1703  if (this%nodes < this%nodesuser) then
1704  nval = this%nodes
1705  dtemp => this%dbuff
1706  do nodeu = 1, this%nodesuser
1707  noder = this%get_nodenumber(nodeu, 0)
1708  if (noder <= 0) then
1709  dtemp(nodeu) = dinact
1710  cycle
1711  end if
1712  dtemp(nodeu) = darray(noder)
1713  end do
1714  else
1715  nval = this%nodes
1716  dtemp => darray
1717  end if
1718  !
1719  ! -- Print to iout if iprint /= 0
1720  if (iprint /= 0) then
1721  istart = 1
1722  do k = 1, nlay
1723  istop = istart + nrow * ncol - 1
1724  call ulaprufw(ncol, nrow, kstp, kper, k, iout, dtemp(istart:istop), &
1725  aname, cdatafmp, nvaluesp, nwidthp, editdesc)
1726  istart = istop + 1
1727  end do
1728  end if
1729  !
1730  ! -- Save array to an external file.
1731  if (idataun > 0) then
1732  ! -- write to binary file by layer
1733  ifirst = 1
1734  istart = 1
1735  do k = 1, nlay
1736  istop = istart + nrow * ncol - 1
1737  if (ifirst == 1) write (iout, fmthsv) &
1738  trim(adjustl(aname)), idataun, &
1739  kstp, kper
1740  ifirst = 0
1741  call ulasav(dtemp(istart:istop), aname, kstp, kper, &
1742  pertim, totim, ncol, nrow, k, idataun)
1743  istart = istop + 1
1744  end do
1745  elseif (idataun < 0) then
1746  !
1747  ! -- write entire array as one record
1748  call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, &
1749  iout, delt, pertim, totim)
1750  end if
1751  !
Here is the call graph for this function:

◆ record_srcdst_list_header()

subroutine disvmodule::record_srcdst_list_header ( class(disvtype 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 1756 of file Disv.f90.

1759  ! -- dummy
1760  class(DisvType) :: this
1761  character(len=16), intent(in) :: text
1762  character(len=16), intent(in) :: textmodel
1763  character(len=16), intent(in) :: textpackage
1764  character(len=16), intent(in) :: dstmodel
1765  character(len=16), intent(in) :: dstpackage
1766  integer(I4B), intent(in) :: naux
1767  character(len=16), dimension(:), intent(in) :: auxtxt
1768  integer(I4B), intent(in) :: ibdchn
1769  integer(I4B), intent(in) :: nlist
1770  integer(I4B), intent(in) :: iout
1771  ! -- local
1772  integer(I4B) :: nlay, nrow, ncol
1773  !
1774  nlay = this%mshape(1)
1775  nrow = 1
1776  ncol = this%mshape(2)
1777  !
1778  ! -- Use ubdsv06 to write list header
1779  call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, &
1780  ibdchn, naux, auxtxt, ncol, nrow, nlay, &
1781  nlist, iout, delt, pertim, totim)
1782  !
Here is the call graph for this function:

◆ source_cell2d()

subroutine disvmodule::source_cell2d ( class(disvtype this)

Definition at line 591 of file Disv.f90.

592  ! -- dummy
593  class(DisvType) :: this
594  ! -- locals
595  integer(I4B), dimension(:), contiguous, pointer :: icell2d => null()
596  integer(I4B), dimension(:), contiguous, pointer :: ncvert => null()
597  integer(I4B), dimension(:), contiguous, pointer :: icvert => null()
598  real(DP), dimension(:), contiguous, pointer :: cell_x => null()
599  real(DP), dimension(:), contiguous, pointer :: cell_y => null()
600  integer(I4B) :: i
601  !
602  ! -- set pointers to input path ncvert and icvert
603  call mem_setptr(icell2d, 'ICELL2D', this%input_mempath)
604  call mem_setptr(ncvert, 'NCVERT', this%input_mempath)
605  call mem_setptr(icvert, 'ICVERT', this%input_mempath)
606  !
607  ! --
608  if (associated(icell2d) .and. associated(ncvert) &
609  .and. associated(icvert)) then
610  call this%define_cellverts(icell2d, ncvert, icvert)
611  else
612  call store_error('Required cell vertex array(s) [ICELL2D, NCVERT, ICVERT] &
613  &not found.')
614  end if
615  !
616  ! -- copy cell center idm sourced values to local arrays
617  call mem_setptr(cell_x, 'XC', this%input_mempath)
618  call mem_setptr(cell_y, 'YC', this%input_mempath)
619  !
620  ! -- set cell centers
621  if (associated(cell_x) .and. associated(cell_y)) then
622  do i = 1, this%ncpl
623  this%cellxy(1, i) = cell_x(i)
624  this%cellxy(2, i) = cell_y(i)
625  end do
626  else
627  call store_error('Required cell center arrays not found.')
628  end if
629  !
630  ! -- log
631  if (this%iout > 0) then
632  write (this%iout, '(1x,a)') 'Discretization Cell2d data loaded'
633  end if
634  !
Here is the call graph for this function:

◆ source_dimensions()

subroutine disvmodule::source_dimensions ( class(disvtype this)
private

Definition at line 263 of file Disv.f90.

264  ! -- dummy
265  class(DisvType) :: this
266  ! -- locals
267  integer(I4B) :: j, k
268  type(DisvFoundType) :: found
269  !
270  ! -- update defaults with idm sourced values
271  call mem_set_value(this%nlay, 'NLAY', this%input_mempath, found%nlay)
272  call mem_set_value(this%ncpl, 'NCPL', this%input_mempath, found%ncpl)
273  call mem_set_value(this%nvert, 'NVERT', this%input_mempath, found%nvert)
274  !
275  ! -- log simulation values
276  if (this%iout > 0) then
277  call this%log_dimensions(found)
278  end if
279  !
280  ! -- verify dimensions were set
281  if (this%nlay < 1) then
282  call store_error( &
283  'NLAY was not specified or was specified incorrectly.')
284  call store_error_filename(this%input_fname)
285  end if
286  if (this%ncpl < 1) then
287  call store_error( &
288  'NCPL was not specified or was specified incorrectly.')
289  call store_error_filename(this%input_fname)
290  end if
291  if (this%nvert < 1) then
292  call store_error( &
293  'NVERT was not specified or was specified incorrectly.')
294  call store_error_filename(this%input_fname)
295  end if
296  !
297  ! -- Calculate nodesuser
298  this%nodesuser = this%nlay * this%ncpl
299  !
300  ! -- Allocate non-reduced vectors for disv
301  call mem_allocate(this%idomain, this%ncpl, this%nlay, 'IDOMAIN', &
302  this%memoryPath)
303  call mem_allocate(this%top1d, this%ncpl, 'TOP1D', this%memoryPath)
304  call mem_allocate(this%bot2d, this%ncpl, this%nlay, 'BOT2D', &
305  this%memoryPath)
306  !
307  ! -- Allocate vertices array
308  call mem_allocate(this%vertices, 2, this%nvert, 'VERTICES', this%memoryPath)
309  call mem_allocate(this%cellxy, 2, this%ncpl, 'CELLXY', this%memoryPath)
310  !
311  ! -- initialize all cells to be active (idomain = 1)
312  do k = 1, this%nlay
313  do j = 1, this%ncpl
314  this%idomain(j, k) = 1
315  end do
316  end do
317  !
Here is the call graph for this function:

◆ source_griddata()

subroutine disvmodule::source_griddata ( class(disvtype this)
private

Definition at line 347 of file Disv.f90.

348  ! -- dummy
349  class(DisvType) :: this
350  ! -- locals
351  type(DisvFoundType) :: found
352  !
353  ! -- update defaults with idm sourced values
354  call mem_set_value(this%top1d, 'TOP', this%input_mempath, found%top)
355  call mem_set_value(this%bot2d, 'BOTM', this%input_mempath, found%botm)
356  call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
357  !
358  ! -- log simulation values
359  if (this%iout > 0) then
360  call this%log_griddata(found)
361  end if
362  !

◆ source_options()

subroutine disvmodule::source_options ( class(disvtype this)
private

Definition at line 203 of file Disv.f90.

204  ! -- dummy
205  class(DisvType) :: this
206  ! -- locals
207  character(len=LENVARNAME), dimension(3) :: lenunits = &
208  &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS']
209  type(disvfoundtype) :: found
210  !
211  ! -- update defaults with idm sourced values
212  call mem_set_value(this%lenuni, 'LENGTH_UNITS', this%input_mempath, &
213  lenunits, found%length_units)
214  call mem_set_value(this%nogrb, 'NOGRB', this%input_mempath, found%nogrb)
215  call mem_set_value(this%xorigin, 'XORIGIN', this%input_mempath, found%xorigin)
216  call mem_set_value(this%yorigin, 'YORIGIN', this%input_mempath, found%yorigin)
217  call mem_set_value(this%angrot, 'ANGROT', this%input_mempath, found%angrot)
218  !
219  ! -- log values to list file
220  if (this%iout > 0) then
221  call this%log_options(found)
222  end if
223  !

◆ source_vertices()

subroutine disvmodule::source_vertices ( class(disvtype this)
private

Definition at line 519 of file Disv.f90.

520  ! -- dummy
521  class(DisvType) :: this
522  ! -- local
523  integer(I4B) :: i
524  real(DP), dimension(:), contiguous, pointer :: vert_x => null()
525  real(DP), dimension(:), contiguous, pointer :: vert_y => null()
526  !
527  ! -- set pointers to memory manager input arrays
528  call mem_setptr(vert_x, 'XV', this%input_mempath)
529  call mem_setptr(vert_y, 'YV', this%input_mempath)
530  !
531  ! -- set vertices 2d array
532  if (associated(vert_x) .and. associated(vert_y)) then
533  do i = 1, this%nvert
534  this%vertices(1, i) = vert_x(i)
535  this%vertices(2, i) = vert_y(i)
536  end do
537  else
538  call store_error('Required Vertex arrays not found.')
539  end if
540  !
541  ! -- log
542  if (this%iout > 0) then
543  write (this%iout, '(1x,a)') 'Discretization Vertex data loaded'
544  end if
545  !
Here is the call graph for this function:

◆ supports_layers()

logical function disvmodule::supports_layers ( class(disvtype this)
private

Definition at line 1404 of file Disv.f90.

1405  ! -- dummy
1406  class(DisvType) :: this
1407  !
1408  supports_layers = .true.
1409  !

◆ write_grb()

subroutine disvmodule::write_grb ( class(disvtype this,
integer(i4b), dimension(:), intent(in)  icelltype 
)
private

Definition at line 711 of file Disv.f90.

712  ! -- modules
713  use openspecmodule, only: access, form
714  use constantsmodule, only: lenbigline
715  ! -- dummy
716  class(DisvType) :: this
717  integer(I4B), dimension(:), intent(in) :: icelltype
718  ! -- local
719  integer(I4B) :: iunit, i, ntxt, version
720  integer(I4B), parameter :: lentxt = 100
721  character(len=50) :: txthdr
722  character(len=lentxt) :: txt
723  character(len=LINELENGTH) :: fname
724  character(len=LENBIGLINE) :: crs
725  logical(LGP) :: found_crs
726  ! -- formats
727  character(len=*), parameter :: fmtgrdsave = &
728  "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', &
729  &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)"
730  !
731  ! -- Initialize
732  version = 1
733  ntxt = 20
734  !
735  call mem_set_value(crs, 'CRS', this%input_mempath, found_crs)
736  !
737  ! -- set version
738  if (found_crs) then
739  ntxt = ntxt + 1
740  version = 2
741  end if
742  !
743  ! -- Open the file
744  fname = trim(this%output_fname)
745  iunit = getunit()
746  write (this%iout, fmtgrdsave) iunit, trim(adjustl(fname))
747  call openfile(iunit, this%iout, trim(adjustl(fname)), 'DATA(BINARY)', &
748  form, access, 'REPLACE')
749  !
750  ! -- write header information
751  write (txthdr, '(a)') 'GRID DISV'
752  txthdr(50:50) = new_line('a')
753  write (iunit) txthdr
754  write (txthdr, '(a, i0)') 'VERSION ', version
755  txthdr(50:50) = new_line('a')
756  write (iunit) txthdr
757  write (txthdr, '(a, i0)') 'NTXT ', ntxt
758  txthdr(50:50) = new_line('a')
759  write (iunit) txthdr
760  write (txthdr, '(a, i0)') 'LENTXT ', lentxt
761  txthdr(50:50) = new_line('a')
762  write (iunit) txthdr
763  !
764  ! -- write variable definitions
765  write (txt, '(3a, i0)') 'NCELLS ', 'INTEGER ', 'NDIM 0 # ', this%nodesuser
766  txt(lentxt:lentxt) = new_line('a')
767  write (iunit) txt
768  write (txt, '(3a, i0)') 'NLAY ', 'INTEGER ', 'NDIM 0 # ', this%nlay
769  txt(lentxt:lentxt) = new_line('a')
770  write (iunit) txt
771  write (txt, '(3a, i0)') 'NCPL ', 'INTEGER ', 'NDIM 0 # ', this%ncpl
772  txt(lentxt:lentxt) = new_line('a')
773  write (iunit) txt
774  write (txt, '(3a, i0)') 'NVERT ', 'INTEGER ', 'NDIM 0 # ', this%nvert
775  txt(lentxt:lentxt) = new_line('a')
776  write (iunit) txt
777  write (txt, '(3a, i0)') 'NJAVERT ', 'INTEGER ', 'NDIM 0 # ', size(this%javert)
778  txt(lentxt:lentxt) = new_line('a')
779  write (iunit) txt
780  write (txt, '(3a, i0)') 'NJA ', 'INTEGER ', 'NDIM 0 # ', this%con%nja
781  txt(lentxt:lentxt) = new_line('a')
782  write (iunit) txt
783  write (txt, '(3a, 1pg25.15e3)') &
784  'XORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%xorigin
785  txt(lentxt:lentxt) = new_line('a')
786  write (iunit) txt
787  write (txt, '(3a, 1pg25.15e3)') &
788  'YORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%yorigin
789  txt(lentxt:lentxt) = new_line('a')
790  write (iunit) txt
791  write (txt, '(3a, 1pg25.15e3)') 'ANGROT ', 'DOUBLE ', 'NDIM 0 # ', this%angrot
792  txt(lentxt:lentxt) = new_line('a')
793  write (iunit) txt
794  write (txt, '(3a, i0)') 'TOP ', 'DOUBLE ', 'NDIM 1 ', this%ncpl
795  txt(lentxt:lentxt) = new_line('a')
796  write (iunit) txt
797  write (txt, '(3a, i0)') 'BOTM ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser
798  txt(lentxt:lentxt) = new_line('a')
799  write (iunit) txt
800  write (txt, '(3a, i0)') 'VERTICES ', 'DOUBLE ', 'NDIM 2 2 ', this%nvert
801  txt(lentxt:lentxt) = new_line('a')
802  write (iunit) txt
803  write (txt, '(3a, i0)') 'CELLX ', 'DOUBLE ', 'NDIM 1 ', this%ncpl
804  txt(lentxt:lentxt) = new_line('a')
805  write (iunit) txt
806  write (txt, '(3a, i0)') 'CELLY ', 'DOUBLE ', 'NDIM 1 ', this%ncpl
807  txt(lentxt:lentxt) = new_line('a')
808  write (iunit) txt
809  write (txt, '(3a, i0)') 'IAVERT ', 'INTEGER ', 'NDIM 1 ', this%ncpl + 1
810  txt(lentxt:lentxt) = new_line('a')
811  write (iunit) txt
812  write (txt, '(3a, i0)') 'JAVERT ', 'INTEGER ', 'NDIM 1 ', size(this%javert)
813  txt(lentxt:lentxt) = new_line('a')
814  write (iunit) txt
815  write (txt, '(3a, i0)') 'IA ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1
816  txt(lentxt:lentxt) = new_line('a')
817  write (iunit) txt
818  write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', size(this%con%jausr)
819  txt(lentxt:lentxt) = new_line('a')
820  write (iunit) txt
821  write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
822  txt(lentxt:lentxt) = new_line('a')
823  write (iunit) txt
824  write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
825  txt(lentxt:lentxt) = new_line('a')
826  write (iunit) txt
827  !
828  ! -- if version 2 write character array headers
829  if (version == 2) then
830  if (found_crs) then
831  write (txt, '(3a, i0)') 'CRS ', 'CHARACTER ', 'NDIM 1 ', &
832  len_trim(crs)
833  txt(lentxt:lentxt) = new_line('a')
834  write (iunit) txt
835  end if
836  end if
837  !
838  ! -- write data
839  write (iunit) this%nodesuser ! ncells
840  write (iunit) this%nlay ! nlay
841  write (iunit) this%ncpl ! ncpl
842  write (iunit) this%nvert ! nvert
843  write (iunit) size(this%javert) ! njavert
844  write (iunit) this%nja ! nja
845  write (iunit) this%xorigin ! xorigin
846  write (iunit) this%yorigin ! yorigin
847  write (iunit) this%angrot ! angrot
848  write (iunit) this%top1d ! top
849  write (iunit) this%bot2d ! botm
850  write (iunit) this%vertices ! vertices
851  write (iunit) (this%cellxy(1, i), i=1, this%ncpl) ! cellx
852  write (iunit) (this%cellxy(2, i), i=1, this%ncpl) ! celly
853  write (iunit) this%iavert ! iavert
854  write (iunit) this%javert ! javert
855  write (iunit) this%con%iausr ! iausr
856  write (iunit) this%con%jausr ! jausr
857  write (iunit) this%idomain ! idomain
858  write (iunit) icelltype ! icelltype
859  !
860  ! -- if version 2 write character array data
861  if (version == 2) then
862  if (found_crs) write (iunit) trim(crs) ! crs user input
863  end if
864  !
865  ! -- Close the file
866  close (iunit)
867  !
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: