From b8bdf6106fdae5c290389fdfbbac30c117b03b67 Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Mon, 15 Apr 2024 14:40:12 -0500 Subject: [PATCH] feat(disv2d): add DISV2D grid support for SWF overland flow (#1726) * feat(disv2d): introduce limited 2D vertex grid (for overland flow) * needs testing * add io guide docs for disv2d * msvs proj and makefile updates * remove unnecessary files --- doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn | 9 + .../dfn/{swf-disv.dfn => swf-disv2d.dfn} | 53 +- .../mf6ivar/examples/swf-disv2d-example.dat | 31 + doc/mf6io/mf6ivar/md/mf6ivar.md | 41 +- doc/mf6io/mf6ivar/mf6ivar.py | 2 +- doc/mf6io/mf6ivar/tex/appendixA.tex | 10 +- doc/mf6io/mf6ivar/tex/swf-dfw-desc.tex | 2 + doc/mf6io/mf6ivar/tex/swf-dfw-options.dat | 1 + doc/mf6io/mf6ivar/tex/swf-disl-cell2d.dat | 5 - doc/mf6io/mf6ivar/tex/swf-disv-griddata.dat | 8 - doc/mf6io/mf6ivar/tex/swf-disv1d-desc.tex | 2 + doc/mf6io/mf6ivar/tex/swf-disv1d-options.dat | 1 + ...-disv-cell2d.dat => swf-disv2d-cell2d.dat} | 0 ...{swf-disv-desc.tex => swf-disv2d-desc.tex} | 16 +- ...mensions.dat => swf-disv2d-dimensions.dat} | 3 +- doc/mf6io/mf6ivar/tex/swf-disv2d-griddata.dat | 6 + ...isv-options.dat => swf-disv2d-options.dat} | 0 ...v-vertices.dat => swf-disv2d-vertices.dat} | 0 doc/mf6io/swf/disv2d.tex | 19 + doc/mf6io/swf/namefile.tex | 1 + doc/mf6io/swf/swf.tex | 4 + make/makefile | 3 +- msvs/mf6core.vfproj | 5 +- src/Idm/selector/IdmSwfDfnSelector.f90 | 20 +- src/Idm/swf-disv1didm.f90 | 19 + .../{swf-disvidm.f90 => swf-disv2didm.f90} | 206 +-- src/Model/Discretization/Disv2d.f90 | 1387 +++++++++++++++++ src/Model/SurfaceWaterFlow/swf-cdb.f90 | 13 +- src/Model/SurfaceWaterFlow/swf-dfw.f90 | 3 + src/Model/SurfaceWaterFlow/swf.f90 | 8 +- src/Utilities/Idm/SourceCommon.f90 | 17 + .../Idm/mf6blockfile/LoadMf6File.f90 | 3 +- src/meson.build | 3 +- utils/idmloader/dfns.txt | 2 +- 34 files changed, 1666 insertions(+), 237 deletions(-) rename doc/mf6io/mf6ivar/dfn/{swf-disv.dfn => swf-disv2d.dfn} (66%) create mode 100644 doc/mf6io/mf6ivar/examples/swf-disv2d-example.dat delete mode 100644 doc/mf6io/mf6ivar/tex/swf-disl-cell2d.dat delete mode 100644 doc/mf6io/mf6ivar/tex/swf-disv-griddata.dat rename doc/mf6io/mf6ivar/tex/{swf-disv-cell2d.dat => swf-disv2d-cell2d.dat} (100%) rename doc/mf6io/mf6ivar/tex/{swf-disv-desc.tex => swf-disv2d-desc.tex} (57%) rename doc/mf6io/mf6ivar/tex/{swf-disv-dimensions.dat => swf-disv2d-dimensions.dat} (63%) create mode 100644 doc/mf6io/mf6ivar/tex/swf-disv2d-griddata.dat rename doc/mf6io/mf6ivar/tex/{swf-disv-options.dat => swf-disv2d-options.dat} (100%) rename doc/mf6io/mf6ivar/tex/{swf-disv-vertices.dat => swf-disv2d-vertices.dat} (100%) create mode 100644 doc/mf6io/swf/disv2d.tex rename src/Idm/{swf-disvidm.f90 => swf-disv2didm.f90} (70%) create mode 100644 src/Model/Discretization/Disv2d.f90 diff --git a/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn b/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn index 0ce04347669..84fcc63e939 100644 --- a/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn +++ b/doc/mf6io/mf6ivar/dfn/swf-disv1d.dfn @@ -40,6 +40,15 @@ optional true longname rotation angle description counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +block options +name export_array_ascii +type keyword +reader urword +optional true +mf6internal export_ascii +longname export array variables to layered ascii files. +description keyword that specifies input grid arrays, which already support the layered keyword, should be written to layered ascii output files. + # --------------------- swf disv1d dimensions --------------------- block dimensions diff --git a/doc/mf6io/mf6ivar/dfn/swf-disv.dfn b/doc/mf6io/mf6ivar/dfn/swf-disv2d.dfn similarity index 66% rename from doc/mf6io/mf6ivar/dfn/swf-disv.dfn rename to doc/mf6io/mf6ivar/dfn/swf-disv2d.dfn index 062fb01f169..e9b41acabc4 100644 --- a/doc/mf6io/mf6ivar/dfn/swf-disv.dfn +++ b/doc/mf6io/mf6ivar/dfn/swf-disv2d.dfn @@ -1,4 +1,4 @@ -# --------------------- gwt disv options --------------------- +# --------------------- swf disv2d options --------------------- block options name length_units @@ -21,16 +21,16 @@ name xorigin type double precision reader urword optional true -longname x-position origin of the model grid coordinate system -description x-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +longname x-position of the model grid origin +description x-position of the lower-left corner of the model grid. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. block options name yorigin type double precision reader urword optional true -longname y-position origin of the model grid coordinate system -description y-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +longname y-position of the model grid origin +description y-position of the lower-left corner of the model grid. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. block options name angrot @@ -38,7 +38,7 @@ type double precision reader urword optional true longname rotation angle -description counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +description counter-clockwise rotation angle (in degrees) of the lower-left corner of the model grid. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. block options name export_array_ascii @@ -49,18 +49,10 @@ mf6internal export_ascii longname export array variables to layered ascii files. description keyword that specifies input grid arrays, which already support the layered keyword, should be written to layered ascii output files. -# --------------------- gwt disv dimensions --------------------- +# --------------------- swf disv2d dimensions --------------------- block dimensions -name nlay -type integer -reader urword -optional false -longname number of layers -description is the number of layers in the model grid. - -block dimensions -name ncpl +name nodes type integer reader urword optional false @@ -75,37 +67,28 @@ optional false longname number of columns description is the total number of (x, y) vertex pairs used to characterize the horizontal configuration of the model grid. -# --------------------- gwt disv griddata --------------------- +# --------------------- swf disv2d griddata --------------------- block griddata -name top +name bottom type double precision -shape (ncpl) +shape (nodes) reader readarray -longname model top elevation -description is the top elevation for each cell in the top model layer. - -block griddata -name botm -type double precision -shape (ncpl, nlay) -reader readarray -layered true +layered false longname model bottom elevation description is the bottom elevation for each cell. block griddata name idomain type integer -shape (ncpl, nlay) +shape (nodes) reader readarray -layered true +layered false optional true longname idomain existence array -description is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. - +description is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1 or greater, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. -# --------------------- gwt disv vertices --------------------- +# --------------------- swf disv2d vertices --------------------- block vertices name vertices @@ -148,12 +131,12 @@ longname y-coordinate for vertex description is the y-coordinate for the vertex. -# --------------------- gwt disv cell2d --------------------- +# --------------------- swf disv cell2d --------------------- block cell2d name cell2d type recarray icell2d xc yc ncvert icvert -shape (ncpl) +shape (nodes) reader urword optional false longname cell2d data diff --git a/doc/mf6io/mf6ivar/examples/swf-disv2d-example.dat b/doc/mf6io/mf6ivar/examples/swf-disv2d-example.dat new file mode 100644 index 00000000000..6551a59ee78 --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/swf-disv2d-example.dat @@ -0,0 +1,31 @@ +BEGIN OPTIONS +END OPTIONS + +BEGIN DIMENSIONS + NODES 4 + NVERT 9 +END DIMENSIONS + +BEGIN GRIDDATA + bottom + CONSTANT 0.0 +END GRIDDATA + +BEGIN VERTICES + 1 0.0 0.0 + 2 1.0 0.0 + 3 2.0 0.0 + 4 0.0 1.0 + 5 1.0 1.0 + 6 2.0 1.0 + 7 0.0 2.0 + 8 1.0 2.0 + 9 2.0 2.0 +END VERTICES + +BEGIN CELL2D + 1 0.5 0.5 4 1 2 5 4 + 2 1.5 0.5 4 2 3 6 5 + 3 0.5 1.5 4 4 5 8 7 + 4 1.5 1.5 4 5 6 9 8 +END CELL2D \ No newline at end of file diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index 4ebf683b069..afc0fb00eec 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -1532,6 +1532,7 @@ | SWF | DISV1D | OPTIONS | XORIGIN | DOUBLE PRECISION | x-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | | SWF | DISV1D | OPTIONS | YORIGIN | DOUBLE PRECISION | y-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | | SWF | DISV1D | OPTIONS | ANGROT | DOUBLE PRECISION | counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | +| SWF | DISV1D | OPTIONS | EXPORT_ARRAY_ASCII | KEYWORD | keyword that specifies input grid arrays, which already support the layered keyword, should be written to layered ascii output files. | | SWF | DISV1D | DIMENSIONS | NODES | INTEGER | is the number of linear cells. | | SWF | DISV1D | DIMENSIONS | NVERT | INTEGER | is the total number of (x, y, z) vertex pairs used to characterize the model grid. | | SWF | DISV1D | GRIDDATA | LENGTH | DOUBLE PRECISION (NODES) | length for each one-dimensional cell | @@ -1557,34 +1558,34 @@ | SWF | DIS2D | GRIDDATA | DELC | DOUBLE PRECISION (NROW) | is the row spacing in the column direction. | | SWF | DIS2D | GRIDDATA | BOTM | DOUBLE PRECISION (NCOL, NROW) | is the bottom elevation for each cell. | | SWF | DIS2D | GRIDDATA | IDOMAIN | INTEGER (NCOL, NROW) | is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. | -| SWF | DISV | OPTIONS | LENGTH_UNITS | STRING | is the length units used for this model. Values can be ``FEET'', ``METERS'', or ``CENTIMETERS''. If not specified, the default is ``UNKNOWN''. | -| SWF | DISV | OPTIONS | NOGRB | KEYWORD | keyword to deactivate writing of the binary grid file. | -| SWF | DISV | OPTIONS | XORIGIN | DOUBLE PRECISION | x-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | -| SWF | DISV | OPTIONS | YORIGIN | DOUBLE PRECISION | y-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | -| SWF | DISV | OPTIONS | ANGROT | DOUBLE PRECISION | counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | -| SWF | DISV | OPTIONS | EXPORT_ARRAY_ASCII | KEYWORD | keyword that specifies input grid arrays, which already support the layered keyword, should be written to layered ascii output files. | -| SWF | DISV | DIMENSIONS | NLAY | INTEGER | is the number of layers in the model grid. | -| SWF | DISV | DIMENSIONS | NCPL | INTEGER | is the number of cells per layer. This is a constant value for the grid and it applies to all layers. | -| SWF | DISV | DIMENSIONS | NVERT | INTEGER | is the total number of (x, y) vertex pairs used to characterize the horizontal configuration of the model grid. | -| SWF | DISV | GRIDDATA | TOP | DOUBLE PRECISION (NCPL) | is the top elevation for each cell in the top model layer. | -| SWF | DISV | GRIDDATA | BOTM | DOUBLE PRECISION (NCPL, NLAY) | is the bottom elevation for each cell. | -| SWF | DISV | GRIDDATA | IDOMAIN | INTEGER (NCPL, NLAY) | is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. | -| SWF | DISV | VERTICES | IV | INTEGER | is the vertex number. Records in the VERTICES block must be listed in consecutive order from 1 to NVERT. | -| SWF | DISV | VERTICES | XV | DOUBLE PRECISION | is the x-coordinate for the vertex. | -| SWF | DISV | VERTICES | YV | DOUBLE PRECISION | is the y-coordinate for the vertex. | -| SWF | DISV | CELL2D | ICELL2D | INTEGER | is the CELL2D number. Records in the CELL2D block must be listed in consecutive order from the first to the last. | -| SWF | DISV | CELL2D | XC | DOUBLE PRECISION | is the x-coordinate for the cell center. | -| SWF | DISV | CELL2D | YC | DOUBLE PRECISION | is the y-coordinate for the cell center. | -| SWF | DISV | CELL2D | NCVERT | INTEGER | is the number of vertices required to define the cell. There may be a different number of vertices for each cell. | -| SWF | DISV | CELL2D | ICVERT | INTEGER (NCVERT) | is an array of integer values containing vertex numbers (in the VERTICES block) used to define the cell. Vertices must be listed in clockwise order. Cells that are connected must share vertices. | +| SWF | DISV2D | OPTIONS | LENGTH_UNITS | STRING | is the length units used for this model. Values can be ``FEET'', ``METERS'', or ``CENTIMETERS''. If not specified, the default is ``UNKNOWN''. | +| SWF | DISV2D | OPTIONS | NOGRB | KEYWORD | keyword to deactivate writing of the binary grid file. | +| SWF | DISV2D | OPTIONS | XORIGIN | DOUBLE PRECISION | x-position of the lower-left corner of the model grid. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | +| SWF | DISV2D | OPTIONS | YORIGIN | DOUBLE PRECISION | y-position of the lower-left corner of the model grid. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | +| SWF | DISV2D | OPTIONS | ANGROT | DOUBLE PRECISION | counter-clockwise rotation angle (in degrees) of the lower-left corner of the model grid. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | +| SWF | DISV2D | OPTIONS | EXPORT_ARRAY_ASCII | KEYWORD | keyword that specifies input grid arrays, which already support the layered keyword, should be written to layered ascii output files. | +| SWF | DISV2D | DIMENSIONS | NODES | INTEGER | is the number of cells per layer. This is a constant value for the grid and it applies to all layers. | +| SWF | DISV2D | DIMENSIONS | NVERT | INTEGER | is the total number of (x, y) vertex pairs used to characterize the horizontal configuration of the model grid. | +| SWF | DISV2D | GRIDDATA | BOTTOM | DOUBLE PRECISION (NODES) | is the bottom elevation for each cell. | +| SWF | DISV2D | GRIDDATA | IDOMAIN | INTEGER (NODES) | is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1 or greater, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. | +| SWF | DISV2D | VERTICES | IV | INTEGER | is the vertex number. Records in the VERTICES block must be listed in consecutive order from 1 to NVERT. | +| SWF | DISV2D | VERTICES | XV | DOUBLE PRECISION | is the x-coordinate for the vertex. | +| SWF | DISV2D | VERTICES | YV | DOUBLE PRECISION | is the y-coordinate for the vertex. | +| SWF | DISV2D | CELL2D | ICELL2D | INTEGER | is the CELL2D number. Records in the CELL2D block must be listed in consecutive order from the first to the last. | +| SWF | DISV2D | CELL2D | XC | DOUBLE PRECISION | is the x-coordinate for the cell center. | +| SWF | DISV2D | CELL2D | YC | DOUBLE PRECISION | is the y-coordinate for the cell center. | +| SWF | DISV2D | CELL2D | NCVERT | INTEGER | is the number of vertices required to define the cell. There may be a different number of vertices for each cell. | +| SWF | DISV2D | CELL2D | ICVERT | INTEGER (NCVERT) | is an array of integer values containing vertex numbers (in the VERTICES block) used to define the cell. Vertices must be listed in clockwise order. Cells that are connected must share vertices. | | SWF | DFW | OPTIONS | CENTRAL_IN_SPACE | KEYWORD | keyword to indicate conductance should be calculated using central-in-space weighting instead of the default upstream weighting approach. This option should be used with caution as it does not work well unless all of the stream reaches are saturated. With this option, there is no way for water to flow into a dry reach from connected reaches. | | SWF | DFW | OPTIONS | LENGTH_CONVERSION | DOUBLE PRECISION | real value that is used to convert user-specified Manning's roughness coefficients from meters to model length units. LENGTH\_CONVERSION should be set to 3.28081, 1.0, and 100.0 when using length units (LENGTH\_UNITS) of feet, meters, or centimeters in the simulation, respectively. LENGTH\_CONVERSION does not need to be specified if LENGTH\_UNITS are meters. | | SWF | DFW | OPTIONS | TIME_CONVERSION | DOUBLE PRECISION | real value that is used to convert user-specified Manning's roughness coefficients from seconds to model time units. TIME\_CONVERSION should be set to 1.0, 60.0, 3,600.0, 86,400.0, and 31,557,600.0 when using time units (TIME\_UNITS) of seconds, minutes, hours, days, or years in the simulation, respectively. TIME\_CONVERSION does not need to be specified if TIME\_UNITS are seconds. | | SWF | DFW | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control. | | SWF | DFW | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that calculated flows between cells will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. This option can produce extremely large list files because all cell-by-cell flows are printed. It should only be used with the DFW Package for models that have a small number of cells. | +| SWF | DFW | OPTIONS | SAVE_VELOCITY | KEYWORD | keyword to indicate that x, y, and z components of velocity will be calculated at cell centers and written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control. If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled). Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable. | | SWF | DFW | OPTIONS | OBS6 | KEYWORD | keyword to specify that record corresponds to an observations file. | | SWF | DFW | OPTIONS | FILEIN | KEYWORD | keyword to specify that an input filename is expected next. | | SWF | DFW | OPTIONS | OBS6_FILENAME | STRING | name of input file to define observations for the DFW package. See the ``Observation utility'' section for instructions for preparing observation input files. Tables \ref{table:gwf-obstypetable} and \ref{table:gwt-obstypetable} lists observation type(s) supported by the DFW package. | +| SWF | DFW | OPTIONS | DEV_SWR_CONDUCTANCE | KEYWORD | use the conductance formulation in the Surface Water Routing (SWR) Process for MODFLOW-2005. | | SWF | DFW | GRIDDATA | MANNINGSN | DOUBLE PRECISION (NODES) | mannings roughness coefficient | | SWF | DFW | GRIDDATA | IDCXS | INTEGER (NODES) | integer value indication the cross section identifier in the Cross Section Package that applies to the reach. If not provided then reach will be treated as hydraulically wide. | | SWF | CXS | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of stream reach information will be written to the listing file immediately after it is read. | diff --git a/doc/mf6io/mf6ivar/mf6ivar.py b/doc/mf6io/mf6ivar/mf6ivar.py index 9b7874c24b9..1d82ce971cb 100644 --- a/doc/mf6io/mf6ivar/mf6ivar.py +++ b/doc/mf6io/mf6ivar/mf6ivar.py @@ -762,7 +762,7 @@ def write_appendix(texdir, allblocks): "swf-nam", "swf-disv1d", "swf-dis2d", - "swf-disv", + "swf-disv2d", "swf-dfw", "swf-cxs", "swf-ic", diff --git a/doc/mf6io/mf6ivar/tex/appendixA.tex b/doc/mf6io/mf6ivar/tex/appendixA.tex index 4fe5ef129c3..1108af6ba01 100644 --- a/doc/mf6io/mf6ivar/tex/appendixA.tex +++ b/doc/mf6io/mf6ivar/tex/appendixA.tex @@ -332,11 +332,11 @@ SWF & DIS2D & DIMENSIONS & yes \\ SWF & DIS2D & GRIDDATA & no \\ \hline -SWF & DISV & OPTIONS & yes \\ -SWF & DISV & DIMENSIONS & yes \\ -SWF & DISV & GRIDDATA & no \\ -SWF & DISV & VERTICES & yes \\ -SWF & DISV & CELL2D & yes \\ +SWF & DISV2D & OPTIONS & yes \\ +SWF & DISV2D & DIMENSIONS & yes \\ +SWF & DISV2D & GRIDDATA & no \\ +SWF & DISV2D & VERTICES & yes \\ +SWF & DISV2D & CELL2D & yes \\ \hline SWF & DFW & OPTIONS & yes \\ SWF & DFW & GRIDDATA & no \\ diff --git a/doc/mf6io/mf6ivar/tex/swf-dfw-desc.tex b/doc/mf6io/mf6ivar/tex/swf-dfw-desc.tex index e2fe42e2e1b..24c3c05eb63 100644 --- a/doc/mf6io/mf6ivar/tex/swf-dfw-desc.tex +++ b/doc/mf6io/mf6ivar/tex/swf-dfw-desc.tex @@ -13,6 +13,8 @@ \item \texttt{PRINT\_FLOWS}---keyword to indicate that calculated flows between cells will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. This option can produce extremely large list files because all cell-by-cell flows are printed. It should only be used with the DFW Package for models that have a small number of cells. +\item \texttt{SAVE\_VELOCITY}---keyword to indicate that x, y, and z components of velocity will be calculated at cell centers and written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control. If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled). Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable. + \item \texttt{OBS6}---keyword to specify that record corresponds to an observations file. \item \texttt{FILEIN}---keyword to specify that an input filename is expected next. diff --git a/doc/mf6io/mf6ivar/tex/swf-dfw-options.dat b/doc/mf6io/mf6ivar/tex/swf-dfw-options.dat index fb2a386b7c5..126bc5c9f27 100644 --- a/doc/mf6io/mf6ivar/tex/swf-dfw-options.dat +++ b/doc/mf6io/mf6ivar/tex/swf-dfw-options.dat @@ -4,5 +4,6 @@ BEGIN OPTIONS [TIME_CONVERSION ] [SAVE_FLOWS] [PRINT_FLOWS] + [SAVE_VELOCITY] [OBS6 FILEIN ] END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/swf-disl-cell2d.dat b/doc/mf6io/mf6ivar/tex/swf-disl-cell2d.dat deleted file mode 100644 index f8ad1f65f50..00000000000 --- a/doc/mf6io/mf6ivar/tex/swf-disl-cell2d.dat +++ /dev/null @@ -1,5 +0,0 @@ -BEGIN CELL2D - - - ... -END CELL2D diff --git a/doc/mf6io/mf6ivar/tex/swf-disv-griddata.dat b/doc/mf6io/mf6ivar/tex/swf-disv-griddata.dat deleted file mode 100644 index e263cb1d7bb..00000000000 --- a/doc/mf6io/mf6ivar/tex/swf-disv-griddata.dat +++ /dev/null @@ -1,8 +0,0 @@ -BEGIN GRIDDATA - TOP - -- READARRAY - BOTM [LAYERED] - -- READARRAY - [IDOMAIN [LAYERED] - -- READARRAY] -END GRIDDATA diff --git a/doc/mf6io/mf6ivar/tex/swf-disv1d-desc.tex b/doc/mf6io/mf6ivar/tex/swf-disv1d-desc.tex index b89b7eeff1f..e20a304428a 100644 --- a/doc/mf6io/mf6ivar/tex/swf-disv1d-desc.tex +++ b/doc/mf6io/mf6ivar/tex/swf-disv1d-desc.tex @@ -13,6 +13,8 @@ \item \texttt{angrot}---counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +\item \texttt{EXPORT\_ARRAY\_ASCII}---keyword that specifies input grid arrays, which already support the layered keyword, should be written to layered ascii output files. + \end{description} \item \textbf{Block: DIMENSIONS} diff --git a/doc/mf6io/mf6ivar/tex/swf-disv1d-options.dat b/doc/mf6io/mf6ivar/tex/swf-disv1d-options.dat index 67e3ed895ae..965538a5ef9 100644 --- a/doc/mf6io/mf6ivar/tex/swf-disv1d-options.dat +++ b/doc/mf6io/mf6ivar/tex/swf-disv1d-options.dat @@ -4,4 +4,5 @@ BEGIN OPTIONS [XORIGIN ] [YORIGIN ] [ANGROT ] + [EXPORT_ARRAY_ASCII] END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/swf-disv-cell2d.dat b/doc/mf6io/mf6ivar/tex/swf-disv2d-cell2d.dat similarity index 100% rename from doc/mf6io/mf6ivar/tex/swf-disv-cell2d.dat rename to doc/mf6io/mf6ivar/tex/swf-disv2d-cell2d.dat diff --git a/doc/mf6io/mf6ivar/tex/swf-disv-desc.tex b/doc/mf6io/mf6ivar/tex/swf-disv2d-desc.tex similarity index 57% rename from doc/mf6io/mf6ivar/tex/swf-disv-desc.tex rename to doc/mf6io/mf6ivar/tex/swf-disv2d-desc.tex index f8680fc27b7..7ccc0793f5e 100644 --- a/doc/mf6io/mf6ivar/tex/swf-disv-desc.tex +++ b/doc/mf6io/mf6ivar/tex/swf-disv2d-desc.tex @@ -7,11 +7,11 @@ \item \texttt{NOGRB}---keyword to deactivate writing of the binary grid file. -\item \texttt{xorigin}---x-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +\item \texttt{xorigin}---x-position of the lower-left corner of the model grid. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. -\item \texttt{yorigin}---y-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +\item \texttt{yorigin}---y-position of the lower-left corner of the model grid. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. -\item \texttt{angrot}---counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. +\item \texttt{angrot}---counter-clockwise rotation angle (in degrees) of the lower-left corner of the model grid. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. \item \texttt{EXPORT\_ARRAY\_ASCII}---keyword that specifies input grid arrays, which already support the layered keyword, should be written to layered ascii output files. @@ -19,9 +19,7 @@ \item \textbf{Block: DIMENSIONS} \begin{description} -\item \texttt{nlay}---is the number of layers in the model grid. - -\item \texttt{ncpl}---is the number of cells per layer. This is a constant value for the grid and it applies to all layers. +\item \texttt{nodes}---is the number of cells per layer. This is a constant value for the grid and it applies to all layers. \item \texttt{nvert}---is the total number of (x, y) vertex pairs used to characterize the horizontal configuration of the model grid. @@ -29,11 +27,9 @@ \item \textbf{Block: GRIDDATA} \begin{description} -\item \texttt{top}---is the top elevation for each cell in the top model layer. - -\item \texttt{botm}---is the bottom elevation for each cell. +\item \texttt{bottom}---is the bottom elevation for each cell. -\item \texttt{idomain}---is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. +\item \texttt{idomain}---is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1 or greater, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. \end{description} \item \textbf{Block: VERTICES} diff --git a/doc/mf6io/mf6ivar/tex/swf-disv-dimensions.dat b/doc/mf6io/mf6ivar/tex/swf-disv2d-dimensions.dat similarity index 63% rename from doc/mf6io/mf6ivar/tex/swf-disv-dimensions.dat rename to doc/mf6io/mf6ivar/tex/swf-disv2d-dimensions.dat index b05791a77b3..b4954718a77 100644 --- a/doc/mf6io/mf6ivar/tex/swf-disv-dimensions.dat +++ b/doc/mf6io/mf6ivar/tex/swf-disv2d-dimensions.dat @@ -1,5 +1,4 @@ BEGIN DIMENSIONS - NLAY - NCPL + NODES NVERT END DIMENSIONS diff --git a/doc/mf6io/mf6ivar/tex/swf-disv2d-griddata.dat b/doc/mf6io/mf6ivar/tex/swf-disv2d-griddata.dat new file mode 100644 index 00000000000..da1f967582c --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/swf-disv2d-griddata.dat @@ -0,0 +1,6 @@ +BEGIN GRIDDATA + BOTTOM + -- READARRAY + [IDOMAIN + -- READARRAY] +END GRIDDATA diff --git a/doc/mf6io/mf6ivar/tex/swf-disv-options.dat b/doc/mf6io/mf6ivar/tex/swf-disv2d-options.dat similarity index 100% rename from doc/mf6io/mf6ivar/tex/swf-disv-options.dat rename to doc/mf6io/mf6ivar/tex/swf-disv2d-options.dat diff --git a/doc/mf6io/mf6ivar/tex/swf-disv-vertices.dat b/doc/mf6io/mf6ivar/tex/swf-disv2d-vertices.dat similarity index 100% rename from doc/mf6io/mf6ivar/tex/swf-disv-vertices.dat rename to doc/mf6io/mf6ivar/tex/swf-disv2d-vertices.dat diff --git a/doc/mf6io/swf/disv2d.tex b/doc/mf6io/swf/disv2d.tex new file mode 100644 index 00000000000..115782fc63c --- /dev/null +++ b/doc/mf6io/swf/disv2d.tex @@ -0,0 +1,19 @@ +Input to the Discretization by Vertices in Two Dimensions (DISV2D) Package is read from the file that has type ``DISV2D6'' in the Name File. Only one discretization package (either DIS2D, DISV1D, or DISV2D) can be specified for a SWF Model. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/swf-disv2d-options.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/swf-disv2d-dimensions.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/swf-disv2d-griddata.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/swf-disv2d-vertices.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/swf-disv2d-cell2d.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/swf-disv2d-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/swf-disv2d-example.dat} diff --git a/doc/mf6io/swf/namefile.tex b/doc/mf6io/swf/namefile.tex index 922e8a23b52..6bae5913123 100644 --- a/doc/mf6io/swf/namefile.tex +++ b/doc/mf6io/swf/namefile.tex @@ -24,6 +24,7 @@ \subsubsection{Explanation of Variables} \hline DISV1D6 & Discretization by Vertices in 1D Input File \\ DIS2D6 & Structured Grid Discretization Input File \\ +DISV2D6 & Discretization by Vertices in 2D Input File \\ DFW6 & Diffusive Wave Package \\ CXS6 & Cross Section Package \\ OC6 & Output Control Option \\ diff --git a/doc/mf6io/swf/swf.tex b/doc/mf6io/swf/swf.tex index 6a4be6f7dbc..c351be0e11b 100644 --- a/doc/mf6io/swf/swf.tex +++ b/doc/mf6io/swf/swf.tex @@ -22,6 +22,10 @@ \subsection{Discretization by Vertices in One Dimension (DISV1D) Package} \subsection{Structured Discretization in Two Dimensions (DIS2D) Package} \input{swf/dis2d} +\newpage +\subsection{Discretization by Vertices in Two Dimensions (DISV2D) Package} +\input{swf/disv2d} + \newpage \subsection{Diffusive Wave (DFW) Package} \input{swf/dfw} diff --git a/make/makefile b/make/makefile index cb20fed9cf0..6f75627e51a 100644 --- a/make/makefile +++ b/make/makefile @@ -107,7 +107,7 @@ $(OBJDIR)/swf-zdgidm.o \ $(OBJDIR)/swf-namidm.o \ $(OBJDIR)/swf-icidm.o \ $(OBJDIR)/swf-flwidm.o \ -$(OBJDIR)/swf-disvidm.o \ +$(OBJDIR)/swf-disv2didm.o \ $(OBJDIR)/swf-disv1didm.o \ $(OBJDIR)/swf-dis2didm.o \ $(OBJDIR)/swf-dfwidm.o \ @@ -384,6 +384,7 @@ $(OBJDIR)/swf-obs.o \ $(OBJDIR)/swf-flw.o \ $(OBJDIR)/swf-dfw.o \ $(OBJDIR)/swf-cdb.o \ +$(OBJDIR)/Disv2d.o \ $(OBJDIR)/Dis2d.o \ $(OBJDIR)/prt-prp.o \ $(OBJDIR)/prt-oc.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index c9d9046c877..9a1a862c927 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -181,7 +181,7 @@ - + @@ -207,7 +207,8 @@ - + + diff --git a/src/Idm/selector/IdmSwfDfnSelector.f90 b/src/Idm/selector/IdmSwfDfnSelector.f90 index 21b1c2b4ca7..95f32f646b8 100644 --- a/src/Idm/selector/IdmSwfDfnSelector.f90 +++ b/src/Idm/selector/IdmSwfDfnSelector.f90 @@ -8,7 +8,7 @@ module IdmSwfDfnSelectorModule use SwfNamInputModule use SwfDisv1DInputModule use SwfDis2DInputModule - use SwfDisvInputModule + use SwfDisv2DInputModule use SwfCxsInputModule use SwfDfwInputModule use SwfIcInputModule @@ -50,8 +50,8 @@ function swf_param_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, swf_disv1d_param_definitions) case ('DIS2D') call set_param_pointer(input_definition, swf_dis2d_param_definitions) - case ('DISV') - call set_param_pointer(input_definition, swf_disv_param_definitions) + case ('DISV2D') + call set_param_pointer(input_definition, swf_disv2d_param_definitions) case ('CXS') call set_param_pointer(input_definition, swf_cxs_param_definitions) case ('DFW') @@ -82,8 +82,8 @@ function swf_aggregate_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, swf_disv1d_aggregate_definitions) case ('DIS2D') call set_param_pointer(input_definition, swf_dis2d_aggregate_definitions) - case ('DISV') - call set_param_pointer(input_definition, swf_disv_aggregate_definitions) + case ('DISV2D') + call set_param_pointer(input_definition, swf_disv2d_aggregate_definitions) case ('CXS') call set_param_pointer(input_definition, swf_cxs_aggregate_definitions) case ('DFW') @@ -114,8 +114,8 @@ function swf_block_definitions(subcomponent) result(input_definition) call set_block_pointer(input_definition, swf_disv1d_block_definitions) case ('DIS2D') call set_block_pointer(input_definition, swf_dis2d_block_definitions) - case ('DISV') - call set_block_pointer(input_definition, swf_disv_block_definitions) + case ('DISV2D') + call set_block_pointer(input_definition, swf_disv2d_block_definitions) case ('CXS') call set_block_pointer(input_definition, swf_cxs_block_definitions) case ('DFW') @@ -145,8 +145,8 @@ function swf_idm_multi_package(subcomponent) result(multi_package) multi_package = swf_disv1d_multi_package case ('DIS2D') multi_package = swf_dis2d_multi_package - case ('DISV') - multi_package = swf_disv_multi_package + case ('DISV2D') + multi_package = swf_disv2d_multi_package case ('CXS') multi_package = swf_cxs_multi_package case ('DFW') @@ -180,7 +180,7 @@ function swf_idm_integrated(subcomponent) result(integrated) integrated = .true. case ('DIS2D') integrated = .true. - case ('DISV') + case ('DISV2D') integrated = .true. case ('CXS') integrated = .true. diff --git a/src/Idm/swf-disv1didm.f90 b/src/Idm/swf-disv1didm.f90 index 4c1b963cf10..3534fe20b3c 100644 --- a/src/Idm/swf-disv1didm.f90 +++ b/src/Idm/swf-disv1didm.f90 @@ -16,6 +16,7 @@ module SwfDisv1DInputModule logical :: xorigin = .false. logical :: yorigin = .false. logical :: angrot = .false. + logical :: export_ascii = .false. logical :: nodes = .false. logical :: nvert = .false. logical :: length = .false. @@ -118,6 +119,23 @@ module SwfDisv1DInputModule .false. & ! timeseries ) + type(InputParamDefinitionType), parameter :: & + swfdisv1d_export_ascii = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'DISV1D', & ! subcomponent + 'OPTIONS', & ! block + 'EXPORT_ARRAY_ASCII', & ! tag name + 'EXPORT_ASCII', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + type(InputParamDefinitionType), parameter :: & swfdisv1d_nodes = InputParamDefinitionType & ( & @@ -347,6 +365,7 @@ module SwfDisv1DInputModule swfdisv1d_xorigin, & swfdisv1d_yorigin, & swfdisv1d_angrot, & + swfdisv1d_export_ascii, & swfdisv1d_nodes, & swfdisv1d_nvert, & swfdisv1d_length, & diff --git a/src/Idm/swf-disvidm.f90 b/src/Idm/swf-disv2didm.f90 similarity index 70% rename from src/Idm/swf-disvidm.f90 rename to src/Idm/swf-disv2didm.f90 index 0896e9b410c..d119a1866a2 100644 --- a/src/Idm/swf-disvidm.f90 +++ b/src/Idm/swf-disv2didm.f90 @@ -1,27 +1,25 @@ ! ** Do Not Modify! MODFLOW 6 system generated file. ** -module SwfDisvInputModule +module SwfDisv2DInputModule use ConstantsModule, only: LENVARNAME use InputDefinitionModule, only: InputParamDefinitionType, & InputBlockDefinitionType private - public swf_disv_param_definitions - public swf_disv_aggregate_definitions - public swf_disv_block_definitions - public SwfDisvParamFoundType - public swf_disv_multi_package + public swf_disv2d_param_definitions + public swf_disv2d_aggregate_definitions + public swf_disv2d_block_definitions + public SwfDisv2dParamFoundType + public swf_disv2d_multi_package - type SwfDisvParamFoundType + type SwfDisv2dParamFoundType logical :: length_units = .false. logical :: nogrb = .false. logical :: xorigin = .false. logical :: yorigin = .false. logical :: angrot = .false. logical :: export_ascii = .false. - logical :: nlay = .false. - logical :: ncpl = .false. + logical :: nodes = .false. logical :: nvert = .false. - logical :: top = .false. - logical :: botm = .false. + logical :: bottom = .false. logical :: idomain = .false. logical :: iv = .false. logical :: xv = .false. @@ -31,15 +29,15 @@ module SwfDisvInputModule logical :: yc = .false. logical :: ncvert = .false. logical :: icvert = .false. - end type SwfDisvParamFoundType + end type SwfDisv2dParamFoundType - logical :: swf_disv_multi_package = .false. + logical :: swf_disv2d_multi_package = .false. type(InputParamDefinitionType), parameter :: & - swfdisv_length_units = InputParamDefinitionType & + swfdisv2d_length_units = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'OPTIONS', & ! block 'LENGTH_UNITS', & ! tag name 'LENGTH_UNITS', & ! fortran variable @@ -53,10 +51,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_nogrb = InputParamDefinitionType & + swfdisv2d_nogrb = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'OPTIONS', & ! block 'NOGRB', & ! tag name 'NOGRB', & ! fortran variable @@ -70,10 +68,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_xorigin = InputParamDefinitionType & + swfdisv2d_xorigin = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'OPTIONS', & ! block 'XORIGIN', & ! tag name 'XORIGIN', & ! fortran variable @@ -87,10 +85,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_yorigin = InputParamDefinitionType & + swfdisv2d_yorigin = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'OPTIONS', & ! block 'YORIGIN', & ! tag name 'YORIGIN', & ! fortran variable @@ -104,10 +102,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_angrot = InputParamDefinitionType & + swfdisv2d_angrot = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'OPTIONS', & ! block 'ANGROT', & ! tag name 'ANGROT', & ! fortran variable @@ -121,10 +119,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_export_ascii = InputParamDefinitionType & + swfdisv2d_export_ascii = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'OPTIONS', & ! block 'EXPORT_ARRAY_ASCII', & ! tag name 'EXPORT_ASCII', & ! fortran variable @@ -138,13 +136,13 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_nlay = InputParamDefinitionType & + swfdisv2d_nodes = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'DIMENSIONS', & ! block - 'NLAY', & ! tag name - 'NLAY', & ! fortran variable + 'NODES', & ! tag name + 'NODES', & ! fortran variable 'INTEGER', & ! type '', & ! shape .true., & ! required @@ -155,27 +153,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_ncpl = InputParamDefinitionType & + swfdisv2d_nvert = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent - 'DIMENSIONS', & ! block - 'NCPL', & ! tag name - 'NCPL', & ! fortran variable - 'INTEGER', & ! type - '', & ! shape - .true., & ! required - .false., & ! multi-record - .false., & ! preserve case - .false., & ! layered - .false. & ! timeseries - ) - - type(InputParamDefinitionType), parameter :: & - swfdisv_nvert = InputParamDefinitionType & - ( & - 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'DIMENSIONS', & ! block 'NVERT', & ! tag name 'NVERT', & ! fortran variable @@ -189,15 +170,15 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_top = InputParamDefinitionType & + swfdisv2d_bottom = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'GRIDDATA', & ! block - 'TOP', & ! tag name - 'TOP', & ! fortran variable + 'BOTTOM', & ! tag name + 'BOTTOM', & ! fortran variable 'DOUBLE1D', & ! type - 'NCPL', & ! shape + 'NODES', & ! shape .true., & ! required .false., & ! multi-record .false., & ! preserve case @@ -206,44 +187,27 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_botm = InputParamDefinitionType & + swfdisv2d_idomain = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent - 'GRIDDATA', & ! block - 'BOTM', & ! tag name - 'BOTM', & ! fortran variable - 'DOUBLE2D', & ! type - 'NCPL NLAY', & ! shape - .true., & ! required - .false., & ! multi-record - .false., & ! preserve case - .true., & ! layered - .false. & ! timeseries - ) - - type(InputParamDefinitionType), parameter :: & - swfdisv_idomain = InputParamDefinitionType & - ( & - 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'GRIDDATA', & ! block 'IDOMAIN', & ! tag name 'IDOMAIN', & ! fortran variable - 'INTEGER2D', & ! type - 'NCPL NLAY', & ! shape + 'INTEGER1D', & ! type + 'NODES', & ! shape .false., & ! required .false., & ! multi-record .false., & ! preserve case - .true., & ! layered + .false., & ! layered .false. & ! timeseries ) type(InputParamDefinitionType), parameter :: & - swfdisv_iv = InputParamDefinitionType & + swfdisv2d_iv = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'VERTICES', & ! block 'IV', & ! tag name 'IV', & ! fortran variable @@ -257,10 +221,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_xv = InputParamDefinitionType & + swfdisv2d_xv = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'VERTICES', & ! block 'XV', & ! tag name 'XV', & ! fortran variable @@ -274,10 +238,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_yv = InputParamDefinitionType & + swfdisv2d_yv = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'VERTICES', & ! block 'YV', & ! tag name 'YV', & ! fortran variable @@ -291,10 +255,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_icell2d = InputParamDefinitionType & + swfdisv2d_icell2d = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'CELL2D', & ! block 'ICELL2D', & ! tag name 'ICELL2D', & ! fortran variable @@ -308,10 +272,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_xc = InputParamDefinitionType & + swfdisv2d_xc = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'CELL2D', & ! block 'XC', & ! tag name 'XC', & ! fortran variable @@ -325,10 +289,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_yc = InputParamDefinitionType & + swfdisv2d_yc = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'CELL2D', & ! block 'YC', & ! tag name 'YC', & ! fortran variable @@ -342,10 +306,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_ncvert = InputParamDefinitionType & + swfdisv2d_ncvert = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'CELL2D', & ! block 'NCVERT', & ! tag name 'NCVERT', & ! fortran variable @@ -359,10 +323,10 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_icvert = InputParamDefinitionType & + swfdisv2d_icvert = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'CELL2D', & ! block 'ICVERT', & ! tag name 'ICVERT', & ! fortran variable @@ -376,35 +340,33 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swf_disv_param_definitions(*) = & + swf_disv2d_param_definitions(*) = & [ & - swfdisv_length_units, & - swfdisv_nogrb, & - swfdisv_xorigin, & - swfdisv_yorigin, & - swfdisv_angrot, & - swfdisv_export_ascii, & - swfdisv_nlay, & - swfdisv_ncpl, & - swfdisv_nvert, & - swfdisv_top, & - swfdisv_botm, & - swfdisv_idomain, & - swfdisv_iv, & - swfdisv_xv, & - swfdisv_yv, & - swfdisv_icell2d, & - swfdisv_xc, & - swfdisv_yc, & - swfdisv_ncvert, & - swfdisv_icvert & + swfdisv2d_length_units, & + swfdisv2d_nogrb, & + swfdisv2d_xorigin, & + swfdisv2d_yorigin, & + swfdisv2d_angrot, & + swfdisv2d_export_ascii, & + swfdisv2d_nodes, & + swfdisv2d_nvert, & + swfdisv2d_bottom, & + swfdisv2d_idomain, & + swfdisv2d_iv, & + swfdisv2d_xv, & + swfdisv2d_yv, & + swfdisv2d_icell2d, & + swfdisv2d_xc, & + swfdisv2d_yc, & + swfdisv2d_ncvert, & + swfdisv2d_icvert & ] type(InputParamDefinitionType), parameter :: & - swfdisv_vertices = InputParamDefinitionType & + swfdisv2d_vertices = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'VERTICES', & ! block 'VERTICES', & ! tag name 'VERTICES', & ! fortran variable @@ -418,15 +380,15 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swfdisv_cell2d = InputParamDefinitionType & + swfdisv2d_cell2d = InputParamDefinitionType & ( & 'SWF', & ! component - 'DISV', & ! subcomponent + 'DISV2D', & ! subcomponent 'CELL2D', & ! block 'CELL2D', & ! tag name 'CELL2D', & ! fortran variable 'RECARRAY ICELL2D XC YC NCVERT ICVERT', & ! type - 'NCPL', & ! shape + 'NODES', & ! shape .true., & ! required .false., & ! multi-record .false., & ! preserve case @@ -435,14 +397,14 @@ module SwfDisvInputModule ) type(InputParamDefinitionType), parameter :: & - swf_disv_aggregate_definitions(*) = & + swf_disv2d_aggregate_definitions(*) = & [ & - swfdisv_vertices, & - swfdisv_cell2d & + swfdisv2d_vertices, & + swfdisv2d_cell2d & ] type(InputBlockDefinitionType), parameter :: & - swf_disv_block_definitions(*) = & + swf_disv2d_block_definitions(*) = & [ & InputBlockDefinitionType( & 'OPTIONS', & ! blockname @@ -476,4 +438,4 @@ module SwfDisvInputModule ) & ] -end module SwfDisvInputModule +end module SwfDisv2DInputModule diff --git a/src/Model/Discretization/Disv2d.f90 b/src/Model/Discretization/Disv2d.f90 new file mode 100644 index 00000000000..63d9490b371 --- /dev/null +++ b/src/Model/Discretization/Disv2d.f90 @@ -0,0 +1,1387 @@ +module Disv2dModule + + use ArrayReadersModule, only: ReadArray + use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: LINELENGTH, LENMEMPATH, LENVARNAME, & + DZERO, DONE, DHALF + use BaseDisModule, only: DisBaseType + use GeomUtilModule, only: get_node, get_ijk, get_jk + use InputOutputModule, only: URWORD, ulasav, & + ulaprufw, ubdsv1, ubdsv06, getunit, openfile + use SimModule, only: count_errors, store_error, store_error_unit, & + store_error_filename + use SimVariablesModule, only: errmsg, idm_context + use DisvGeom, only: DisvGeomType, line_unit_vector + use MemoryManagerModule, only: mem_allocate, mem_deallocate, mem_setptr + use MemoryManagerExtModule, only: mem_set_value, memorylist_remove + use TdisModule, only: kstp, kper, pertim, totim, delt + + implicit none + private + public disv2d_cr, Disv2dType + + !> @brief Vertex grid discretization + type, extends(DisBaseType) :: Disv2dType + integer(I4B), pointer :: nvert => null() !< number of x,y vertices + real(DP), dimension(:, :), pointer, contiguous :: vertices => null() !< cell vertices stored as 2d array of x and y + real(DP), dimension(:, :), pointer, contiguous :: cellxy => null() !< cell center stored as 2d array of x and y + integer(I4B), dimension(:), pointer, contiguous :: iavert => null() !< cell vertex pointer ia array + integer(I4B), dimension(:), pointer, contiguous :: javert => null() !< cell vertex pointer ja array + real(DP), dimension(:), pointer, contiguous :: bottom => null() !< bottom elevations for each cell (nodes) + integer(I4B), dimension(:), pointer, contiguous :: idomain => null() !< idomain (nodes) + + contains + + procedure :: dis_df => disv2d_df + procedure :: dis_da => disv2d_da + procedure :: disv2d_load + procedure :: get_dis_type => get_dis_type + procedure, public :: record_array + procedure, public :: record_srcdst_list_header + ! -- helper functions + procedure :: get_nodenumber_idx1 + procedure :: nodeu_to_string + procedure :: nodeu_to_array + procedure :: nodeu_from_string + procedure :: nodeu_from_cellid + procedure :: connection_normal + procedure :: connection_vector + procedure :: get_polyverts + ! -- private + procedure :: source_options + procedure :: source_dimensions + procedure :: source_griddata + procedure :: log_options + procedure :: log_dimensions + procedure :: log_griddata + procedure :: source_vertices + procedure :: source_cell2d + procedure :: define_cellverts + procedure :: grid_finalize + procedure :: connect + procedure :: write_grb + procedure :: allocate_scalars + procedure :: allocate_arrays + procedure :: get_cell2d_area + + end type Disv2dType + + type DisvFoundType + logical :: length_units = .false. + logical :: nogrb = .false. + logical :: xorigin = .false. + logical :: yorigin = .false. + logical :: angrot = .false. + logical :: nodes = .false. + logical :: nvert = .false. + logical :: bottom = .false. + logical :: idomain = .false. + logical :: iv = .false. + logical :: xv = .false. + logical :: yv = .false. + logical :: icell2d = .false. + logical :: xc = .false. + logical :: yc = .false. + logical :: ncvert = .false. + logical :: icvert = .false. + end type DisvFoundType + +contains + + !> @brief Create a new discretization by vertices object + !< + subroutine disv2d_cr(dis, name_model, input_mempath, inunit, iout) + ! -- dummy + 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 + ! -- local + type(Disv2dType), pointer :: disnew + ! -- formats + character(len=*), parameter :: fmtheader = & + "(1X, /1X, 'DISV -- VERTEX GRID DISCRETIZATION PACKAGE,', & + &' VERSION 1 : 12/23/2015 - INPUT READ FROM MEMPATH: ', A, //)" + ! + allocate (disnew) + dis => disnew + call disnew%allocate_scalars(name_model, input_mempath) + dis%inunit = inunit + dis%iout = iout + ! + ! -- If disv enabled + if (inunit > 0) then + ! + ! -- Identify package + if (iout > 0) then + write (iout, fmtheader) dis%input_mempath + end if + ! + ! -- load disv + call disnew%disv2d_load() + end if + ! + end subroutine disv2d_cr + + !> @brief Transfer IDM data into this discretization object + !< + subroutine disv2d_load(this) + ! -- dummy + class(Disv2dType) :: this + ! + ! -- source input data + call this%source_options() + call this%source_dimensions() + call this%source_griddata() + call this%source_vertices() + call this%source_cell2d() + ! + end subroutine disv2d_load + + !> @brief Define the discretization + !< + subroutine disv2d_df(this) + ! -- dummy + class(Disv2dType) :: this + ! + call this%grid_finalize() + ! + end subroutine disv2d_df + + subroutine disv2d_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + use MemoryManagerExtModule, only: memorylist_remove + use SimVariablesModule, only: idm_context + ! -- dummy + class(Disv2dType) :: this + ! -- local + + ! -- Deallocate idm memory + call memorylist_remove(this%name_model, 'DISV2D', idm_context) + + ! -- scalars + call mem_deallocate(this%nvert) + + ! -- arrays + call mem_deallocate(this%nodeuser) + call mem_deallocate(this%nodereduced) + call mem_deallocate(this%bottom) + call mem_deallocate(this%idomain) + + ! -- cdl hack for arrays for vertices and cell2d blocks + call mem_deallocate(this%vertices) + call mem_deallocate(this%cellxy) + call mem_deallocate(this%iavert) + call mem_deallocate(this%javert) + ! + ! -- DisBaseType deallocate + call this%DisBaseType%dis_da() + ! + ! -- Return + return + end subroutine disv2d_da + + ! !> @brief Deallocate variables + ! !< + ! subroutine disv2d_da(this) + ! ! -- dummy + ! class(Disv2dType) :: this + ! ! + ! ! -- Deallocate idm memory + ! call memorylist_remove(this%name_model, 'DISV', idm_context) + ! call memorylist_remove(component=this%name_model, & + ! context=idm_context) + ! ! + ! ! -- DisBaseType deallocate + ! call this%DisBaseType%dis_da() + ! ! + ! ! -- Deallocate scalars + ! call mem_deallocate(this%nvert) + ! ! + ! ! -- Deallocate Arrays + ! call mem_deallocate(this%nodereduced) + ! call mem_deallocate(this%nodeuser) + ! call mem_deallocate(this%vertices) + ! call mem_deallocate(this%cellxy) + ! call mem_deallocate(this%iavert) + ! call mem_deallocate(this%javert) + ! call mem_deallocate(this%bottom) + ! call mem_deallocate(this%idomain) + ! ! + ! end subroutine disv2d_da + + !> @brief Copy options from IDM into package + !< + subroutine source_options(this) + ! -- dummy + class(Disv2dType) :: this + ! -- locals + character(len=LENVARNAME), dimension(3) :: lenunits = & + &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS'] + type(DisvFoundType) :: found + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%lenuni, 'LENGTH_UNITS', this%input_mempath, & + lenunits, found%length_units) + call mem_set_value(this%nogrb, 'NOGRB', this%input_mempath, found%nogrb) + call mem_set_value(this%xorigin, 'XORIGIN', this%input_mempath, found%xorigin) + call mem_set_value(this%yorigin, 'YORIGIN', this%input_mempath, found%yorigin) + call mem_set_value(this%angrot, 'ANGROT', this%input_mempath, found%angrot) + ! + ! -- log values to list file + if (this%iout > 0) then + call this%log_options(found) + end if + ! + end subroutine source_options + + !> @brief Write user options to list file + !< + subroutine log_options(this, found) + ! -- dummy + class(Disv2dType) :: this + type(DisvFoundType), intent(in) :: found + ! + write (this%iout, '(1x,a)') 'Setting Discretization Options' + ! + if (found%length_units) then + write (this%iout, '(4x,a,i0)') 'Model length unit [0=UND, 1=FEET, & + &2=METERS, 3=CENTIMETERS] set as ', this%lenuni + end if + ! + if (found%nogrb) then + write (this%iout, '(4x,a,i0)') 'Binary grid file [0=GRB, 1=NOGRB] & + &set as ', this%nogrb + end if + ! + if (found%xorigin) then + write (this%iout, '(4x,a,G0)') 'XORIGIN = ', this%xorigin + end if + ! + if (found%yorigin) then + write (this%iout, '(4x,a,G0)') 'YORIGIN = ', this%yorigin + end if + ! + if (found%angrot) then + write (this%iout, '(4x,a,G0)') 'ANGROT = ', this%angrot + end if + ! + write (this%iout, '(1x,a,/)') 'End Setting Discretization Options' + ! + end subroutine log_options + + !> @brief Copy dimensions from IDM into package + !< + subroutine source_dimensions(this) + ! -- dummy + class(Disv2dType) :: this + ! -- locals + integer(I4B) :: j + type(DisvFoundType) :: found + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%nodes, 'NODES', this%input_mempath, found%nodes) + call mem_set_value(this%nvert, 'NVERT', this%input_mempath, found%nvert) + ! + ! -- log simulation values + if (this%iout > 0) then + call this%log_dimensions(found) + end if + ! + ! -- verify dimensions were set + if (this%nodes < 1) then + call store_error( & + 'NODES was not specified or was specified incorrectly.') + call store_error_filename(this%input_fname) + end if + if (this%nvert < 1) then + call store_error( & + 'NVERT was not specified or was specified incorrectly.') + call store_error_filename(this%input_fname) + end if + ! + ! -- Calculate nodesuser + this%nodesuser = this%nodes + ! + ! -- Allocate non-reduced vectors for disv + call mem_allocate(this%idomain, this%nodes, 'IDOMAIN', & + this%memoryPath) + call mem_allocate(this%bottom, this%nodes, 'BOTTOM', & + this%memoryPath) + ! + ! -- Allocate vertices array + call mem_allocate(this%vertices, 2, this%nvert, 'VERTICES', this%memoryPath) + call mem_allocate(this%cellxy, 2, this%nodes, 'CELLXY', this%memoryPath) + ! + ! -- initialize all cells to be active (idomain = 1) + do j = 1, this%nodes + this%idomain(j) = 1 + end do + ! + end subroutine source_dimensions + + !> @brief Write dimensions to list file + !< + subroutine log_dimensions(this, found) + ! -- dummy + class(Disv2dType) :: this + type(DisvFoundType), intent(in) :: found + ! + write (this%iout, '(1x,a)') 'Setting Discretization Dimensions' + ! + if (found%nodes) then + write (this%iout, '(4x,a,i0)') 'NODES = ', this%nodes + end if + ! + if (found%nvert) then + write (this%iout, '(4x,a,i0)') 'NVERT = ', this%nvert + end if + ! + write (this%iout, '(1x,a,/)') 'End Setting Discretization Dimensions' + ! + end subroutine log_dimensions + + !> @brief Copy grid data from IDM into package + !< + subroutine source_griddata(this) + ! -- dummy + class(Disv2dType) :: this + ! -- locals + type(DisvFoundType) :: found + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%bottom, 'BOTTOM', this%input_mempath, found%bottom) + call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain) + ! + ! -- log simulation values + if (this%iout > 0) then + call this%log_griddata(found) + end if + ! + end subroutine source_griddata + + !> @brief Write griddata found to list file + !< + subroutine log_griddata(this, found) + ! -- dummy + class(Disv2dType) :: this + type(DisvFoundType), intent(in) :: found + ! + write (this%iout, '(1x,a)') 'Setting Discretization Griddata' + ! + if (found%bottom) then + write (this%iout, '(4x,a)') 'BOTTOM set from input file' + end if + ! + if (found%idomain) then + write (this%iout, '(4x,a)') 'IDOMAIN set from input file' + end if + ! + write (this%iout, '(1x,a,/)') 'End Setting Discretization Griddata' + ! + end subroutine log_griddata + + !> @brief Finalize grid (check properties, allocate arrays, compute connections) + !< + subroutine grid_finalize(this) + ! -- dummy + class(Disv2dType) :: this + ! -- locals + integer(I4B) :: node, noder, j, ncell_count + ! -- formats + character(len=*), parameter :: fmtnr = & + "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',& + &/1x, 'Number of user nodes: ',I0,& + &/1X, 'Number of nodes in solution: ', I0, //)" + ! + ! -- count active cells + ncell_count = 0 + do j = 1, this%nodes + if (this%idomain(j) > 0) ncell_count = ncell_count + 1 + end do + ! + ! -- Check to make sure nodes is a valid number + if (ncell_count == 0) then + call store_error('Model does not have any active nodes. & + &Ensure IDOMAIN array has some values greater & + &than zero.') + call store_error_filename(this%input_fname) + end if + + if (count_errors() > 0) then + call store_error_filename(this%input_fname) + end if + ! + ! -- Array size is now known, so allocate + call this%allocate_arrays() + ! + ! -- Fill the nodereduced array with the reduced nodenumber, or + ! a negative number to indicate it is a pass-through cell, or + ! a zero to indicate that the cell is excluded from the + ! solution. + if (this%nodes < this%nodesuser) then + node = 1 + noder = 1 + do j = 1, this%nodes + if (this%idomain(j) > 0) then + this%nodereduced(node) = noder + noder = noder + 1 + else + this%nodereduced(node) = 0 + end if + node = node + 1 + end do + end if + ! + ! -- allocate and fill nodeuser if a reduced grid + if (this%nodes < this%nodesuser) then + node = 1 + noder = 1 + do j = 1, this%nodes + if (this%idomain(j) > 0) then + this%nodeuser(noder) = node + noder = noder + 1 + end if + node = node + 1 + end do + end if + + ! Copy bottom into bot + do node = 1, this%nodesuser + this%bot(node) = this%bottom(node) + end do + + ! -- Move bottom into bot + ! and set x and y center coordinates + node = 0 + do j = 1, this%nodes + node = node + 1 + noder = node + if (this%nodes < this%nodesuser) noder = this%nodereduced(node) + if (noder <= 0) cycle + this%bot(noder) = this%bottom(j) + this%xc(noder) = this%cellxy(1, j) + this%yc(noder) = this%cellxy(2, j) + end do + ! + ! -- Build connections + call this%connect() + ! + end subroutine grid_finalize + + !> @brief Load grid vertices from IDM into package + !< + subroutine source_vertices(this) + ! -- dummy + class(Disv2dType) :: this + ! -- local + integer(I4B) :: i + real(DP), dimension(:), contiguous, pointer :: vert_x => null() + real(DP), dimension(:), contiguous, pointer :: vert_y => null() + ! + ! -- set pointers to memory manager input arrays + call mem_setptr(vert_x, 'XV', this%input_mempath) + call mem_setptr(vert_y, 'YV', this%input_mempath) + ! + ! -- set vertices 2d array + if (associated(vert_x) .and. associated(vert_y)) then + do i = 1, this%nvert + this%vertices(1, i) = vert_x(i) + this%vertices(2, i) = vert_y(i) + end do + else + call store_error('Required Vertex arrays not found.') + end if + ! + ! -- log + if (this%iout > 0) then + write (this%iout, '(1x,a)') 'Discretization Vertex data loaded' + end if + ! + end subroutine source_vertices + + !> @brief Build data structures to hold cell vertex info + !< + subroutine define_cellverts(this, icell2d, ncvert, icvert) + ! -- modules + use SparseModule, only: sparsematrix + ! -- dummy + class(Disv2dType) :: this + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icell2d + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: ncvert + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icvert + ! -- locals + type(sparsematrix) :: vert_spm + integer(I4B) :: i, j, ierr + integer(I4B) :: icv_idx, startvert, maxnnz = 5 + ! + ! -- initialize sparse matrix + call vert_spm%init(this%nodes, this%nvert, maxnnz) + ! + ! -- add sparse matrix connections from input memory paths + icv_idx = 1 + do i = 1, this%nodes + if (icell2d(i) /= i) call store_error('ICELL2D input sequence violation.') + do j = 1, ncvert(i) + call vert_spm%addconnection(i, icvert(icv_idx), 0) + if (j == 1) then + startvert = icvert(icv_idx) + elseif (j == ncvert(i) .and. (icvert(icv_idx) /= startvert)) then + call vert_spm%addconnection(i, startvert, 0) + end if + icv_idx = icv_idx + 1 + end do + end do + ! + ! -- allocate and fill iavert and javert + call mem_allocate(this%iavert, this%nodes + 1, 'IAVERT', this%memoryPath) + call mem_allocate(this%javert, vert_spm%nnz, 'JAVERT', this%memoryPath) + call vert_spm%filliaja(this%iavert, this%javert, ierr) + call vert_spm%destroy() + ! + end subroutine define_cellverts + + !> @brief Copy cell2d data from IDM into package + !< + subroutine source_cell2d(this) + ! -- dummy + class(Disv2dType) :: this + ! -- locals + integer(I4B), dimension(:), contiguous, pointer :: icell2d => null() + integer(I4B), dimension(:), contiguous, pointer :: ncvert => null() + integer(I4B), dimension(:), contiguous, pointer :: icvert => null() + real(DP), dimension(:), contiguous, pointer :: cell_x => null() + real(DP), dimension(:), contiguous, pointer :: cell_y => null() + integer(I4B) :: i + ! + ! -- set pointers to input path ncvert and icvert + call mem_setptr(icell2d, 'ICELL2D', this%input_mempath) + call mem_setptr(ncvert, 'NCVERT', this%input_mempath) + call mem_setptr(icvert, 'ICVERT', this%input_mempath) + ! + ! -- + if (associated(icell2d) .and. associated(ncvert) & + .and. associated(icvert)) then + call this%define_cellverts(icell2d, ncvert, icvert) + else + call store_error('Required cell vertex array(s) [ICELL2D, NCVERT, ICVERT] & + ¬ found.') + end if + ! + ! -- copy cell center idm sourced values to local arrays + call mem_setptr(cell_x, 'XC', this%input_mempath) + call mem_setptr(cell_y, 'YC', this%input_mempath) + ! + ! -- set cell centers + if (associated(cell_x) .and. associated(cell_y)) then + do i = 1, this%nodes + this%cellxy(1, i) = cell_x(i) + this%cellxy(2, i) = cell_y(i) + end do + else + call store_error('Required cell center arrays not found.') + end if + ! + ! -- log + if (this%iout > 0) then + write (this%iout, '(1x,a)') 'Discretization Cell2d data loaded' + end if + ! + end subroutine source_cell2d + + !> @brief Build grid connections + !< + subroutine connect(this) + ! -- dummy + class(Disv2dType) :: this + ! -- local + integer(I4B) :: j + integer(I4B) :: noder, nrsize + integer(I4B) :: narea_eq_zero + integer(I4B) :: narea_lt_zero + real(DP) :: area + ! + ! -- Initialize + narea_eq_zero = 0 + narea_lt_zero = 0 + ! + ! -- Assign the cell area + do j = 1, this%nodes + area = this%get_cell2d_area(j) + noder = this%get_nodenumber(j, 0) + if (noder > 0) this%area(noder) = area + if (area < DZERO) then + narea_lt_zero = narea_lt_zero + 1 + write (errmsg, '(a,i0,a)') & + &'Calculated CELL2D area less than zero for cell ', j, '.' + call store_error(errmsg) + end if + if (area == DZERO) then + narea_eq_zero = narea_eq_zero + 1 + write (errmsg, '(a,i0,a)') & + 'Calculated CELL2D area is zero for cell ', j, '.' + call store_error(errmsg) + end if + end do + ! + ! -- check for errors + if (count_errors() > 0) then + if (narea_lt_zero > 0) then + write (errmsg, '(i0,a)') narea_lt_zero, & + ' cell(s) have an area less than zero. Calculated cell & + &areas must be greater than zero. Negative areas often & + &mean vertices are not listed in clockwise order.' + call store_error(errmsg) + end if + if (narea_eq_zero > 0) then + write (errmsg, '(i0,a)') narea_eq_zero, & + ' cell(s) have an area equal to zero. Calculated cell & + &areas must be greater than zero. Calculated cell & + &areas equal to zero indicate that the cell is not defined & + &by a valid polygon.' + call store_error(errmsg) + end if + call store_error_filename(this%input_fname) + end if + ! + ! -- create and fill the connections object + nrsize = 0 + if (this%nodes < this%nodesuser) nrsize = this%nodes + allocate (this%con) + call this%con%disvconnections(this%name_model, this%nodes, & + this%nodes, 1, nrsize, & + this%nvert, this%vertices, this%iavert, & + this%javert, this%cellxy, & + this%bot, this%bot, & + this%nodereduced, this%nodeuser) + this%nja = this%con%nja + this%njas = this%con%njas + ! + end subroutine connect + + !> @brief Write a binary grid file + !< + subroutine write_grb(this, icelltype) + ! -- modules + use OpenSpecModule, only: access, form + ! -- dummy + class(Disv2dType) :: this + integer(I4B), dimension(:), intent(in) :: icelltype + ! -- local + integer(I4B) :: iunit, i, ntxt + integer(I4B), parameter :: lentxt = 100 + character(len=50) :: txthdr + character(len=lentxt) :: txt + character(len=LINELENGTH) :: fname + ! -- formats + character(len=*), parameter :: fmtgrdsave = & + "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', & + &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)" + ! + ! -- Initialize + ntxt = 18 + ! + ! -- Open the file + fname = trim(this%input_fname)//'.grb' + iunit = getunit() + write (this%iout, fmtgrdsave) iunit, trim(adjustl(fname)) + call openfile(iunit, this%iout, trim(adjustl(fname)), 'DATA(BINARY)', & + form, access, 'REPLACE') + ! + ! -- write header information + write (txthdr, '(a)') 'GRID DISV2D' + txthdr(50:50) = new_line('a') + write (iunit) txthdr + write (txthdr, '(a)') 'VERSION 1' + txthdr(50:50) = new_line('a') + write (iunit) txthdr + write (txthdr, '(a, i0)') 'NTXT ', ntxt + txthdr(50:50) = new_line('a') + write (iunit) txthdr + write (txthdr, '(a, i0)') 'LENTXT ', lentxt + txthdr(50:50) = new_line('a') + write (iunit) txthdr + ! + ! -- write variable definitions + write (txt, '(3a, i0)') 'NCELLS ', 'INTEGER ', 'NDIM 0 # ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'NODES ', 'INTEGER ', 'NDIM 0 # ', this%nodes + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'NVERT ', 'INTEGER ', 'NDIM 0 # ', this%nvert + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'NJAVERT ', 'INTEGER ', 'NDIM 0 # ', size(this%javert) + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'NJA ', 'INTEGER ', 'NDIM 0 # ', this%con%nja + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, 1pg25.15e3)') & + 'XORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%xorigin + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, 1pg25.15e3)') & + 'YORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%yorigin + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, 1pg25.15e3)') 'ANGROT ', 'DOUBLE ', 'NDIM 0 # ', this%angrot + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'BOTM ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'VERTICES ', 'DOUBLE ', 'NDIM 2 2 ', this%nvert + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'CELLX ', 'DOUBLE ', 'NDIM 1 ', this%nodes + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'CELLY ', 'DOUBLE ', 'NDIM 1 ', this%nodes + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'IAVERT ', 'INTEGER ', 'NDIM 1 ', this%nodes + 1 + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'JAVERT ', 'INTEGER ', 'NDIM 1 ', size(this%javert) + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'IA ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1 + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', size(this%con%jausr) + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + ! + ! -- write data + write (iunit) this%nodesuser ! ncells + write (iunit) this%nodes ! nodes + write (iunit) this%nvert ! nvert + write (iunit) size(this%javert) ! njavert + write (iunit) this%nja ! nja + write (iunit) this%xorigin ! xorigin + write (iunit) this%yorigin ! yorigin + write (iunit) this%angrot ! angrot + write (iunit) this%bottom ! botm + write (iunit) this%vertices ! vertices + write (iunit) (this%cellxy(1, i), i=1, this%nodes) ! cellx + write (iunit) (this%cellxy(2, i), i=1, this%nodes) ! celly + write (iunit) this%iavert ! iavert + write (iunit) this%javert ! javert + write (iunit) this%con%iausr ! iausr + write (iunit) this%con%jausr ! jausr + write (iunit) this%idomain ! idomain + write (iunit) icelltype ! icelltype + ! + ! -- Close the file + close (iunit) + ! + end subroutine write_grb + + !> @brief Convert a user nodenumber to a string (nodenumber) or (k,j) + !< + subroutine nodeu_to_string(this, nodeu, str) + ! -- dummy + class(Disv2dType) :: this + integer(I4B), intent(in) :: nodeu + character(len=*), intent(inout) :: str + ! -- local + integer(I4B) :: i, j, k + character(len=10) :: jstr + ! + call get_ijk(nodeu, 1, this%nodes, 1, i, j, k) + write (jstr, '(i10)') j + str = '('//trim(adjustl(jstr))//')' + ! + end subroutine nodeu_to_string + + !> @brief Convert a user nodenumber to an array (nodenumber) or (k,j) + !< + subroutine nodeu_to_array(this, nodeu, arr) + ! -- dummy + class(Disv2dType) :: this + integer(I4B), intent(in) :: nodeu + integer(I4B), dimension(:), intent(inout) :: arr + ! -- local + integer(I4B) :: isize + integer(I4B) :: i, j, k + ! + ! -- check the size of arr + isize = size(arr) + if (isize /= this%ndim) then + write (errmsg, '(a,i0,a,i0,a)') & + 'Program error: nodeu_to_array size of array (', isize, & + ') is not equal to the discretization dimension (', this%ndim, ').' + call store_error(errmsg, terminate=.TRUE.) + end if + ! + ! -- get k, i, j + call get_ijk(nodeu, 1, this%nodes, 1, i, j, k) + ! + ! -- fill array + arr(1) = j + ! + end subroutine nodeu_to_array + + !> @brief Get reduced node number from user node number + !< + function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) + ! -- return + integer(I4B) :: nodenumber + ! -- dummy + class(Disv2dType), intent(in) :: this + integer(I4B), intent(in) :: nodeu + integer(I4B), intent(in) :: icheck + ! -- local + ! + ! -- check the node number if requested + if (icheck /= 0) then + ! + ! -- If within valid range, convert to reduced nodenumber + if (nodeu < 1 .or. nodeu > this%nodesuser) then + nodenumber = 0 + write (errmsg, '(a,i0,a,i0,a)') & + 'Node number (', nodeu, ') is less than 1 or greater than nodes (', & + this%nodesuser, ').' + call store_error(errmsg) + else + nodenumber = nodeu + if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu) + end if + else + nodenumber = nodeu + if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu) + end if + ! + end function get_nodenumber_idx1 + + !> @brief Get normal vector components between the cell and a given neighbor + !! + !! The normal points outward from the shared face between noden and nodem. + !< + subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & + ipos) + ! -- dummy + class(Disv2dType) :: this + integer(I4B), intent(in) :: noden !< cell (reduced nn) + integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) + integer(I4B), intent(in) :: ihc !< horizontal connection flag + real(DP), intent(inout) :: xcomp + real(DP), intent(inout) :: ycomp + real(DP), intent(inout) :: zcomp + integer(I4B), intent(in) :: ipos + ! -- local + real(DP) :: angle, dmult + ! + ! -- Set vector components based on ihc + if (ihc == 0) then + xcomp = DZERO + ycomp = DZERO + if (nodem < noden) then + ! + ! -- nodem must be above noden, so upward connection + zcomp = DONE + else + ! + ! -- nodem must be below noden, so downward connection + zcomp = -DONE + end if + else + ! -- find from anglex, since anglex is symmetric, need to flip vector + ! for lower triangle (nodem < noden) + !ipos = this%con%getjaindex(noden, nodem) + angle = this%con%anglex(this%con%jas(ipos)) + dmult = DONE + if (nodem < noden) dmult = -DONE + xcomp = cos(angle) * dmult + ycomp = sin(angle) * dmult + zcomp = DZERO + end if + ! + end subroutine connection_normal + + !> @brief Get unit vector components between the cell and a given neighbor + !! + !! Saturation must be provided to compute cell center vertical coordinates. + !! Also return the straight-line connection length. + !< + subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & + xcomp, ycomp, zcomp, conlen) + ! -- dummy + class(Disv2dType) :: this + integer(I4B), intent(in) :: noden !< cell (reduced nn) + integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) + logical, intent(in) :: nozee !< do not use z in calculations + real(DP), intent(in) :: satn !< not used for disv1d + real(DP), intent(in) :: satm !< not used for disv1d + integer(I4B), intent(in) :: ihc !< horizontal connection flag + real(DP), intent(inout) :: xcomp !< x component of connection vector + real(DP), intent(inout) :: ycomp !< y component of connection vector + real(DP), intent(inout) :: zcomp !< z component of connection vector + real(DP), intent(inout) :: conlen !< calculated straight-line distance between cell centers + ! -- local + integer(I4B) :: nodeun, nodeum + real(DP) :: xn, xm, yn, ym, zn, zm + + ! horizontal connection, with possible z component due to cell offsets + ! and/or water table conditions + if (nozee) then + zn = DZERO + zm = DZERO + else + zn = this%bot(noden) + zm = this%bot(nodem) + end if + nodeun = this%get_nodeuser(noden) + nodeum = this%get_nodeuser(nodem) + xn = this%cellxy(1, nodeun) + yn = this%cellxy(2, nodeun) + xm = this%cellxy(1, nodeum) + ym = this%cellxy(2, nodeum) + call line_unit_vector(xn, yn, zn, xm, ym, zm, xcomp, ycomp, zcomp, & + conlen) + + end subroutine connection_vector + + !> @brief Get the discretization type + !< + subroutine get_dis_type(this, dis_type) + ! -- dummy + class(Disv2dType), intent(in) :: this + character(len=*), intent(out) :: dis_type + ! + dis_type = "DISV2D" + ! + end subroutine get_dis_type + + !> @brief Allocate and initialize scalars + !< + subroutine allocate_scalars(this, name_model, input_mempath) + ! -- dummy + class(Disv2dType) :: this + character(len=*), intent(in) :: name_model + character(len=*), intent(in) :: input_mempath + ! + ! -- Allocate parent scalars + call this%DisBaseType%allocate_scalars(name_model, input_mempath) + ! + ! -- Allocate + call mem_allocate(this%nvert, 'NVERT', this%memoryPath) + ! + ! -- Initialize + this%nvert = 0 + this%ndim = 1 + ! + end subroutine allocate_scalars + + !> @brief Allocate and initialize arrays + !< + subroutine allocate_arrays(this) + ! -- dummy + class(Disv2dType) :: this + ! + ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) + call this%DisBaseType%allocate_arrays() + ! + ! -- Allocate arrays for DisvType + if (this%nodes < this%nodesuser) then + call mem_allocate(this%nodeuser, this%nodes, 'NODEUSER', this%memoryPath) + call mem_allocate(this%nodereduced, this%nodesuser, 'NODEREDUCED', & + this%memoryPath) + else + call mem_allocate(this%nodeuser, 1, 'NODEUSER', this%memoryPath) + call mem_allocate(this%nodereduced, 1, 'NODEREDUCED', this%memoryPath) + end if + ! -- Initialize + this%mshape(1) = this%nodes + ! + end subroutine allocate_arrays + + !> @brief Get the signed area of the cell + !! + !! 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)] + !< + function get_cell2d_area(this, icell2d) result(area) + ! -- dummy + class(Disv2dType) :: this + integer(I4B), intent(in) :: icell2d + ! -- return + real(DP) :: area + ! -- local + integer(I4B) :: ivert + integer(I4B) :: nvert + integer(I4B) :: icount + integer(I4B) :: iv1 + real(DP) :: x + real(DP) :: y + real(DP) :: x1 + real(DP) :: y1 + ! + area = DZERO + nvert = this%iavert(icell2d + 1) - this%iavert(icell2d) + icount = 1 + iv1 = this%javert(this%iavert(icell2d)) + x1 = this%vertices(1, iv1) + y1 = this%vertices(2, iv1) + do ivert = this%iavert(icell2d), this%iavert(icell2d + 1) - 1 + x = this%vertices(1, this%javert(ivert)) + if (icount < nvert) then + y = this%vertices(2, this%javert(ivert + 1)) + else + y = this%vertices(2, this%javert(this%iavert(icell2d))) + end if + area = area + (x - x1) * (y - y1) + icount = icount + 1 + end do + ! + icount = 1 + do ivert = this%iavert(icell2d), this%iavert(icell2d + 1) - 1 + y = this%vertices(2, this%javert(ivert)) + if (icount < nvert) then + x = this%vertices(1, this%javert(ivert + 1)) + else + x = this%vertices(1, this%javert(this%iavert(icell2d))) + end if + area = area - (x - x1) * (y - y1) + icount = icount + 1 + end do + ! + area = -DONE * area * DHALF + ! + end function get_cell2d_area + + !> @brief Convert a string to a user nodenumber + !! + !! 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. + !< + function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & + flag_string, allow_zero) result(nodeu) + ! -- dummy + class(Disv2dType) :: 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, optional, intent(in) :: flag_string + logical, optional, intent(in) :: allow_zero + integer(I4B) :: nodeu + ! -- local + integer(I4B) :: j, nodes + integer(I4B) :: lloclocal, ndum, istat, n + real(DP) :: r + ! + if (present(flag_string)) then + if (flag_string) then + ! Check to see if first token in line can be read as an integer. + lloclocal = lloc + call urword(line, lloclocal, istart, istop, 1, ndum, r, iout, in) + read (line(istart:istop), *, iostat=istat) n + if (istat /= 0) then + ! First token in line is not an integer; return flag to this effect. + nodeu = -2 + return + end if + end if + end if + ! + nodes = this%mshape(1) + ! + call urword(line, lloc, istart, istop, 2, j, r, iout, in) + ! + if (j == 0) then + if (present(allow_zero)) then + if (allow_zero) then + nodeu = 0 + return + end if + end if + end if + ! + errmsg = '' + ! + if (j < 1 .or. j > nodes) then + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Cell number in list (', j, & + ') is outside of the grid.' + end if + ! + nodeu = get_node(1, 1, j, 1, 1, nodes) + ! + if (nodeu < 1 .or. nodeu > this%nodesuser) then + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), & + "Node number in list (", nodeu, ") is outside of the grid. "// & + "Cell number cannot be determined in line '"// & + trim(adjustl(line))//"'." + end if + ! + if (len_trim(adjustl(errmsg)) > 0) then + call store_error(errmsg) + call store_error_unit(in) + end if + ! + end function nodeu_from_string + + !> @brief Convert a cellid string to a 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. + !! + !! 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. + !< + function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & + allow_zero) result(nodeu) + ! -- return + integer(I4B) :: nodeu + ! -- dummy + class(Disv2dType) :: this + character(len=*), intent(inout) :: cellid + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + logical, optional, intent(in) :: flag_string + logical, optional, intent(in) :: allow_zero + ! -- local + integer(I4B) :: j, nodes + integer(I4B) :: lloclocal, ndum, istat, n + integer(I4B) :: istart, istop + real(DP) :: r + ! + if (present(flag_string)) then + if (flag_string) then + ! Check to see if first token in cellid can be read as an integer. + lloclocal = 1 + call urword(cellid, lloclocal, istart, istop, 1, ndum, r, iout, inunit) + read (cellid(istart:istop), *, iostat=istat) n + if (istat /= 0) then + ! First token in cellid is not an integer; return flag to this effect. + nodeu = -2 + return + end if + end if + end if + ! + nodes = this%mshape(1) + ! + lloclocal = 1 + call urword(cellid, lloclocal, istart, istop, 2, j, r, iout, inunit) + ! + if (j == 0) then + if (present(allow_zero)) then + if (allow_zero) then + nodeu = 0 + return + end if + end if + end if + ! + errmsg = '' + ! + if (j < 1 .or. j > nodes) then + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Cell2d number in list (', j, & + ') is outside of the grid.' + end if + ! + nodeu = get_node(1, 1, j, 1, 1, nodes) + ! + if (nodeu < 1 .or. nodeu > this%nodesuser) then + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), & + "Cell number cannot be determined for cellid ("// & + trim(adjustl(cellid))//") and results in a user "// & + "node number (", nodeu, ") that is outside of the grid." + end if + ! + if (len_trim(adjustl(errmsg)) > 0) then + call store_error(errmsg) + call store_error_unit(inunit) + end if + ! + end function nodeu_from_cellid + + !> @brief Get a 2D array of polygon vertices, listed in clockwise order + !! beginning with the lower left corner + !< + subroutine get_polyverts(this, ic, polyverts, closed) + ! -- dummy + class(Disv2dType), intent(inout) :: this + integer(I4B), intent(in) :: ic !< cell number (reduced) + real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices (column-major indexing) + logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex (default false) + ! -- local + integer(I4B) :: icu, icu2d, iavert, nverts, m, j + logical(LGP) :: lclosed + ! + ! count vertices + icu = this%get_nodeuser(ic) + icu2d = icu - ((icu - 1) / this%nodes) * this%nodes + nverts = this%iavert(icu2d + 1) - this%iavert(icu2d) - 1 + if (nverts .le. 0) nverts = nverts + size(this%javert) + ! + ! check closed option + if (.not. (present(closed))) then + lclosed = .false. + else + lclosed = closed + end if + ! + ! allocate vertices array + if (lclosed) then + allocate (polyverts(2, nverts + 1)) + else + allocate (polyverts(2, nverts)) + end if + ! + ! set vertices + iavert = this%iavert(icu2d) + do m = 1, nverts + j = this%javert(iavert - 1 + m) + polyverts(:, m) = (/this%vertices(1, j), this%vertices(2, j)/) + end do + ! + ! close if enabled + if (lclosed) & + polyverts(:, nverts + 1) = polyverts(:, 1) + ! + end subroutine + + !> @brief Record a double precision array + !! + !! The array is written to a formatted or unformatted external file depending + !! on the arguments. + !< + subroutine record_array(this, darray, iout, iprint, idataun, aname, & + cdatafmp, nvaluesp, nwidthp, editdesc, dinact) + ! -- dummy + class(Disv2dType), intent(inout) :: this + real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record + integer(I4B), intent(in) :: iout !< ascii output unit number + integer(I4B), intent(in) :: iprint !< whether to print the array + integer(I4B), intent(in) :: idataun !< binary output unit number, if negative don't write by layers, write entire array + character(len=*), intent(in) :: aname !< text descriptor + character(len=*), intent(in) :: cdatafmp !< write format + integer(I4B), intent(in) :: nvaluesp !< values per line + integer(I4B), intent(in) :: nwidthp !< number width + character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E) + real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain + ! -- local + integer(I4B) :: k, ifirst + integer(I4B) :: nlay + integer(I4B) :: nrow + integer(I4B) :: ncol + integer(I4B) :: nval + integer(I4B) :: nodeu, noder + integer(I4B) :: istart, istop + real(DP), dimension(:), pointer, contiguous :: dtemp + ! -- formats + character(len=*), parameter :: fmthsv = & + "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, & + &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)" + ! + ! -- set variables + nlay = 1 + nrow = 1 + ncol = this%mshape(1) + ! + ! -- If this is a reduced model, then copy the values from darray into + ! dtemp. + if (this%nodes < this%nodesuser) then + nval = this%nodes + dtemp => this%dbuff + do nodeu = 1, this%nodesuser + noder = this%get_nodenumber(nodeu, 0) + if (noder <= 0) then + dtemp(nodeu) = dinact + cycle + end if + dtemp(nodeu) = darray(noder) + end do + else + nval = this%nodes + dtemp => darray + end if + ! + ! -- Print to iout if iprint /= 0 + if (iprint /= 0) then + istart = 1 + do k = 1, nlay + istop = istart + nrow * ncol - 1 + call ulaprufw(ncol, nrow, kstp, kper, k, iout, dtemp(istart:istop), & + aname, cdatafmp, nvaluesp, nwidthp, editdesc) + istart = istop + 1 + end do + end if + ! + ! -- Save array to an external file. + if (idataun > 0) then + ! -- write to binary file by layer + ifirst = 1 + istart = 1 + do k = 1, nlay + istop = istart + nrow * ncol - 1 + if (ifirst == 1) write (iout, fmthsv) & + trim(adjustl(aname)), idataun, & + kstp, kper + ifirst = 0 + call ulasav(dtemp(istart:istop), aname, kstp, kper, & + pertim, totim, ncol, nrow, k, idataun) + istart = istop + 1 + end do + elseif (idataun < 0) then + ! + ! -- write entire array as one record + call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, & + iout, delt, pertim, totim) + end if + ! + end subroutine record_array + + !> @brief Record list header for imeth=6 + !< + subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & + dstmodel, dstpackage, naux, auxtxt, & + ibdchn, nlist, iout) + ! -- dummy + class(Disv2dType) :: 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 + ! -- local + integer(I4B) :: nlay, nrow, ncol + ! + nlay = 1 + nrow = 1 + ncol = this%mshape(1) + ! + ! -- Use ubdsv06 to write list header + call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, & + ibdchn, naux, auxtxt, ncol, nrow, nlay, & + nlist, iout, delt, pertim, totim) + ! + end subroutine record_srcdst_list_header + +end module Disv2dModule diff --git a/src/Model/SurfaceWaterFlow/swf-cdb.f90 b/src/Model/SurfaceWaterFlow/swf-cdb.f90 index fff076a427b..6a56c42f28c 100644 --- a/src/Model/SurfaceWaterFlow/swf-cdb.f90 +++ b/src/Model/SurfaceWaterFlow/swf-cdb.f90 @@ -23,7 +23,6 @@ module SwfCdbModule use InputOutputModule, only: GetUnit, openfile use MatrixBaseModule use BaseDisModule, only: DisBaseType - use Disv1dModule, only: Disv1dType use SwfCxsModule, only: SwfCxsType ! implicit none @@ -41,7 +40,6 @@ module SwfCdbModule real(DP), pointer :: gravconv => null() !< conversion factor gravity in m/s^2 to model units ! -- pointers other objects - type(Disv1dType), pointer :: disv1d type(SwfCxsType), pointer :: cxs contains @@ -111,13 +109,10 @@ subroutine cdb_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, & packobj%ncolbnd = 1 packobj%iscloc = 1 packobj%ictMemPath = create_mem_path(namemodel, 'DFW') - ! - ! -- store pointer to disv1d - select type (dis) - type is (Disv1dType) - cdbobj%disv1d => dis - end select - ! + + ! -- store pointer to dis + cdbobj%dis => dis + ! -- store pointer to cxs cdbobj%cxs => cxs ! diff --git a/src/Model/SurfaceWaterFlow/swf-dfw.f90 b/src/Model/SurfaceWaterFlow/swf-dfw.f90 index 8071e46a5eb..969e6b00068 100644 --- a/src/Model/SurfaceWaterFlow/swf-dfw.f90 +++ b/src/Model/SurfaceWaterFlow/swf-dfw.f90 @@ -184,6 +184,9 @@ subroutine dfw_df(this, dis) if (distype == "DIS2D") then this%is2d = 1 end if + if (distype == "DISV2D") then + this%is2d = 1 + end if ! -- check if dfw is enabled ! this will need to become if (.not. present(dfw_options)) then diff --git a/src/Model/SurfaceWaterFlow/swf.f90 b/src/Model/SurfaceWaterFlow/swf.f90 index b4825144e0e..d01a3d02bfc 100644 --- a/src/Model/SurfaceWaterFlow/swf.f90 +++ b/src/Model/SurfaceWaterFlow/swf.f90 @@ -116,7 +116,7 @@ module SwfModule !< integer(I4B), parameter :: SWF_NBASEPKG = 9 character(len=LENPACKAGETYPE), dimension(SWF_NBASEPKG) :: & - SWF_BASEPKG = ['DISV1D6', 'DIS2D6 ', 'DISV6 ', & + SWF_BASEPKG = ['DISV1D6', 'DIS2D6 ', 'DISV2D6', & 'DFW6 ', 'CXS6 ', 'OC6 ', & 'IC6 ', 'OBS6 ', 'STO6 '] @@ -1144,7 +1144,7 @@ subroutine create_packages(this) use SimVariablesModule, only: idm_context use Disv1dModule, only: disv1d_cr use Dis2dModule, only: dis2d_cr - use DisvModule, only: disv_cr + use Disv2dModule, only: disv2d_cr use SwfDfWModule, only: dfw_cr use SwfCxsModule, only: cxs_cr use SwfStoModule, only: sto_cr @@ -1198,9 +1198,9 @@ subroutine create_packages(this) case ('DIS2D6') indis = 1 call dis2d_cr(this%dis, this%name, mempath, indis, this%iout) - case ('DISV6') + case ('DISV2D6') indis = 1 - call disv_cr(this%dis, this%name, mempath, indis, this%iout) + call disv2d_cr(this%dis, this%name, mempath, indis, this%iout) case ('DFW6') this%indfw = 1 mempathdfw = mempath diff --git a/src/Utilities/Idm/SourceCommon.f90 b/src/Utilities/Idm/SourceCommon.f90 index ad8ecc29268..a61419dd137 100644 --- a/src/Utilities/Idm/SourceCommon.f90 +++ b/src/Utilities/Idm/SourceCommon.f90 @@ -333,6 +333,23 @@ subroutine set_model_shape(ftype, fname, model_mempath, dis_mempath, & else call store_error_filename(fname) end if + case ('DISV2D6') + ! + call get_isize('NODES', dis_mempath, dim1_size) + ! + if (dim1_size <= 0) then + write (errmsg, '(a)') & + 'Required input dimension "NODES" not found.' + call store_error(errmsg) + end if + ! + if (dim1_size >= 1) then + call mem_allocate(model_shape, 1, 'MODEL_SHAPE', model_mempath) + call mem_setptr(ndim1, 'NODES', dis_mempath) + model_shape = [ndim1] + else + call store_error_filename(fname) + end if case ('DISU6', 'DISV1D6') ! call get_isize('NODES', dis_mempath, dim1_size) diff --git a/src/Utilities/Idm/mf6blockfile/LoadMf6File.f90 b/src/Utilities/Idm/mf6blockfile/LoadMf6File.f90 index bd06ac99dc9..429ace441ba 100644 --- a/src/Utilities/Idm/mf6blockfile/LoadMf6File.f90 +++ b/src/Utilities/Idm/mf6blockfile/LoadMf6File.f90 @@ -279,7 +279,8 @@ recursive subroutine parse_block(this, iblk, recursive_call) ! ! -- disu vertices/cell2d blocks are contingent on NVERT dimension if (this%mf6_input%pkgtype == 'DISU6' .or. & - this%mf6_input%pkgtype == 'DISV1D6') then + this%mf6_input%pkgtype == 'DISV1D6' .or. & + this%mf6_input%pkgtype == 'DISV2D6') then if (this%mf6_input%block_dfns(iblk)%blockname == 'VERTICES' .or. & this%mf6_input%block_dfns(iblk)%blockname == 'CELL2D') then call get_from_memorylist('NVERT', this%mf6_input%mempath, mt, found, & diff --git a/src/meson.build b/src/meson.build index 92e13ea4e4c..005762bd07d 100644 --- a/src/meson.build +++ b/src/meson.build @@ -77,7 +77,7 @@ modflow_sources = files( 'Idm' / 'gwt-namidm.f90', 'Idm' / 'swf-disv1didm.f90', 'Idm' / 'swf-dis2didm.f90', - 'Idm' / 'swf-disvidm.f90', + 'Idm' / 'swf-disv2didm.f90', 'Idm' / 'swf-namidm.f90', 'Idm' / 'swf-cxsidm.f90', 'Idm' / 'swf-dfwidm.f90', @@ -118,6 +118,7 @@ modflow_sources = files( 'Model' / 'Discretization' / 'Disu.f90', 'Model' / 'Discretization' / 'Disv.f90', 'Model' / 'Discretization' / 'Disv1d.f90', + 'Model' / 'Discretization' / 'Disv2d.f90', 'Model' / 'Geometry' / 'BaseGeometry.f90', 'Model' / 'Geometry' / 'CircularGeometry.f90', 'Model' / 'Geometry' / 'RectangularGeometry.f90', diff --git a/utils/idmloader/dfns.txt b/utils/idmloader/dfns.txt index c1545782f1e..c2e09f3af85 100644 --- a/utils/idmloader/dfns.txt +++ b/utils/idmloader/dfns.txt @@ -45,7 +45,7 @@ gwe-dis.dfn swf-nam.dfn swf-disv1d.dfn swf-dis2d.dfn -swf-disv.dfn +swf-disv2d.dfn swf-cxs.dfn swf-dfw.dfn swf-ic.dfn