From ccfe0e9f8b4ab6fdecba41f1d02b6ee64d16932b Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Fri, 24 Jan 2025 19:38:49 -0500 Subject: [PATCH] feat(fmi): support binary grid file (#2153) Support an optional GWFGRID entry in FMI for the flow model's binary grid file. If provided, check the grids have the same size and idomain like we do in EXG as of #2149. This makes for the same error-checking in the exg (coupled) and fmi (post-processing) cases. Tentative idea is to make the grid file entry optional for some time (1-2 releases), then begin to require it. At which point the reader could be modified to load the full dis, and we could consolidate these haphazard checks into a robust grid comparison. The grb file reader in zonebudget is moved to a new GridFileReader module in the main src dir. The reader now loads headers and builds an index at init time, after which variables can be read as necessary. --- autotest/test_prt_fmi.py | 2 + autotest/test_prt_voronoi2.py | 2 + doc/ReleaseNotes/develop.tex | 3 +- doc/mf6io/gwe/fmi.tex | 1 + doc/mf6io/gwt/fmi.tex | 1 + doc/mf6io/mf6ivar/dfn/gwe-fmi.dfn | 2 +- doc/mf6io/mf6ivar/dfn/gwt-fmi.dfn | 2 +- doc/mf6io/mf6ivar/dfn/prt-fmi.dfn | 2 +- doc/mf6io/prt/fmi.tex | 1 + make/makefile | 17 +- msvs/mf6core.vfproj | 1 + .../ModelUtilities/FlowModelInterface.f90 | 200 +++++++--- src/Model/TransportModel/tsp-fmi.f90 | 1 + src/Utilities/GridFileReader.f90 | 350 ++++++++++++++++++ src/meson.build | 1 + utils/mf5to6/make/makefile | 2 +- utils/zonebudget/make/makefile | 5 +- utils/zonebudget/meson.build | 1 - utils/zonebudget/msvs/zonebudget.vfproj | 3 +- utils/zonebudget/pymake/extrafiles.txt | 2 + utils/zonebudget/src/grb.f90 | 168 --------- utils/zonebudget/src/zbud6.f90 | 9 +- 22 files changed, 545 insertions(+), 231 deletions(-) create mode 100644 src/Utilities/GridFileReader.f90 delete mode 100644 utils/zonebudget/src/grb.f90 diff --git a/autotest/test_prt_fmi.py b/autotest/test_prt_fmi.py index 1664cf13bf4..85f2969a789 100644 --- a/autotest/test_prt_fmi.py +++ b/autotest/test_prt_fmi.py @@ -135,9 +135,11 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): gwf_name = get_model_name(name, "gwf") gwf_budget_file = gwf_ws / f"{gwf_name}.bud" gwf_head_file = gwf_ws / f"{gwf_name}.hds" + grb_file = gwf_ws / f"{gwf_name}.dis.grb" flopy.mf6.ModflowPrtfmi( prt, packagedata=[ + ("GWFGRID", grb_file), ("GWFHEAD", gwf_head_file), ("GWFBUDGET", gwf_budget_file), ], diff --git a/autotest/test_prt_voronoi2.py b/autotest/test_prt_voronoi2.py index 59ff282a2bb..d5004e8cfad 100644 --- a/autotest/test_prt_voronoi2.py +++ b/autotest/test_prt_voronoi2.py @@ -175,9 +175,11 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets, cell_ids): ) gwf_budget_file = gwf_ws / f"{gwf_name}.bud" gwf_head_file = gwf_ws / f"{gwf_name}.hds" + grb_file = gwf_ws / f"{gwf_name}.disv.grb" flopy.mf6.ModflowPrtfmi( prt, packagedata=[ + ("GWFGRID", grb_file), ("GWFHEAD", gwf_head_file), ("GWFBUDGET", gwf_budget_file), ], diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index f45a91dc2fc..f3407df54de 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -20,7 +20,8 @@ \textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\ \underline{BASIC FUNCTIONALITY} \begin{itemize} - \item Previously, GWF Model exchanges with other model types (GWT, GWE, PRT) verified that the flow model and the coupled model had the same number of active nodes, but did not check that the models' IDOMAIN arrays selected precisely the same set. Exchanges will now check that model IDOMAIN arrays are identical. + \item GWT, GWE and PRT Models require a grid identical to the grid of the corresponding GWF Model. Previously, GWT, GWE and PRT exchanges performed a minimal check that the grids have the same number of active nodes. Exchanges will now also perform an additional check that grid IDOMAIN arrays are identical. + \item GWT, GWE and PRT FMI Packages can now read a GWF Model's binary grid file via a new GWFGRID entry. This allows FMI Packages to perform the same grid equivalence checks as exchanges, which guarantees identical error-checking behavior whether a GWT, GWE or PRT Model is coupled to a GWF Model or running as a post-processor. The GWFGRID file entry is optional but recommended. A future version may make the GWFGRID entry mandatory. \item A regression was recently introduced into the PRT model's generalized tracking method, in which a coordinate transformation was carried out prematurely. This could result in incorrect particle positions in and near quad-refined cells. This bug has been fixed. \item The PRT model previously allowed particles to be released at any time. Release times falling outside the bounds of the simulation's time discretization could produce undefined behavior. Any release times occurring before the simulation begins (i.e. negative times) will now be skipped with a warning message. If EXTEND\_TRACKING is not enabled, release times occurring after the end of the simulation will now be skipped with a warning message as well. If EXTEND\_TRACKING is enabled, release times after the end of the simulation are allowed. \end{itemize} diff --git a/doc/mf6io/gwe/fmi.tex b/doc/mf6io/gwe/fmi.tex index 898d93f1cc4..fbcaadf1818 100644 --- a/doc/mf6io/gwe/fmi.tex +++ b/doc/mf6io/gwe/fmi.tex @@ -31,6 +31,7 @@ \item If the binary budget and head files have more than one time step for a single stress period, then the budget and head information must be contained within the binary file for every time step in the simulation stress period. \item The binary budget and head files must correspond in terms of information stored for each time step and stress period. \item If the binary budget and head files have information provided for only the first time step of a given stress period, this information will be used for all time steps in that stress period in the GWE simulation. If the final stress period (which may be the only stress period) in the binary budget and head files has information provided for only one time step, this information will be used for any subsequent time steps and stress periods in the GWE simulation. This makes it possible to provide flows, for example, from a steady-state GWF stress period and have those flows used for all GWE time steps in that stress period, for all remaining time steps in the GWE simulation, or for all time steps throughout the entire GWE simulation. With this option, it is possible to have smaller time steps in the GWE simulation than the time steps used in the GWF simulation. Note that this cannot be done when the GWF and GWE models are run in the same simulation, because in that case, both models are solved for each time step in the stress period, as listed in the TDIS Package. This option for reading flows from a previous GWF simulation may offer an efficient alternative to running both models in the same simulation, but it comes at the cost of having potentially very large budget files. +\item The binary grid file is optional but recommended, as it allows \mf to verify that the GWF and GWE model grids are identical. \end{itemize} \end{itemize} diff --git a/doc/mf6io/gwt/fmi.tex b/doc/mf6io/gwt/fmi.tex index f951253b44d..7108f1ea301 100644 --- a/doc/mf6io/gwt/fmi.tex +++ b/doc/mf6io/gwt/fmi.tex @@ -31,6 +31,7 @@ \item If the binary budget and head files have more than one time step for a single stress period, then the budget and head information must be contained within the binary file for every time step in the simulation stress period. \item The binary budget and head files must correspond in terms of information stored for each time step and stress period. \item If the binary budget and head files have information provided for only the first time step of a given stress period, this information will be used for all time steps in that stress period in the GWT simulation. If the final stress period (which may be the only stress period) in the binary budget and head files has information provided for only one time step, this information will be used for any subsequent time steps and stress periods in the GWT simulation. This makes it possible to provide flows, for example, from a steady-state GWF stress period and have those flows used for all GWT time steps in that stress period, for all remaining time steps in the GWT simulation, or for all time steps throughout the entire GWT simulation. With this option, it is possible to have smaller time steps in the GWT simulation than the time steps used in the GWF simulation. Note that this cannot be done when the GWF and GWT models are run in the same simulation, because in that case, both models are solved for each time step in the stress period, as listed in the TDIS Package. This option for reading flows from a previous GWF simulation may offer an efficient alternative to running both models in the same simulation, but it comes at the cost of having potentially very large budget files. +\item The binary grid file is optional but recommended, as it allows \mf to verify that the GWF and GWT model grids are identical. \end{itemize} \end{itemize} diff --git a/doc/mf6io/mf6ivar/dfn/gwe-fmi.dfn b/doc/mf6io/mf6ivar/dfn/gwe-fmi.dfn index fb71131e605..0a4dfad871a 100644 --- a/doc/mf6io/mf6ivar/dfn/gwe-fmi.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwe-fmi.dfn @@ -33,7 +33,7 @@ type string tagged false reader urword longname flow type -description is the word GWFBUDGET, GWFHEAD, GWFMOVER or the name of an advanced GWF stress package. If GWFBUDGET is specified, then the corresponding file must be a budget file from a previous GWF Model run. If an advanced GWF stress package name appears then the corresponding file must be the budget file saved by a LAK, SFR, MAW or UZF Package. +description is the word GWFBUDGET, GWFHEAD, GWFGRID, GWFMOVER or the name of an advanced GWF stress package from a previous model run. If GWFBUDGET is specified, then the corresponding file must be a budget file. If GWFHEAD is specified, the file must be a head file. If GWFGRID is specified, the file must be a binary grid file. If GWFMOVER is specified, the file must be a mover file. If an advanced GWF stress package name appears then the corresponding file must be the budget file saved by a LAK, SFR, MAW or UZF Package. block packagedata name filein diff --git a/doc/mf6io/mf6ivar/dfn/gwt-fmi.dfn b/doc/mf6io/mf6ivar/dfn/gwt-fmi.dfn index b1293717932..5fe314ff5d7 100644 --- a/doc/mf6io/mf6ivar/dfn/gwt-fmi.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwt-fmi.dfn @@ -33,7 +33,7 @@ type string tagged false reader urword longname flow type -description is the word GWFBUDGET, GWFHEAD, GWFMOVER or the name of an advanced GWF stress package. If GWFBUDGET is specified, then the corresponding file must be a budget file from a previous GWF Model run. If an advanced GWF stress package name appears then the corresponding file must be the budget file saved by a LAK, SFR, MAW or UZF Package. +description is the word GWFBUDGET, GWFHEAD, GWFMOVER or the name of an advanced GWF stress package. If GWFBUDGET is specified, then the corresponding file must be a budget file. If GWFHEAD is specified, the file must be a head file. If GWFGRID is specified, the file must be a binary grid file. If an advanced GWF stress package name appears then the corresponding file must be the budget file saved by a LAK, SFR, MAW or UZF Package. block packagedata name filein diff --git a/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn b/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn index 3deb36646fd..37895ad3e70 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn +++ b/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn @@ -25,7 +25,7 @@ type string tagged false reader urword longname flow type -description is the word GWFBUDGET or GWFHEAD. If GWFBUDGET is specified, then the corresponding file must be a budget file from a previous GWF Model run. +description is the word GWFBUDGET, GWFHEAD, or GWFGRID. If GWFBUDGET is specified, then the corresponding file must be a budget file. If GWFHEAD is specified, the file must be a head file. If GWFGRID is specified, the file must be a binary grid file. block packagedata name filein diff --git a/doc/mf6io/prt/fmi.tex b/doc/mf6io/prt/fmi.tex index 76ec7ef112d..236a1d51f5b 100644 --- a/doc/mf6io/prt/fmi.tex +++ b/doc/mf6io/prt/fmi.tex @@ -28,6 +28,7 @@ \item If the binary budget and head files have more than one time step for a single stress period, then the budget and head information must be contained within the binary file for every time step in the simulation stress period. \item The binary budget and head files must correspond in terms of information stored for each time step and stress period. \item If the binary budget and head files have information provided for only the first time step of a given stress period, this information will be used for all time steps in that stress period in the PRT simulation. If the final (or only) stress period in the binary budget and head files contains data for only one time step, this information will be used for any subsequent time steps and stress periods in the PRT simulation. This makes it possible to provide flows, for example, from a steady-state GWF stress period and have those flows used for all PRT time steps in that stress period, for all remaining time steps in the PRT simulation, or for all time steps throughout the entire PRT simulation. With this option, it is possible to have smaller time steps in the PRT simulation than the time steps used in the GWF simulation. Note that this cannot be done when the GWF and PRT models are run in the same simulation, because in that case, both models are solved over the same sequence of time steps and stress periods, as listed in the TDIS Package. The option to read flows from a previous GWF simulation via Flow Model Interface may offer an efficient alternative to running both models in the same simulation, but comes at the cost of having potentially very large budget files. +\item The binary grid file is optional but recommended, as it allows \mf to verify that the GWF and PRT model grids are identical. \end{itemize} \end{itemize} diff --git a/make/makefile b/make/makefile index c2502cc85d0..89b829517ba 100644 --- a/make/makefile +++ b/make/makefile @@ -1,4 +1,4 @@ -# makefile created by pymake (version 1.2.11.dev0) for the 'mf6' executable. +# makefile created by pymake (version 1.4.0.dev0) for the 'mf6' executable. include ./makedefaults @@ -249,6 +249,7 @@ $(OBJDIR)/BudgetFileReader.o \ $(OBJDIR)/TimeArraySeriesLink.o \ $(OBJDIR)/ObsUtility.o \ $(OBJDIR)/ObsContainer.o \ +$(OBJDIR)/Disv1dGeom.o \ $(OBJDIR)/BudgetTerm.o \ $(OBJDIR)/Budget.o \ $(OBJDIR)/TimeArraySeriesManager.o \ @@ -258,6 +259,13 @@ $(OBJDIR)/NumericalPackage.o \ $(OBJDIR)/Particle.o \ $(OBJDIR)/PackageBudget.o \ $(OBJDIR)/HeadFileReader.o \ +$(OBJDIR)/GridFileReader.o \ +$(OBJDIR)/Disv.o \ +$(OBJDIR)/Disv2d.o \ +$(OBJDIR)/Disv1d.o \ +$(OBJDIR)/Disu.o \ +$(OBJDIR)/Dis.o \ +$(OBJDIR)/Dis2d.o \ $(OBJDIR)/BudgetObject.o \ $(OBJDIR)/BoundaryPackage.o \ $(OBJDIR)/CellDefn.o \ @@ -332,7 +340,6 @@ $(OBJDIR)/MethodCellTernary.o \ $(OBJDIR)/MethodCellPollockQuad.o \ $(OBJDIR)/MethodCellPollock.o \ $(OBJDIR)/MethodCellPassToBot.o \ -$(OBJDIR)/Disv1dGeom.o \ $(OBJDIR)/VectorInterpolation.o \ $(OBJDIR)/swf-cxs.o \ $(OBJDIR)/CellWithNbrs.o \ @@ -342,9 +349,6 @@ $(OBJDIR)/tsp-oc.o \ $(OBJDIR)/tsp-obs.o \ $(OBJDIR)/tsp-mvt.o \ $(OBJDIR)/tsp-adv.o \ -$(OBJDIR)/Disv.o \ -$(OBJDIR)/Disu.o \ -$(OBJDIR)/Dis.o \ $(OBJDIR)/gwf-uzf.o \ $(OBJDIR)/tsp-apt.o \ $(OBJDIR)/gwt-mst.o \ @@ -371,7 +375,6 @@ $(OBJDIR)/Double2dReader.o \ $(OBJDIR)/Double1dReader.o \ $(OBJDIR)/MethodCellPool.o \ $(OBJDIR)/CellUtil.o \ -$(OBJDIR)/Disv1d.o \ $(OBJDIR)/swf-dfw.o \ $(OBJDIR)/swf-ic.o \ $(OBJDIR)/VirtualExchange.o \ @@ -426,8 +429,6 @@ $(OBJDIR)/swf-oc.o \ $(OBJDIR)/swf-obs.o \ $(OBJDIR)/swf-flw.o \ $(OBJDIR)/swf-cdb.o \ -$(OBJDIR)/Disv2d.o \ -$(OBJDIR)/Dis2d.o \ $(OBJDIR)/GridConnection.o \ $(OBJDIR)/DistributedVariable.o \ $(OBJDIR)/gwt.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index c97b8e6638a..c003b23875c 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -628,6 +628,7 @@ + diff --git a/src/Model/ModelUtilities/FlowModelInterface.f90 b/src/Model/ModelUtilities/FlowModelInterface.f90 index 60a4729e257..00636e8bc8c 100644 --- a/src/Model/ModelUtilities/FlowModelInterface.f90 +++ b/src/Model/ModelUtilities/FlowModelInterface.f90 @@ -10,6 +10,7 @@ module FlowModelInterfaceModule use ListModule, only: ListType use BudgetFileReaderModule, only: BudgetFileReaderType use HeadFileReaderModule, only: HeadFileReaderType + use GridFileReaderModule, only: GridFileReaderType use PackageBudgetModule, only: PackageBudgetType use BudgetObjectModule, only: BudgetObjectType, budgetobject_cr_bfr @@ -23,7 +24,7 @@ module FlowModelInterfaceModule logical, pointer :: flows_from_file => null() !< if .false., then flows come from GWF through GWF-Model exg type(ListType), pointer :: gwfbndlist => null() !< list of gwf stress packages integer(I4B), pointer :: iflowsupdated => null() !< flows were updated for this time step - integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to Model ibound + integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to this model ibound real(DP), dimension(:), pointer, contiguous :: gwfflowja => null() !< pointer to the GWF flowja array real(DP), dimension(:, :), pointer, contiguous :: gwfspdis => null() !< pointer to npf specific discharge array real(DP), dimension(:), pointer, contiguous :: gwfhead => null() !< pointer to the GWF head array @@ -37,12 +38,14 @@ module FlowModelInterfaceModule integer(I4B), pointer :: iubud => null() !< unit number GWF budget file integer(I4B), pointer :: iuhds => null() !< unit number GWF head file integer(I4B), pointer :: iumvr => null() !< unit number GWF mover budget file + integer(I4B), pointer :: iugrb => null() !< unit number binary grid file integer(I4B), pointer :: nflowpack => null() !< number of GWF flow packages integer(I4B), dimension(:), pointer, contiguous :: igwfmvrterm => null() !< flag to indicate that gwf package is a mover term type(BudgetFileReaderType) :: bfr !< budget file reader type(HeadFileReaderType) :: hfr !< head file reader + type(GridFileReaderType) :: gfr !< grid file reader type(PackageBudgetType), dimension(:), allocatable :: gwfpackages !< used to get flows between a package and gwf - type(BudgetObjectType), pointer :: mvrbudobj => null() !< pointer to the mover budget budget object + type(BudgetObjectType), pointer :: mvrbudobj => null() !< pointer to the mover budget object character(len=16), dimension(:), allocatable :: flowpacknamearray !< array of boundary package names (e.g. LAK-1, SFR-3, etc.) character(len=LENVARNAME) :: depvartype = '' @@ -181,6 +184,7 @@ subroutine fmi_da(this) call mem_deallocate(this%iubud) call mem_deallocate(this%iuhds) call mem_deallocate(this%iumvr) + call mem_deallocate(this%iugrb) call mem_deallocate(this%nflowpack) call mem_deallocate(this%idryinactive) ! @@ -208,6 +212,7 @@ subroutine allocate_scalars(this) call mem_allocate(this%iubud, 'IUBUD', this%memoryPath) call mem_allocate(this%iuhds, 'IUHDS', this%memoryPath) call mem_allocate(this%iumvr, 'IUMVR', this%memoryPath) + call mem_allocate(this%iugrb, 'IUGRB', this%memoryPath) call mem_allocate(this%nflowpack, 'NFLOWPACK', this%memoryPath) call mem_allocate(this%idryinactive, "IDRYINACTIVE", this%memoryPath) ! @@ -220,6 +225,7 @@ subroutine allocate_scalars(this) this%iubud = 0 this%iuhds = 0 this%iumvr = 0 + this%iugrb = 0 this%nflowpack = 0 this%idryinactive = 1 end subroutine allocate_scalars @@ -331,6 +337,12 @@ subroutine read_packagedata(this) use ConstantsModule, only: LINELENGTH, DEM6, LENPACKAGENAME use InputOutputModule, only: getunit, openfile, urdaux use SimModule, only: store_error, store_error_unit + use DisModule, only: DisType + use DisvModule, only: DisvType + use DisuModule, only: DisuType + use Dis2dModule, only: Dis2dType + use Disv2dModule, only: Disv2dType + use Disv1dModule, only: Disv1dType ! -- dummy class(FlowModelInterfaceType) :: this ! -- local @@ -338,20 +350,34 @@ subroutine read_packagedata(this) integer(I4B) :: ierr integer(I4B) :: inunit integer(I4B) :: iapt + integer(I4B) :: nodes logical :: isfound, endOfBlock logical :: blockrequired logical :: exist + integer(I4B), allocatable :: idomain1d(:), idomain2d(:, :), idomain3d(:, :, :) + ! -- formats + character(len=*), parameter :: fmtdiserr = & + "('Models do not have the same discretization& + & ',a,'.& + & GWF model has ', i0, ' user nodes and ', i0, ' reduced nodes.& + & This model has ', i0, ' user nodes and ', i0, ' reduced nodes.& + & Ensure discretization packages, including IDOMAIN, are identical.')" + character(len=*), parameter :: fmtidomerr = & + "('Models do not have the same discretization& + & ',a,'.& + & Models have different IDOMAIN arrays.& + & Ensure discretization packages, including IDOMAIN, are identical.')" ! ! -- initialize iapt = 0 blockrequired = .true. ! - ! -- get options block + ! -- get packagedata block call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, & blockRequired=blockRequired, & supportOpenClose=.true.) ! - ! -- parse options block if detected + ! -- parse packagedata block if detected if (isfound) then write (this%iout, '(1x,a)') 'PROCESSING FMI PACKAGEDATA' do @@ -410,6 +436,126 @@ subroutine read_packagedata(this) call budgetobject_cr_bfr(this%mvrbudobj, 'MVT', this%iumvr, & this%iout) call this%mvrbudobj%fill_from_bfr(this%dis, this%iout) + case ('GWFGRID') + call this%parser%GetStringCaps(keyword) + if (keyword /= 'FILEIN') then + call store_error('GWFGRID KEYWORD MUST BE FOLLOWED BY '// & + '"FILEIN" then by filename.') + call this%parser%StoreErrorUnit() + end if + call this%parser%GetString(fname) + inunit = getunit() + call openfile(inunit, this%iout, fname, 'DATA(BINARY)', & + FORM, ACCESS, 'UNKNOWN') + this%iugrb = inunit + call this%gfr%initialize(this%iugrb) + + ! check grid equivalence + select case (this%gfr%grid_type) + case ('DIS') + select type (dis => this%dis) + type is (DisType) + nodes = this%gfr%read_int("NCELLS") + if (nodes /= this%dis%nodes) then + write (errmsg, fmtdiserr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + idomain1d = this%gfr%read_int_1d("IDOMAIN") + idomain3d = reshape(idomain1d, [ & + this%gfr%read_int("NCOL"), & + this%gfr%read_int("NROW"), & + this%gfr%read_int("NLAY") & + ]) + if (.not. all(dis%idomain == idomain3d)) then + write (errmsg, fmtidomerr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + end select + case ('DISV') + select type (dis => this%dis) + type is (DisvType) + nodes = this%gfr%read_int("NCELLS") + if (nodes /= this%dis%nodes) then + write (errmsg, fmtdiserr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + idomain1d = this%gfr%read_int_1d("IDOMAIN") + idomain2d = reshape(idomain1d, [ & + this%gfr%read_int("NCPL"), & + this%gfr%read_int("NLAY") & + ]) + if (.not. all(dis%idomain == idomain2d)) then + write (errmsg, fmtidomerr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + end select + case ('DISU') + select type (dis => this%dis) + type is (DisuType) + nodes = this%gfr%read_int("NODES") + if (nodes /= this%dis%nodes) then + write (errmsg, fmtdiserr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + idomain1d = this%gfr%read_int_1d("IDOMAIN") + if (.not. all(dis%idomain == idomain1d)) then + write (errmsg, fmtidomerr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + end select + case ('DIS2D') + select type (dis => this%dis) + type is (Dis2dType) + nodes = this%gfr%read_int("NCELLS") + if (nodes /= this%dis%nodes) then + write (errmsg, fmtdiserr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + idomain1d = this%gfr%read_int_1d("IDOMAIN") + idomain2d = reshape(idomain1d, [ & + this%gfr%read_int("NCOL"), & + this%gfr%read_int("NROW") & + ]) + if (.not. all(dis%idomain == idomain2d)) then + write (errmsg, fmtidomerr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + end select + case ('DISV2D') + select type (dis => this%dis) + type is (Disv2dType) + nodes = this%gfr%read_int("NODES") + if (nodes /= this%dis%nodes) then + write (errmsg, fmtdiserr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + idomain1d = this%gfr%read_int_1d("IDOMAIN") + if (.not. all(dis%idomain == idomain1d)) then + write (errmsg, fmtidomerr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + end select + case ('DISV1D') + select type (dis => this%dis) + type is (Disv1dType) + nodes = this%gfr%read_int("NCELLS") + if (nodes /= this%dis%nodes) then + write (errmsg, fmtdiserr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + idomain1d = this%gfr%read_int_1d("IDOMAIN") + if (.not. all(dis%idomain == idomain1d)) then + write (errmsg, fmtidomerr) trim(this%text) + call store_error(errmsg, terminate=.TRUE.) + end if + end select + end select + + if (allocated(idomain3d)) deallocate (idomain3d) + if (allocated(idomain2d)) deallocate (idomain2d) + if (allocated(idomain1d)) deallocate (idomain1d) + + call this%gfr%finalize() case default write (errmsg, '(a,3(1x,a))') & 'UNKNOWN', trim(adjustl(this%text)), 'PACKAGEDATA:', trim(keyword) @@ -421,25 +567,18 @@ subroutine read_packagedata(this) end subroutine read_packagedata !> @brief Initialize the budget file reader - !< subroutine initialize_bfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy integer(I4B) :: ncrbud - ! - ! -- Initialize the budget file reader call this%bfr%initialize(this%iubud, this%iout, ncrbud) - ! - ! -- todo: need to run through the budget terms - ! and do some checking + ! todo: need to run through the budget terms + ! and do some checking end subroutine initialize_bfr !> @brief Advance the budget file reader !! !! Advance the budget file reader by reading the next chunk !! of information for the current time step and stress period. - !! !< subroutine advance_bfr(this) ! -- modules @@ -586,35 +725,22 @@ subroutine advance_bfr(this) end subroutine advance_bfr !> @brief Finalize the budget file reader - !< subroutine finalize_bfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy - ! - ! -- Finalize the budget file reader call this%bfr%finalize() - ! end subroutine finalize_bfr !> @brief Initialize the head file reader - !< subroutine initialize_hfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy - ! - ! -- Initialize the budget file reader call this%hfr%initialize(this%iuhds, this%iout) - ! - ! -- todo: need to run through the head terms - ! and do some checking + ! todo: need to run through the head terms + ! and do some checking end subroutine initialize_hfr !> @brief Advance the head file reader - !< subroutine advance_hfr(this) - ! -- modules + ! modules use TdisModule, only: kstp, kper class(FlowModelInterfaceType) :: this integer(I4B) :: nu, nr, i, ilay @@ -705,22 +831,15 @@ subroutine advance_hfr(this) end subroutine advance_hfr !> @brief Finalize the head file reader - !< subroutine finalize_hfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy - ! - ! -- Finalize the head file reader close (this%iuhds) - ! end subroutine finalize_hfr !> @brief Initialize gwf terms from budget file !! !! initialize terms and figure out how many !! different terms and packages are contained within the file - !! !< subroutine initialize_gwfterms_from_bfr(this) ! -- modules @@ -818,7 +937,6 @@ subroutine initialize_gwfterms_from_bfr(this) end subroutine initialize_gwfterms_from_bfr !> @brief Initialize gwf terms from a GWF exchange - !< subroutine initialize_gwfterms_from_gwfbndlist(this) ! -- modules use BndModule, only: BndType, GetBndFromList @@ -918,22 +1036,16 @@ subroutine allocate_gwfpackages(this, ngwfterms) end subroutine allocate_gwfpackages !> @brief Deallocate memory in the gwfpackages array - !< subroutine deallocate_gwfpackages(this) - ! -- modules - ! -- dummy class(FlowModelInterfaceType) :: this - ! -- local integer(I4B) :: n - ! - ! -- initialize + do n = 1, this%nflowpack call this%gwfpackages(n)%da() end do end subroutine deallocate_gwfpackages !> @brief Find the package index for the package with the given name - !< subroutine get_package_index(this, name, idx) use BndModule, only: BndType, GetBndFromList class(FlowModelInterfaceType) :: this diff --git a/src/Model/TransportModel/tsp-fmi.f90 b/src/Model/TransportModel/tsp-fmi.f90 index bd696f98fab..caed1abf753 100644 --- a/src/Model/TransportModel/tsp-fmi.f90 +++ b/src/Model/TransportModel/tsp-fmi.f90 @@ -351,6 +351,7 @@ subroutine gwtfmi_da(this) call mem_deallocate(this%iubud) call mem_deallocate(this%iuhds) call mem_deallocate(this%iumvr) + call mem_deallocate(this%iugrb) call mem_deallocate(this%nflowpack) call mem_deallocate(this%idryinactive) ! diff --git a/src/Utilities/GridFileReader.f90 b/src/Utilities/GridFileReader.f90 new file mode 100644 index 00000000000..6ac8bcd9424 --- /dev/null +++ b/src/Utilities/GridFileReader.f90 @@ -0,0 +1,350 @@ +module GridFileReaderModule + + use KindModule + use SimModule, only: store_error, store_error_unit + use SimVariablesModule, only: errmsg + use ConstantsModule, only: LINELENGTH + use InputOutputModule, only: urword, openfile + use HashTableModule, only: HashTableType, hash_table_cr, hash_table_da + use ArrayHandlersModule, only: ExpandArray + + implicit none + + public :: GridFileReaderType + + type :: GridFileReaderType + private + integer(I4B), public :: inunit !< file unit + ! header + character(len=10), public :: grid_type !< DIS, DISV, DISU, etc + integer(I4B), public :: version !< binary grid file format version + integer(I4B) :: ntxt !< number of variables + integer(I4B) :: lentxt !< header line length per variable + ! index + type(HashTableType), pointer :: dim !< map variable name to number of dims + type(HashTableType), pointer :: pos !< map variable name to position in file + type(HashTableType), pointer :: typ !< map variable name to type (1=int, 2=dbl) + type(HashTableType), pointer :: shp_idx !< map variable name to index in shp + integer(I4B), allocatable :: shp(:) !< flat array of variable shapes + character(len=10), allocatable, public :: keys(:) !< variable names + contains + procedure, public :: initialize + procedure, public :: finalize + procedure, public :: read_int + procedure, public :: read_dbl + procedure, public :: read_int_1d + procedure, public :: read_dbl_1d + procedure, public :: read_grid_shape + procedure, private :: read_header + procedure, private :: read_header_meta + procedure, private :: read_header_body + end type GridFileReaderType + +contains + + !> @Brief Initialize the grid file reader. + subroutine initialize(this, iu) + class(GridFileReaderType) :: this + integer(I4B), intent(in) :: iu + + this%inunit = iu + call hash_table_cr(this%dim) + call hash_table_cr(this%pos) + call hash_table_cr(this%typ) + call hash_table_cr(this%shp_idx) + allocate (this%shp(0)) + call this%read_header() + + end subroutine initialize + + !> @brief Finalize the grid file reader. + subroutine finalize(this) + class(GridFileReaderType) :: this + + close (this%inunit) + call hash_table_da(this%dim) + call hash_table_da(this%pos) + call hash_table_da(this%typ) + call hash_table_da(this%shp_idx) + deallocate (this%shp) + + end subroutine finalize + + !> @brief Read the file's self-describing header. Internal use only. + subroutine read_header(this) + class(GridFileReaderType) :: this + call this%read_header_meta() + call this%read_header_body() + end subroutine read_header + + !> @brief Read self-describing metadata (first four lines). Internal use only. + subroutine read_header_meta(this) + ! dummy + class(GridFileReaderType) :: this + ! local + character(len=50) :: line + integer(I4B) :: lloc, istart, istop + integer(I4B) :: ival + real(DP) :: rval + + ! grid type + read (this%inunit) line + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, 0, 0) + if (line(istart:istop) /= 'GRID') then + call store_error('Binary grid file must begin with "GRID". '//& + &'Found: '//line(istart:istop)) + call store_error_unit(this%inunit) + end if + call urword(line, lloc, istart, istop, 1, ival, rval, 0, 0) + this%grid_type = line(istart:istop) + + ! version + read (this%inunit) line + lloc = 1 + call urword(line, lloc, istart, istop, 0, ival, rval, 0, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, 0, 0) + this%version = ival + + ! ntxt + read (this%inunit) line + lloc = 1 + call urword(line, lloc, istart, istop, 0, ival, rval, 0, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, 0, 0) + this%ntxt = ival + + ! lentxt + read (this%inunit) line + lloc = 1 + call urword(line, lloc, istart, istop, 0, ival, rval, 0, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, 0, 0) + this%lentxt = ival + + end subroutine read_header_meta + + !> @brief Read the header body section (text following first + !< four "meta" lines) and build an index. Internal use only. + subroutine read_header_body(this) + ! dummy + class(GridFileReaderType) :: this + ! local + character(len=:), allocatable :: body + character(len=:), allocatable :: line + character(len=10) :: key, dtype + real(DP) :: rval + integer(I4B) :: i, lloc, istart, istop, ival, pos + integer(I4B) :: nvars, ndim, dim, ishp + integer(I4B), allocatable :: shp(:) + + allocate (this%keys(this%ntxt)) + allocate (character(len=this%lentxt*this%ntxt) :: body) + allocate (character(len=this%lentxt) :: line) + + nvars = 0 + read (this%inunit) body + inquire (this%inunit, pos=pos) + do i = 1, this%lentxt * this%ntxt, this%lentxt + line = body(i:i + this%lentxt - 1) + lloc = 1 + + ! key + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, 0, 0) + key = line(istart:istop) + nvars = nvars + 1 + this%keys(nvars) = key + + ! type + call urword(line, lloc, istart, istop, 1, ival, rval, 0, 0) + dtype = line(istart:istop) + if (dtype == "INTEGER") then + call this%typ%add(key, 1) + else if (dtype == "DOUBLE") then + call this%typ%add(key, 2) + end if + + ! dims + call urword(line, lloc, istart, istop, 0, ival, rval, 0, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, 0, 0) + ndim = ival + call this%dim%add(key, ndim) + + ! shape + if (allocated(shp)) deallocate (shp) + allocate (shp(ndim)) + if (ndim > 0) then + do dim = 1, ndim + call urword(line, lloc, istart, istop, 2, ival, rval, 0, 0) + shp(dim) = ival + end do + ishp = size(this%shp) + call ExpandArray(this%shp, increment=ndim) + this%shp(ishp + 1:ishp + ndim) = shp + call this%shp_idx%add(key, ishp + 1) + end if + + ! position + call this%pos%add(key, pos) + if (ndim == 0) then + if (dtype == "INTEGER") then + pos = pos + 4 + else if (dtype == "DOUBLE") then + pos = pos + 8 + end if + else + if (dtype == "INTEGER") then + pos = pos + (product(shp) * 4) + else if (dtype == "DOUBLE") then + pos = pos + (product(shp) * 8) + end if + end if + end do + + rewind (this%inunit) + + end subroutine read_header_body + + !> @brief Read an integer scalar from a grid file. + function read_int(this, key) result(v) + class(GridFileReaderType), intent(inout) :: this + character(len=*), intent(in) :: key + integer(I4B) :: v + ! local + integer(I4B) :: ndim, pos, typ + character(len=:), allocatable :: msg + + msg = 'Variable '//trim(key)//' is not an integer scalar' + ndim = this%dim%get(key) + if (ndim /= 0) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + typ = this%typ%get(key) + if (typ /= 1) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + pos = this%pos%get(key) + read (this%inunit, pos=pos) v + rewind (this%inunit) + + end function read_int + + !> @brief Read a double precision scalar from a grid file. + function read_dbl(this, key) result(v) + class(GridFileReaderType), intent(inout) :: this + character(len=*), intent(in) :: key + real(DP) :: v + ! local + integer(I4B) :: ndim, pos, typ + character(len=:), allocatable :: msg + + msg = 'Variable '//trim(key)//' is not a double precision scalar' + ndim = this%dim%get(key) + if (ndim /= 0) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + typ = this%typ%get(key) + if (typ /= 2) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + pos = this%pos%get(key) + read (this%inunit, pos=pos) v + rewind (this%inunit) + + end function read_dbl + + !> @brief Read a 1D integer array from a grid file. + function read_int_1d(this, key) result(v) + class(GridFileReaderType), intent(inout) :: this + character(len=*), intent(in) :: key + integer(I4B), allocatable :: v(:) + ! local + integer(I4B) :: idx, ndim, nvals, pos, typ + character(len=:), allocatable :: msg + + msg = 'Variable '//trim(key)//' is not a 1D integer array' + ndim = this%dim%get(key) + if (ndim /= 1) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + typ = this%typ%get(key) + if (typ /= 1) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + idx = this%shp_idx%get(key) + pos = this%pos%get(key) + nvals = this%shp(idx) + allocate (v(nvals)) + read (this%inunit, pos=pos) v + rewind (this%inunit) + + end function read_int_1d + + !> @brief Read a 1D double array from a grid file. + function read_dbl_1d(this, key) result(v) + class(GridFileReaderType), intent(inout) :: this + character(len=*), intent(in) :: key + real(DP), allocatable :: v(:) + ! local + integer(I4B) :: idx, ndim, nvals, pos, typ + character(len=:), allocatable :: msg + + msg = 'Variable '//trim(key)//' is not a 1D double array' + ndim = this%dim%get(key) + if (ndim /= 1) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + typ = this%typ%get(key) + if (typ /= 2) then + write (errmsg, '(a)') msg + call store_error(errmsg, terminate=.TRUE.) + end if + idx = this%shp_idx%get(key) + pos = this%pos%get(key) + nvals = this%shp(idx) + allocate (v(nvals)) + read (this%inunit, pos=pos) v + rewind (this%inunit) + + end function read_dbl_1d + + !> @brief Read the grid shape from a grid file. + function read_grid_shape(this) result(v) + ! dummy + class(GridFileReaderType) :: this + integer(I4B), allocatable :: v(:) + + select case (this%grid_type) + case ("DIS") + allocate (v(3)) + v(1) = this%read_int("NLAY") + v(2) = this%read_int("NROW") + v(3) = this%read_int("NCOL") + case ("DISV") + allocate (v(2)) + v(1) = this%read_int("NLAY") + v(2) = this%read_int("NCPL") + case ("DISU") + allocate (v(1)) + v(1) = this%read_int("NODES") + case ("DIS2D") + allocate (v(2)) + v(1) = this%read_int("NROW") + v(2) = this%read_int("NCOL") + case ("DISV2D") + allocate (v(1)) + v(1) = this%read_int("NODES") + case ("DISV1D") + allocate (v(1)) + v(1) = this%read_int("NCELLS") + end select + + end function read_grid_shape + +end module GridFileReaderModule diff --git a/src/meson.build b/src/meson.build index 3049d11692b..bf059f34b6c 100644 --- a/src/meson.build +++ b/src/meson.build @@ -387,6 +387,7 @@ modflow_sources = files( 'Utilities' / 'DevFeature.f90', 'Utilities' / 'ErrorUtil.f90', 'Utilities' / 'GeomUtil.f90', + 'Utilities' / 'GridFileReader.f90', 'Utilities' / 'HashTable.f90', 'Utilities' / 'HeadFileReader.f90', 'Utilities' / 'HGeoUtil.f90', diff --git a/utils/mf5to6/make/makefile b/utils/mf5to6/make/makefile index 0d16924fd33..90f577019cb 100644 --- a/utils/mf5to6/make/makefile +++ b/utils/mf5to6/make/makefile @@ -1,4 +1,4 @@ -# makefile created by pymake (version 1.2.11.dev0) for the 'mf5to6' executable. +# makefile created by pymake (version 1.4.0.dev0) for the 'mf5to6' executable. include ./makedefaults diff --git a/utils/zonebudget/make/makefile b/utils/zonebudget/make/makefile index 4cb973e9759..c3eb5789bc3 100644 --- a/utils/zonebudget/make/makefile +++ b/utils/zonebudget/make/makefile @@ -1,4 +1,4 @@ -# makefile created by pymake (version 1.2.10.dev0) for the 'zbud6' executable. +# makefile created by pymake (version 1.4.0.dev0) for the 'zbud6' executable. include ./makedefaults @@ -34,8 +34,9 @@ $(OBJDIR)/BlockParser.o \ $(OBJDIR)/ArrayReaders.o \ $(OBJDIR)/zone.o \ $(OBJDIR)/Budget.o \ +$(OBJDIR)/HashTable.o \ $(OBJDIR)/zoneoutput.o \ -$(OBJDIR)/grb.o \ +$(OBJDIR)/GridFileReader.o \ $(OBJDIR)/budgetdata.o \ $(OBJDIR)/MathUtil.o \ $(OBJDIR)/GeomUtil.o \ diff --git a/utils/zonebudget/meson.build b/utils/zonebudget/meson.build index 1c9d4fb3e43..b147c5260f9 100644 --- a/utils/zonebudget/meson.build +++ b/utils/zonebudget/meson.build @@ -1,6 +1,5 @@ zonebudget_sources = files( 'src' / 'budgetdata.f90', - 'src' / 'grb.f90', 'src' / 'zbud6.f90', 'src' / 'zone.f90', 'src' / 'zoneoutput.f90', diff --git a/utils/zonebudget/msvs/zonebudget.vfproj b/utils/zonebudget/msvs/zonebudget.vfproj index 89a64a5816f..52f0b30bd5c 100644 --- a/utils/zonebudget/msvs/zonebudget.vfproj +++ b/utils/zonebudget/msvs/zonebudget.vfproj @@ -44,6 +44,8 @@ + + @@ -58,7 +60,6 @@ - diff --git a/utils/zonebudget/pymake/extrafiles.txt b/utils/zonebudget/pymake/extrafiles.txt index d33ec55b1a4..58d26f0ca9e 100644 --- a/utils/zonebudget/pymake/extrafiles.txt +++ b/utils/zonebudget/pymake/extrafiles.txt @@ -19,3 +19,5 @@ ../../../src/Utilities/version.f90 ../../../src/Utilities/DevFeature.f90 ../../../src/Utilities/Message.f90 +../../../src/Utilities/GridFileReader.f90 +../../../src/Utilities/HashTable.f90 diff --git a/utils/zonebudget/src/grb.f90 b/utils/zonebudget/src/grb.f90 deleted file mode 100644 index ec205703e03..00000000000 --- a/utils/zonebudget/src/grb.f90 +++ /dev/null @@ -1,168 +0,0 @@ -module GrbModule - - use KindModule - use SimVariablesModule, only: iout - use SimModule, only: store_error, store_error_unit, ustop - implicit none - private - public :: read_grb - -contains - - subroutine read_grb(inunit, ia, ja, mshape) -! ****************************************************************************** -! read_grb -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use InputOutputModule, only: urword - ! -- dummy - integer(I4B), intent(in) :: inunit - integer(I4B), allocatable, dimension(:), intent(out) :: ia - integer(I4B), allocatable, dimension(:), intent(out) :: ja - integer(I4B), allocatable, dimension(:), intent(out) :: mshape - ! -- local - character(len=50) :: hdrtxt - integer(I4B) :: lloc, istart, istop - character(len=50) :: dataname - character(len=50) :: datatype - integer(I4B) :: ntxt, lentxt, ndim, i, j, n, nval - integer(I4B) :: nja, ncells - real(DP) :: r, d - character(len=:), allocatable :: line - character(len=:), allocatable :: dfntxt - integer(I4B), dimension(:), allocatable :: ishape - integer(I4B), dimension(:), allocatable :: itmp - real(DP), dimension(:), allocatable :: dtmp -! ------------------------------------------------------------------------------ - ! - ! -- message - write (iout, '(/,a)') 'Processing Binary Grid File' - ! -- grid keyword - read (inunit) hdrtxt - lloc = 1 - call urword(hdrtxt, lloc, istart, istop, 0, i, r, iout, inunit) - if (hdrtxt(istart:istop) /= 'GRID') then - call store_error('GRB FILE MUST BEGIN WITH WORD GRID. FOUND: '// & - hdrtxt(istart:istop)) - call store_error_unit(inunit) - end if - ! - ! -- grid type, allocate mshape accordingly - call urword(hdrtxt, lloc, istart, istop, 0, i, r, iout, inunit) - if (hdrtxt(istart:istop) == 'DIS') then - write (iout, '(2x, a)') 'Detected regular MODFLOW grid (DIS)' - allocate (mshape(3)) - elseif (hdrtxt(istart:istop) == 'DISV') then - write (iout, '(2x, a)') 'Detected Discretization by Vertices grid (DISV)' - allocate (mshape(2)) - elseif (hdrtxt(istart:istop) == 'DISU') then - write (iout, '(2x, a)') 'Detected unstructured grid (DISU)' - allocate (mshape(1)) - else - call store_error('UNKNOWN GRID TYPE IN GRB FILE: '//hdrtxt(istart:istop)) - call store_error_unit(inunit) - end if - mshape(:) = 0 - ! - ! -- version - read (inunit) hdrtxt - write (iout, '(2x, a, a)') 'Detected ', trim(hdrtxt(1:49)) - ! - ! -- ntxt - read (inunit) hdrtxt - write (iout, '(2x, a, a)') 'Detected ', trim(hdrtxt(1:49)) - lloc = 1 - call urword(hdrtxt, lloc, istart, istop, 0, i, r, iout, inunit) - call urword(hdrtxt, lloc, istart, istop, 2, ntxt, r, iout, inunit) - - ! -- lentxt - read (inunit) hdrtxt - write (iout, '(2x, a, a)') 'Detected ', trim(hdrtxt(1:49)) - lloc = 1 - call urword(hdrtxt, lloc, istart, istop, 0, i, r, iout, inunit) - call urword(hdrtxt, lloc, istart, istop, 2, lentxt, r, iout, inunit) - ! - ! -- read txt definitions - allocate (character(len=lentxt) :: line) - allocate (character(len=lentxt*ntxt) :: dfntxt) - read (inunit) dfntxt - ! -- read each data record - do n = 1, lentxt * ntxt, lentxt - line = dfntxt(n:n + lentxt - 1) - lloc = 1 - call urword(line, lloc, istart, istop, 0, i, r, iout, inunit) - dataname = line(istart:istop) - call urword(line, lloc, istart, istop, 0, i, r, iout, inunit) - datatype = line(istart:istop) - call urword(line, lloc, istart, istop, 0, i, r, iout, inunit) - call urword(line, lloc, istart, istop, 2, ndim, r, iout, inunit) - allocate (ishape(ndim)) - do j = 1, ndim - call urword(line, lloc, istart, istop, 2, ishape(j), r, iout, inunit) - end do - select case (trim(datatype)) - case ('INTEGER') - if (ndim == 0) then - read (inunit) i - write (iout, '(2x, a, a, a, i0)') & - 'Detected ', trim(dataname), ' = ', i - if (trim(dataname) == 'NLAY') mshape(1) = i - if (trim(dataname) == 'NROW') mshape(2) = i - if (trim(dataname) == 'NCOL') mshape(3) = i - if (trim(dataname) == 'NCPL') mshape(2) = i - if (trim(dataname) == 'NJA') nja = i - if (trim(dataname) == 'NCELLS') ncells = i - if (trim(dataname) == 'NODES') then - ncells = i - mshape(1) = i - end if - else - write (iout, '(2x, a, a)') 'Detected integer array ', trim(dataname) - nval = 1 - do j = 1, ndim - nval = nval * ishape(j) - end do - allocate (itmp(nval)) - read (inunit) itmp - if (trim(dataname) == 'IA') then - allocate (ia(ncells + 1)) - ia = itmp - elseif (trim(dataname) == 'JA') then - allocate (ja(nja)) - ja = itmp - end if - deallocate (itmp) - end if - case ('DOUBLE') - if (ndim == 0) then - read (inunit) d - write (iout, '(2x, a, a, a, G0)') & - 'Detected ', trim(dataname), ' = ', d - else - write (iout, '(2x, a, a)') 'Detected double array ', trim(dataname) - nval = 1 - do j = 1, ndim - nval = nval * ishape(j) - end do - allocate (dtmp(nval)) - read (inunit) dtmp - deallocate (dtmp) - end if - end select - deallocate (ishape) - end do - close (inunit) - write (iout, '(a)') 'Done processing Binary Grid File' - ! - ! -- deallocate local storage - deallocate (line) - deallocate (dfntxt) - ! - ! -- return - return - end subroutine read_grb - -end module GrbModule diff --git a/utils/zonebudget/src/zbud6.f90 b/utils/zonebudget/src/zbud6.f90 index 29a5c98e85e..784cb850fd4 100644 --- a/utils/zonebudget/src/zbud6.f90 +++ b/utils/zonebudget/src/zbud6.f90 @@ -158,7 +158,7 @@ subroutine process_budget(iunit_csv, iunit_bud, iunit_zon, iunit_grb) zone_finalize, nmznfl, vbvl, vbznfl, maxzone use ZoneOutputModule, only: zoneoutput_init, zoneoutput_write, & zoneoutput_finalize - use GrbModule, only: read_grb + use GridFileReaderModule, only: GridFileReaderType use MessageModule, only: write_message, write_message_centered implicit none ! -- dummy @@ -181,6 +181,7 @@ subroutine process_budget(iunit_csv, iunit_bud, iunit_zon, iunit_grb) logical :: success logical :: hasiaja = .false. logical :: foundchd = .false. + type(GridFileReaderType) :: gfr ! ------------------------------------------------------------------------------ ! ! -- initialize local variables @@ -195,7 +196,11 @@ subroutine process_budget(iunit_csv, iunit_bud, iunit_zon, iunit_grb) inquire (unit=iunit_grb, opened=opengrb) if (opengrb) then hasiaja = .true. - call read_grb(iunit_grb, ia, ja, mshape) + call gfr%initialize(iunit_grb) + mshape = gfr%read_grid_shape() + ia = gfr%read_int_1d("IA") + ja = gfr%read_int_1d("JA") + call gfr%finalize() ncrgrb = size(ia) - 1 else errmsg = 'BUDGET FILE HAS "FLOW-JA-FACE" RECORD BUT NO GRB FILE SPECIFIED.'