From efafa5c3bd900bad23d4cb3e9259f0e03bd114b8 Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Fri, 8 Sep 2023 15:29:53 -0400 Subject: [PATCH] refactor(fmi): introduce FlowModelInterface module/type, use for gwt (#1332) --- make/makefile | 9 +- msvs/mf6core.vfproj | 1 + src/Model/Connection/GwtInterfaceModel.f90 | 2 +- src/Model/GroundWaterTransport/gwt1.f90 | 225 ++-- src/Model/GroundWaterTransport/gwt1fmi1.f90 | 650 +----------- .../ModelUtilities/FlowModelInterface.f90 | 965 ++++++++++++++++++ src/meson.build | 1 + 7 files changed, 1098 insertions(+), 755 deletions(-) create mode 100644 src/Model/ModelUtilities/FlowModelInterface.f90 diff --git a/make/makefile b/make/makefile index 8f85b9062aa..f2e1a58dd0b 100644 --- a/make/makefile +++ b/make/makefile @@ -125,17 +125,19 @@ $(OBJDIR)/PackageMover.o \ $(OBJDIR)/Obs3.o \ $(OBJDIR)/NumericalPackage.o \ $(OBJDIR)/Budget.o \ +$(OBJDIR)/BudgetTerm.o \ $(OBJDIR)/sort.o \ $(OBJDIR)/SfrCrossSectionUtils.o \ -$(OBJDIR)/BudgetTerm.o \ $(OBJDIR)/VirtualBase.o \ $(OBJDIR)/STLVecInt.o \ $(OBJDIR)/BoundaryPackage.o \ $(OBJDIR)/BaseModel.o \ $(OBJDIR)/InputDefinition.o \ +$(OBJDIR)/PackageBudget.o \ +$(OBJDIR)/HeadFileReader.o \ +$(OBJDIR)/BudgetObject.o \ $(OBJDIR)/SfrCrossSectionManager.o \ $(OBJDIR)/dag_module.o \ -$(OBJDIR)/BudgetObject.o \ $(OBJDIR)/VirtualDataLists.o \ $(OBJDIR)/VirtualDataContainer.o \ $(OBJDIR)/SimStages.o \ @@ -151,8 +153,7 @@ $(OBJDIR)/gwf3idm.o \ $(OBJDIR)/gwf3disv8idm.o \ $(OBJDIR)/gwf3disu8idm.o \ $(OBJDIR)/gwf3dis8idm.o \ -$(OBJDIR)/PackageBudget.o \ -$(OBJDIR)/HeadFileReader.o \ +$(OBJDIR)/FlowModelInterface.o \ $(OBJDIR)/PrintSaveManager.o \ $(OBJDIR)/Xt3dAlgorithm.o \ $(OBJDIR)/gwf3tvbase8.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 829cdc0728d..ebde0a6be87 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -187,6 +187,7 @@ + diff --git a/src/Model/Connection/GwtInterfaceModel.f90 b/src/Model/Connection/GwtInterfaceModel.f90 index 43c1d25865c..81d08b1a064 100644 --- a/src/Model/Connection/GwtInterfaceModel.f90 +++ b/src/Model/Connection/GwtInterfaceModel.f90 @@ -127,7 +127,7 @@ subroutine gwtifmod_df(this) ! define DISU disPtr => this%dis call this%gridConnection%getDiscretization(CastAsDisuType(disPtr)) - call this%fmi%fmi_df(this%dis, 0) + call this%fmi%fmi_df(this%dis) if (this%inadv > 0) then call this%adv%adv_df(adv_options) diff --git a/src/Model/GroundWaterTransport/gwt1.f90 b/src/Model/GroundWaterTransport/gwt1.f90 index 0e391664630..29bc7bbadce 100644 --- a/src/Model/GroundWaterTransport/gwt1.f90 +++ b/src/Model/GroundWaterTransport/gwt1.f90 @@ -71,7 +71,6 @@ module GwtModule procedure :: model_ot => gwt_ot procedure :: model_da => gwt_da procedure :: model_bdentry => gwt_bdentry - procedure :: allocate_scalars procedure, private :: package_create procedure, private :: ftype_check @@ -89,13 +88,8 @@ module GwtModule contains + !> @brief Create a new groundwater transport model object subroutine gwt_cr(filename, id, modelname) -! ****************************************************************************** -! gwt_cr -- Create a new groundwater transport model object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ListsModule, only: basemodellist use BaseModelModule, only: AddBaseModelToList @@ -115,15 +109,14 @@ subroutine gwt_cr(filename, id, modelname) character(len=LENMEMPATH) :: input_mempath character(len=LINELENGTH) :: lst_fname type(GwfNamParamFoundType) :: found - ! -- format -! ------------------------------------------------------------------------------ ! - ! -- Allocate a new GWT Model (this) and add it to basemodellist + ! -- Allocate a new GWT Model (this) allocate (this) ! ! -- Set memory path before allocation in memory manager can be done this%memoryPath = create_mem_path(modelname) ! + ! -- Allocate scalars and add model to basemodellist call this%allocate_scalars(modelname) model => this call AddBaseModelToList(basemodellist, model) @@ -168,27 +161,25 @@ subroutine gwt_cr(filename, id, modelname) return end subroutine gwt_cr + !> @brief Define packages of the model + ! + ! (1) call df routines for each package + ! (2) set variables and pointers + ! + !< subroutine gwt_df(this) -! ****************************************************************************** -! gwt_df -- Define packages of the model -! Subroutine: (1) call df routines for each package -! (2) set variables and pointers -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ModelPackageInputsModule, only: NIUNIT_GWT + use SimModule, only: store_error ! -- dummy class(GwtModelType) :: this ! -- local integer(I4B) :: ip class(BndType), pointer :: packobj -! ------------------------------------------------------------------------------ ! ! -- Define packages and utility objects call this%dis%dis_df() - call this%fmi%fmi_df(this%dis, this%inssm) + call this%fmi%fmi_df(this%dis) if (this%inmvt > 0) call this%mvt%mvt_df(this%dis) if (this%inadv > 0) call this%adv%adv_df() if (this%indsp > 0) call this%dsp%dsp_df(this%dis) @@ -196,6 +187,15 @@ subroutine gwt_df(this) call this%oc%oc_df() call this%budget%budget_df(NIUNIT_GWT, 'MASS', 'M') ! + ! -- Check for SSM package + if (this%inssm == 0) then + if (this%fmi%nflowpack > 0) then + call store_error('Flow model has boundary packages, but there & + &is no SSM package. The SSM package must be activated.', & + terminate=.TRUE.) + end if + end if + ! ! -- Assign or point model members to dis members this%neq = this%dis%nodes this%nja = this%dis%nja @@ -220,13 +220,9 @@ subroutine gwt_df(this) return end subroutine gwt_df + !> @brief Add the internal connections of this model to the sparse matrix + !< subroutine gwt_ac(this, sparse) -! ****************************************************************************** -! gwt_ac -- Add the internal connections of this model to the sparse matrix -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SparseModule, only: sparsematrix ! -- dummy @@ -235,7 +231,6 @@ subroutine gwt_ac(this, sparse) ! -- local class(BndType), pointer :: packobj integer(I4B) :: ip -! ------------------------------------------------------------------------------ ! ! -- Add the internal connections of this model to sparse call this%dis%dis_ac(this%moffset, sparse) @@ -252,21 +247,14 @@ subroutine gwt_ac(this, sparse) return end subroutine gwt_ac + !> @brief Map connection positions in numerical solution coefficient matrix. subroutine gwt_mc(this, matrix_sln) -! ****************************************************************************** -! gwt_mc -- Map the positions of this models connections in the -! numerical solution coefficient matrix. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwtModelType) :: this class(MatrixBaseType), pointer :: matrix_sln !< global system matrix ! -- local class(BndType), pointer :: packobj integer(I4B) :: ip -! ------------------------------------------------------------------------------ ! ! -- Find the position of each connection in the global ia, ja structure ! and store them in idxglo. @@ -283,15 +271,13 @@ subroutine gwt_mc(this, matrix_sln) return end subroutine gwt_mc + !> @brief Allocate and Read + ! + ! (1) allocates and reads packages part of this model, + ! (2) allocates memory for arrays part of this model object + ! + !< subroutine gwt_ar(this) -! ****************************************************************************** -! gwt_ar -- GroundWater Transport Model Allocate and Read -! Subroutine: (1) allocates and reads packages part of this model, -! (2) allocates memory for arrays part of this model object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: DHNOFLO ! -- dummy @@ -299,7 +285,6 @@ subroutine gwt_ar(this) ! -- locals integer(I4B) :: ip class(BndType), pointer :: packobj -! ------------------------------------------------------------------------------ ! ! -- Allocate and read modules attached to model call this%fmi%fmi_ar(this%ibound) @@ -331,14 +316,9 @@ subroutine gwt_ar(this) return end subroutine gwt_ar + !> @brief Read and prepare (calls package read and prepare routines) + !< subroutine gwt_rp(this) -! ****************************************************************************** -! gwt_rp -- GroundWater Transport Model Read and Prepare -! Subroutine: (1) calls package read and prepare routines -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use TdisModule, only: readnewdata ! -- dummy @@ -346,7 +326,6 @@ subroutine gwt_rp(this) ! -- local class(BndType), pointer :: packobj integer(I4B) :: ip -! ------------------------------------------------------------------------------ ! ! -- In fmi, check for mvt and mvrbudobj consistency call this%fmi%fmi_rp(this%inmvt) @@ -368,14 +347,9 @@ subroutine gwt_rp(this) return end subroutine gwt_rp + !> @brief Time step advance (calls package advance subroutines) + !< subroutine gwt_ad(this) -! ****************************************************************************** -! gwt_ad -- GroundWater Transport Model Time Step Advance -! Subroutine: (1) calls package advance subroutines -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SimVariablesModule, only: isimcheck, iFailedStepRetry ! -- dummy @@ -384,7 +358,6 @@ subroutine gwt_ad(this) ! -- local integer(I4B) :: irestore integer(I4B) :: ip, n -! ------------------------------------------------------------------------------ ! ! -- Reset state variable irestore = 0 @@ -429,21 +402,15 @@ subroutine gwt_ad(this) return end subroutine gwt_ad + !> @brief Calculate coefficients + !< subroutine gwt_cf(this, kiter) -! ****************************************************************************** -! gwt_cf -- GroundWater Transport Model calculate coefficients -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwtModelType) :: this integer(I4B), intent(in) :: kiter ! -- local class(BndType), pointer :: packobj integer(I4B) :: ip -! ------------------------------------------------------------------------------ ! ! -- Call package cf routines do ip = 1, this%bndlist%Count() @@ -455,14 +422,9 @@ subroutine gwt_cf(this, kiter) return end subroutine gwt_cf + !> @brief Fill coefficients + !< subroutine gwt_fc(this, kiter, matrix_sln, inwtflag) -! ****************************************************************************** -! gwt_fc -- GroundWater Transport Model fill coefficients -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwtModelType) :: this integer(I4B), intent(in) :: kiter @@ -471,7 +433,6 @@ subroutine gwt_fc(this, kiter, matrix_sln, inwtflag) ! -- local class(BndType), pointer :: packobj integer(I4B) :: ip -! ------------------------------------------------------------------------------ ! ! -- call fc routines call this%fmi%fmi_fc(this%dis%nodes, this%xold, this%nja, matrix_sln, & @@ -505,14 +466,9 @@ subroutine gwt_fc(this, kiter, matrix_sln, inwtflag) return end subroutine gwt_fc + !> @brief Final convergence check (calls package cc routines) + !< subroutine gwt_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) -! ****************************************************************************** -! gwt_cc -- GroundWater Transport Model Final Convergence Check -! Subroutine: (1) calls package cc routines -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwtModelType) :: this integer(I4B), intent(in) :: innertot @@ -523,32 +479,26 @@ subroutine gwt_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak) integer(I4B), intent(inout) :: ipak real(DP), intent(inout) :: dpak ! -- local - !class(BndType), pointer :: packobj - !integer(I4B) :: ip + ! class(BndType), pointer :: packobj + ! integer(I4B) :: ip ! -- formats -! ------------------------------------------------------------------------------ ! ! -- If mover is on, then at least 2 outers required if (this%inmvt > 0) call this%mvt%mvt_cc(kiter, iend, icnvgmod, cpak, dpak) ! ! -- Call package cc routines - !do ip = 1, this%bndlist%Count() - ! packobj => GetBndFromList(this%bndlist, ip) - ! call packobj%bnd_cc(iend, icnvg, hclose, rclose) - !enddo + ! do ip = 1, this%bndlist%Count() + ! packobj => GetBndFromList(this%bndlist, ip) + ! call packobj%bnd_cc(iend, icnvg, hclose, rclose) + ! enddo ! ! -- return return end subroutine gwt_cc + !> @brief Calculate intercell flows (flowja) + !< subroutine gwt_cq(this, icnvg, isuppress_output) -! ****************************************************************************** -! gwt_cq --Groundwater transport model calculate flow -! Subroutine: (1) Calculate intercell flows (flowja) -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SparseModule, only: csr_diagsum ! -- dummy @@ -559,7 +509,6 @@ subroutine gwt_cq(this, icnvg, isuppress_output) integer(I4B) :: i integer(I4B) :: ip class(BndType), pointer :: packobj -! ------------------------------------------------------------------------------ ! ! -- Construct the flowja array. Flowja is calculated each time, even if ! output is suppressed. (flowja is positive into a cell.) The diagonal @@ -594,15 +543,13 @@ subroutine gwt_cq(this, icnvg, isuppress_output) return end subroutine gwt_cq + !> @brief Model budget + ! + ! (1) Calculate intercell flows (flowja) + ! (2) Calculate package contributions to model budget + ! + !< subroutine gwt_bd(this, icnvg, isuppress_output) -! ****************************************************************************** -! gwt_bd --GroundWater Transport Model Budget -! Subroutine: (1) Calculate intercell flows (flowja) -! (2) Calculate package contributions to model budget -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ use ConstantsModule, only: DZERO ! -- dummy class(GwtModelType) :: this @@ -611,7 +558,6 @@ subroutine gwt_bd(this, icnvg, isuppress_output) ! -- local integer(I4B) :: ip class(BndType), pointer :: packobj -! ------------------------------------------------------------------------------ ! ! -- Save the solution convergence flag this%icnvg = icnvg @@ -635,13 +581,9 @@ subroutine gwt_bd(this, icnvg, isuppress_output) return end subroutine gwt_bd + !> @brief Print and/or save model output + !< subroutine gwt_ot(this) -! ****************************************************************************** -! gwt_ot -- GroundWater Transport Model Output -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use TdisModule, only: kstp, kper, tdis_ot, endofperiod ! -- dummy @@ -657,7 +599,6 @@ subroutine gwt_ot(this) character(len=*), parameter :: fmtnocnvg = & "(1X,/9X,'****FAILED TO MEET SOLVER CONVERGENCE CRITERIA IN TIME STEP ', & &I0,' OF STRESS PERIOD ',I0,'****')" -! ------------------------------------------------------------------------------ ! ! -- Set write and print flags idvsave = 0 @@ -700,6 +641,8 @@ subroutine gwt_ot(this) return end subroutine gwt_ot + !> @brief Calculate and save observations + !< subroutine gwt_ot_obs(this) class(GwtModelType) :: this class(BndType), pointer :: packobj @@ -718,6 +661,8 @@ subroutine gwt_ot_obs(this) end subroutine gwt_ot_obs + !> @brief Save flows + !< subroutine gwt_ot_flow(this, icbcfl, ibudfl, icbcun) class(GwtModelType) :: this integer(I4B), intent(in) :: icbcfl @@ -770,13 +715,9 @@ subroutine gwt_ot_flow(this, icbcfl, ibudfl, icbcun) end subroutine gwt_ot_flow + !> @brief Write intercell flows + !< subroutine gwt_ot_flowja(this, nja, flowja, icbcfl, icbcun) -! ****************************************************************************** -! gwt_ot_flowja -- Write intercell flows -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwtModelType) :: this integer(I4B), intent(in) :: nja @@ -786,7 +727,6 @@ subroutine gwt_ot_flowja(this, nja, flowja, icbcfl, icbcun) ! -- local integer(I4B) :: ibinun ! -- formats -! ------------------------------------------------------------------------------ ! ! -- Set unit number for binary output if (this%ipakcb < 0) then @@ -807,6 +747,8 @@ subroutine gwt_ot_flowja(this, nja, flowja, icbcfl, icbcun) return end subroutine gwt_ot_flowja + !> @brief Print dependent variables + !< subroutine gwt_ot_dv(this, idvsave, idvprint, ipflag) class(GwtModelType) :: this integer(I4B), intent(in) :: idvsave @@ -826,6 +768,8 @@ subroutine gwt_ot_dv(this, idvsave, idvprint, ipflag) end subroutine gwt_ot_dv + !> @brief Print budget summary + !< subroutine gwt_ot_bdsummary(this, ibudfl, ipflag) use TdisModule, only: kstp, kper, totim class(GwtModelType) :: this @@ -857,13 +801,9 @@ subroutine gwt_ot_bdsummary(this, ibudfl, ipflag) end subroutine gwt_ot_bdsummary + !> @brief Deallocate + !< subroutine gwt_da(this) -! ****************************************************************************** -! gwt_da -- Deallocate -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use MemoryManagerModule, only: mem_deallocate use MemoryManagerExtModule, only: memorylist_remove @@ -873,7 +813,6 @@ subroutine gwt_da(this) ! -- local integer(I4B) :: ip class(BndType), pointer :: packobj -! ------------------------------------------------------------------------------ ! ! -- Deallocate idm memory call memorylist_remove(this%name, 'NAM', idm_context) @@ -947,7 +886,6 @@ subroutine gwt_bdentry(this, budterm, budtxt, rowlabel) real(DP), dimension(:, :), intent(in) :: budterm character(len=LENBUDTXT), dimension(:), intent(in) :: budtxt character(len=*), intent(in) :: rowlabel -! ------------------------------------------------------------------------------ ! call this%budget%addentry(budterm, delt, budtxt, rowlabel=rowlabel) ! @@ -988,19 +926,14 @@ function gwt_get_iasym(this) result(iasym) return end function gwt_get_iasym + !> @brief Allocate memory for non-allocatable members + !< subroutine allocate_scalars(this, modelname) -! ****************************************************************************** -! allocate_scalars -- Allocate memory for non-allocatable members -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use MemoryManagerModule, only: mem_allocate ! -- dummy class(GwtModelType) :: this character(len=*), intent(in) :: modelname -! ------------------------------------------------------------------------------ ! ! -- allocate members from parent class call this%NumericalModelType%allocate_scalars(modelname) @@ -1030,14 +963,10 @@ subroutine allocate_scalars(this, modelname) return end subroutine allocate_scalars + !> @brief Create boundary condition packages for this model + !< subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & iout) -! ****************************************************************************** -! package_create -- Create boundary condition packages for this model -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: LINELENGTH use SimModule, only: store_error @@ -1062,7 +991,6 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & class(BndType), pointer :: packobj class(BndType), pointer :: packobj2 integer(I4B) :: ip -! ------------------------------------------------------------------------------ ! ! -- This part creates the package object select case (filtyp) @@ -1109,13 +1037,9 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & return end subroutine package_create + !> @brief Make sure required input files have been specified + !< subroutine ftype_check(this, indis) -! ****************************************************************************** -! ftype_check -- Check to make sure required input files have been specified -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: LINELENGTH use SimModule, only: store_error, count_errors, store_error_filename @@ -1124,7 +1048,6 @@ subroutine ftype_check(this, indis) integer(I4B), intent(in) :: indis ! -- local character(len=LINELENGTH) :: errmsg -! ------------------------------------------------------------------------------ ! ! -- Check for IC6, DIS(u), and MST. Stop if not present. if (this%inic == 0) then @@ -1284,7 +1207,7 @@ subroutine create_packages(this) mempath = mempaths(n) inunit => inunits(n) ! - ! -- create dis package as it is a prerequisite for other packages + ! -- create dis package first as it is a prerequisite for other packages select case (pkgtype) case ('DIS6') indis = 1 diff --git a/src/Model/GroundWaterTransport/gwt1fmi1.f90 b/src/Model/GroundWaterTransport/gwt1fmi1.f90 index cf5680e1868..cd1f5ef6bac 100644 --- a/src/Model/GroundWaterTransport/gwt1fmi1.f90 +++ b/src/Model/GroundWaterTransport/gwt1fmi1.f90 @@ -5,7 +5,7 @@ module GwtFmiModule LENPACKAGENAME use SimModule, only: store_error, store_error_unit use SimVariablesModule, only: errmsg - use NumericalPackageModule, only: NumericalPackageType + use FlowModelInterfaceModule, only: FlowModelInterfaceType use BaseDisModule, only: DisBaseType use ListModule, only: ListType use BudgetFileReaderModule, only: BudgetFileReaderType @@ -19,6 +19,8 @@ module GwtFmiModule public :: GwtFmiType public :: fmi_cr + character(len=LENPACKAGENAME) :: text = ' GWTFMI' + integer(I4B), parameter :: NBDITEMS = 2 character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt data budtxt/' FLOW-ERROR', ' FLOW-CORRECTION'/ @@ -32,97 +34,56 @@ module GwtFmiModule type(BudgetObjectType), pointer :: ptr end type BudObjPtrArray - type, extends(NumericalPackageType) :: GwtFmiType + type, extends(FlowModelInterfaceType) :: GwtFmiType - logical, pointer :: flows_from_file => null() !< if .false., then flows come from GWF through GWF-GWT exg integer(I4B), dimension(:), pointer, contiguous :: iatp => null() !< advanced transport package applied to gwfpackages - type(ListType), pointer :: gwfbndlist => null() !< list of gwf stress packages - integer(I4B), pointer :: iflowsupdated => null() !< flows were updated for this time step integer(I4B), pointer :: iflowerr => null() !< add the flow error correction real(DP), dimension(:), pointer, contiguous :: flowcorrect => null() !< mass flow correction - integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to GWT 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 - real(DP), dimension(:), pointer, contiguous :: gwfsat => null() !< pointer to the GWF saturation array - integer(I4B), dimension(:), pointer, contiguous :: ibdgwfsat0 => null() !< mark cells with saturation = 0 to exclude from dispersion - real(DP), dimension(:), pointer, contiguous :: gwfstrgss => null() !< pointer to flow model QSTOSS - real(DP), dimension(:), pointer, contiguous :: gwfstrgsy => null() !< pointer to flow model QSTOSY - integer(I4B), pointer :: igwfstrgss => null() !< indicates if gwfstrgss is available - integer(I4B), pointer :: igwfstrgsy => null() !< indicates if gwfstrgsy is available - 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 :: 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(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(DataAdvancedPackageType), & dimension(:), pointer, contiguous :: datp => null() - character(len=16), dimension(:), allocatable :: flowpacknamearray !< array of boundary package names (e.g. LAK-1, SFR-3, etc.) type(BudObjPtrArray), dimension(:), allocatable :: aptbudobj !< flow budget objects for the advanced packages contains - procedure :: fmi_df - procedure :: fmi_ar + procedure :: allocate_arrays => gwtfmi_allocate_arrays + procedure :: allocate_gwfpackages => gwtfmi_allocate_gwfpackages + procedure :: allocate_scalars => gwtfmi_allocate_scalars + procedure :: deallocate_gwfpackages => gwtfmi_deallocate_gwfpackages procedure :: fmi_rp procedure :: fmi_ad procedure :: fmi_fc procedure :: fmi_cq procedure :: fmi_bd procedure :: fmi_ot_flow - procedure :: fmi_da - procedure :: allocate_scalars - procedure :: allocate_arrays + procedure :: fmi_da => gwtfmi_da procedure :: gwfsatold - procedure :: read_options - procedure :: read_packagedata - procedure :: initialize_bfr - procedure :: advance_bfr - procedure :: finalize_bfr - procedure :: initialize_hfr - procedure :: advance_hfr - procedure :: finalize_hfr procedure :: initialize_gwfterms_from_bfr procedure :: initialize_gwfterms_from_gwfbndlist - procedure :: allocate_gwfpackages - procedure :: deallocate_gwfpackages - procedure :: get_package_index + procedure :: read_options => gwtfmi_read_options + procedure :: read_packagedata => gwtfmi_read_packagedata procedure :: set_aptbudobj_pointer end type GwtFmiType contains + !> @brief Create a new FMI object subroutine fmi_cr(fmiobj, name_model, inunit, iout) -! ****************************************************************************** -! fmi_cr -- Create a new FMI object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy type(GwtFmiType), pointer :: fmiobj character(len=*), intent(in) :: name_model integer(I4B), intent(in) :: inunit integer(I4B), intent(in) :: iout -! ------------------------------------------------------------------------------ ! ! -- Create the object allocate (fmiobj) ! ! -- create name and memory path call fmiobj%set_names(1, name_model, 'FMI', 'FMI') + fmiobj%text = text ! ! -- Allocate scalars call fmiobj%allocate_scalars() ! - ! -- if inunit == 0, then there is no file to read, but it still needs - ! to be active in order to manage pointers to gwf model - !if (inunit == 0) inunit = 1 - ! ! -- Set variables fmiobj%inunit = inunit fmiobj%iout = iout @@ -134,108 +95,8 @@ subroutine fmi_cr(fmiobj, name_model, inunit, iout) return end subroutine fmi_cr - subroutine fmi_df(this, dis, inssm) -! ****************************************************************************** -! fmi_df -- Define -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use SimModule, only: store_error - ! -- dummy - class(GwtFmiType) :: this - class(DisBaseType), pointer, intent(in) :: dis - integer(I4B), intent(in) :: inssm - ! -- local - ! -- formats - character(len=*), parameter :: fmtfmi = & - "(1x,/1x,'FMI -- FLOW MODEL INTERFACE, VERSION 1, 8/29/2017', & - &' INPUT READ FROM UNIT ', i0, //)" - character(len=*), parameter :: fmtfmi0 = & - &"(1x,/1x,'FMI -- FLOW MODEL INTERFACE, VERSION 1, 8/29/2017')" -! ------------------------------------------------------------------------------ - ! - ! --print a message identifying the FMI package. - if (this%iout > 0) then - if (this%inunit /= 0) then - write (this%iout, fmtfmi) this%inunit - else - write (this%iout, fmtfmi0) - if (this%flows_from_file) then - write (this%iout, '(a)') ' FLOWS ARE ASSUMED TO BE ZERO.' - else - write (this%iout, '(a)') ' FLOWS PROVIDED BY A GWF MODEL IN THIS & - &SIMULATION' - end if - end if - end if - ! - ! -- store pointers to arguments that were passed in - this%dis => dis - ! - ! -- Read fmi options - if (this%inunit /= 0) then - call this%read_options() - end if - ! - ! -- Read packagedata options - if (this%inunit /= 0 .and. this%flows_from_file) then - call this%read_packagedata() - call this%initialize_gwfterms_from_bfr() - end if - ! - ! -- If GWF-GWT exchange is active, then setup gwfterms from bndlist - if (.not. this%flows_from_file) then - call this%initialize_gwfterms_from_gwfbndlist() - end if - ! - ! -- Make sure that ssm is on if there are any boundary packages - if (inssm == 0) then - if (this%nflowpack > 0) then - call store_error('Flow model has boundary packages, but there & - &is no SSM package. The SSM package must be activated.', & - terminate=.TRUE.) - end if - end if - ! - ! -- Return - return - end subroutine fmi_df - - subroutine fmi_ar(this, ibound) -! ****************************************************************************** -! fmi_ar -- Allocate and Read -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use SimModule, only: store_error - ! -- dummy - class(GwtFmiType) :: this - integer(I4B), dimension(:), pointer, contiguous :: ibound - ! -- local - ! -- formats -! ------------------------------------------------------------------------------ - ! - ! -- store pointers to arguments that were passed in - this%ibound => ibound - ! - ! -- Allocate arrays - call this%allocate_arrays(this%dis%nodes) - ! - ! -- Return - return - end subroutine fmi_ar - + !> @brief Read and prepare subroutine fmi_rp(this, inmvr) -! ****************************************************************************** -! fmi_rp -- Read and prepare -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use TdisModule, only: kper, kstp ! -- dummy @@ -243,7 +104,6 @@ subroutine fmi_rp(this, inmvr) integer(I4B), intent(in) :: inmvr ! -- local ! -- formats -! ------------------------------------------------------------------------------ ! ! --Check to make sure MVT Package is active if mvr flows are available. ! This cannot be checked until RP because exchange doesn't set a pointer @@ -266,13 +126,8 @@ subroutine fmi_rp(this, inmvr) return end subroutine fmi_rp + !> @brief Advance subroutine fmi_ad(this, cnew) -! ****************************************************************************** -! fmi_ad -- advance -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: DHDRY ! -- dummy @@ -290,7 +145,6 @@ subroutine fmi_ad(this, cnew) character(len=*), parameter :: fmtrewet = & &"(/1X,'DRY CELL REACTIVATED AT ', a,& &' WITH STARTING CONCENTRATION =',G13.5)" -! ------------------------------------------------------------------------------ ! ! -- Set flag to indicated that flows are being updated. For the case where ! flows may be reused (only when flows are read from a file) then set @@ -526,7 +380,7 @@ subroutine fmi_ot_flow(this, icbcfl, icbcun) return end subroutine fmi_ot_flow - subroutine fmi_da(this) + subroutine gwtfmi_da(this) ! ****************************************************************************** ! fmi_da -- Deallocate variables ! ****************************************************************************** @@ -582,21 +436,15 @@ subroutine fmi_da(this) ! ! -- Return return - end subroutine fmi_da + end subroutine gwtfmi_da - subroutine allocate_scalars(this) -! ****************************************************************************** -! allocate_scalars -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Allocate scalars + subroutine gwtfmi_allocate_scalars(this) ! -- modules use MemoryManagerModule, only: mem_allocate, mem_setptr ! -- dummy class(GwtFmiType) :: this ! -- local -! ------------------------------------------------------------------------------ ! ! -- allocate scalars in NumericalPackageType call this%NumericalPackageType%allocate_scalars() @@ -629,15 +477,10 @@ subroutine allocate_scalars(this) ! ! -- Return return - end subroutine allocate_scalars + end subroutine gwtfmi_allocate_scalars - subroutine allocate_arrays(this, nodes) -! ****************************************************************************** -! allocate_arrays -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Allocate arrays + subroutine gwtfmi_allocate_arrays(this, nodes) use MemoryManagerModule, only: mem_allocate !modules use ConstantsModule, only: DZERO @@ -646,7 +489,6 @@ subroutine allocate_arrays(this, nodes) integer(I4B), intent(in) :: nodes ! -- local integer(I4B) :: n -! ------------------------------------------------------------------------------ ! ! -- Allocate variables needed for all cases if (this%iflowerr == 0) then @@ -707,7 +549,7 @@ subroutine allocate_arrays(this, nodes) ! ! -- Return return - end subroutine allocate_arrays + end subroutine gwtfmi_allocate_arrays function gwfsatold(this, n, delt) result(satold) ! ****************************************************************************** @@ -742,14 +584,8 @@ function gwfsatold(this, n, delt) result(satold) return end function gwfsatold - subroutine read_options(this) -! ****************************************************************************** -! read_options -- Read Options -! Subroutine: (1) read options from input file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Read options from input file + subroutine gwtfmi_read_options(this) ! -- modules use ConstantsModule, only: LINELENGTH, DEM6 use InputOutputModule, only: getunit, openfile, urdaux @@ -765,7 +601,6 @@ subroutine read_options(this) &WHENEVER ICBCFL IS NOT ZERO AND FLOW IMBALANCE CORRECTION ACTIVE.')" character(len=*), parameter :: fmtifc = & &"(4x,'MASS WILL BE ADDED OR REMOVED TO COMPENSATE FOR FLOW IMBALANCE.')" -! ------------------------------------------------------------------------------ ! ! -- get options block call this%parser%GetBlock('OPTIONS', isfound, ierr, blockRequired=.false., & @@ -797,16 +632,10 @@ subroutine read_options(this) ! ! -- return return - end subroutine read_options + end subroutine gwtfmi_read_options - subroutine read_packagedata(this) -! ****************************************************************************** -! read_packagedata -- Read PACKAGEDATA block -! Subroutine: (1) read packagedata block from input file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Read packagedata block from input file + subroutine gwtfmi_read_packagedata(this) ! -- modules use OpenSpecModule, only: ACCESS, FORM use ConstantsModule, only: LINELENGTH, DEM6, LENPACKAGENAME @@ -826,7 +655,6 @@ subroutine read_packagedata(this) logical :: blockrequired logical :: exist type(BudObjPtrArray), dimension(:), allocatable :: tmpbudobj -! ------------------------------------------------------------------------------ ! ! -- initialize iapt = 0 @@ -934,18 +762,17 @@ subroutine read_packagedata(this) ! ! -- return return - end subroutine read_packagedata - + end subroutine gwtfmi_read_packagedata + + !> @brief Set the pointer to a budget object + !! + !! An advanced transport can pass in a name and a + !! pointer budget object, and this routine will look through the budget + !! objects managed by FMI and point to the one with the same name, such as + !! LAK-1, SFR-1, etc. + !! + !< subroutine set_aptbudobj_pointer(this, name, budobjptr) -! ****************************************************************************** -! set_aptbudobj_pointer -- an advanced transport can pass in a name and a -! pointer budget object, and this routine will look through the budget -! objects managed by FMI and point to the one with the same name, such as -! LAK-1, SFR-1, etc. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules class(GwtFmiType) :: this ! -- dumm @@ -953,7 +780,6 @@ subroutine set_aptbudobj_pointer(this, name, budobjptr) type(BudgetObjectType), pointer :: budobjptr ! -- local integer(I4B) :: i -! ------------------------------------------------------------------------------ ! ! -- find and set the pointer do i = 1, size(this%aptbudobj) @@ -967,333 +793,8 @@ subroutine set_aptbudobj_pointer(this, name, budobjptr) return end subroutine set_aptbudobj_pointer - subroutine initialize_bfr(this) -! ****************************************************************************** -! initialize_bfr -- initalize the budget file reader -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - class(GwtFmiType) :: 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 - end subroutine initialize_bfr - - subroutine advance_bfr(this) -! ****************************************************************************** -! advance_bfr -- advance the budget file reader by reading the next chunk -! of information for the current time step and stress period -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use TdisModule, only: kstp, kper - ! -- dummy - class(GwtFmiType) :: this - ! -- local - logical :: success - integer(I4B) :: n - integer(I4B) :: ipos - integer(I4B) :: nu, nr - integer(I4B) :: ip, i - logical :: readnext - ! -- format - character(len=*), parameter :: fmtkstpkper = & - &"(1x,/1x,'FMI READING BUDGET TERMS FOR KSTP ', i0, ' KPER ', i0)" - character(len=*), parameter :: fmtbudkstpkper = & - "(1x,/1x, 'FMI SETTING BUDGET TERMS FOR KSTP ', i0, ' AND KPER ', & - &i0, ' TO BUDGET FILE TERMS FROM KSTP ', i0, ' AND KPER ', i0)" -! ------------------------------------------------------------------------------ - ! - ! -- If the latest record read from the budget file is from a stress - ! -- period with only one time step, reuse that record (do not read a - ! -- new record) if the GWT model is still in that same stress period, - ! -- or if that record is the last one in the budget file. - readnext = .true. - if (kstp * kper > 1) then - if (this%bfr%kstp == 1) then - if (this%bfr%kpernext == kper + 1) then - readnext = .false. - else if (this%bfr%endoffile) then - readnext = .false. - end if - else if (this%bfr%endoffile) then - write (errmsg, '(a)') 'Reached end of GWF budget & - &file before reading sufficient budget information for this & - &GWT simulation.' - call store_error(errmsg) - call store_error_unit(this%iubud) - end if - end if - ! - ! -- Read the next record - if (readnext) then - ! - ! -- Write the current time step and stress period - write (this%iout, fmtkstpkper) kstp, kper - ! - ! -- loop through the budget terms for this stress period - ! i is the counter for gwf flow packages - ip = 1 - do n = 1, this%bfr%nbudterms - call this%bfr%read_record(success, this%iout) - if (.not. success) then - write (errmsg, '(a)') 'GWF budget read not successful' - call store_error(errmsg) - call store_error_unit(this%iubud) - end if - ! - ! -- Ensure kper is same between model and budget file - if (kper /= this%bfr%kper) then - write (errmsg, '(a)') 'Period number in budget file & - &does not match period number in transport model. If there & - &is more than one time step in the budget file for a given stress & - &period, budget file time steps must match GWT model time steps & - &one-for-one in that stress period.' - call store_error(errmsg) - call store_error_unit(this%iubud) - end if - ! - ! -- if budget file kstp > 1, then kstp must match - if (this%bfr%kstp > 1 .and. (kstp /= this%bfr%kstp)) then - write (errmsg, '(a)') 'Time step number in budget file & - &does not match time step number in transport model. If there & - &is more than one time step in the budget file for a given stress & - &period, budget file time steps must match gwt model time steps & - &one-for-one in that stress period.' - call store_error(errmsg) - call store_error_unit(this%iubud) - end if - ! - ! -- parse based on the type of data, and compress all user node - ! numbers into reduced node numbers - select case (trim(adjustl(this%bfr%budtxt))) - case ('FLOW-JA-FACE') - ! - ! -- bfr%flowja contains only reduced connections so there is - ! a one-to-one match with this%gwfflowja - do ipos = 1, size(this%bfr%flowja) - this%gwfflowja(ipos) = this%bfr%flowja(ipos) - end do - case ('DATA-SPDIS') - do i = 1, this%bfr%nlist - nu = this%bfr%nodesrc(i) - nr = this%dis%get_nodenumber(nu, 0) - if (nr <= 0) cycle - this%gwfspdis(1, nr) = this%bfr%auxvar(1, i) - this%gwfspdis(2, nr) = this%bfr%auxvar(2, i) - this%gwfspdis(3, nr) = this%bfr%auxvar(3, i) - end do - case ('DATA-SAT') - do i = 1, this%bfr%nlist - nu = this%bfr%nodesrc(i) - nr = this%dis%get_nodenumber(nu, 0) - if (nr <= 0) cycle - this%gwfsat(nr) = this%bfr%auxvar(1, i) - end do - case ('STO-SS') - do nu = 1, this%dis%nodesuser - nr = this%dis%get_nodenumber(nu, 0) - if (nr <= 0) cycle - this%gwfstrgss(nr) = this%bfr%flow(nu) - end do - case ('STO-SY') - do nu = 1, this%dis%nodesuser - nr = this%dis%get_nodenumber(nu, 0) - if (nr <= 0) cycle - this%gwfstrgsy(nr) = this%bfr%flow(nu) - end do - case default - call this%gwfpackages(ip)%copy_values( & - this%bfr%nlist, & - this%bfr%nodesrc, & - this%bfr%flow, & - this%bfr%auxvar) - do i = 1, this%gwfpackages(ip)%nbound - nu = this%gwfpackages(ip)%nodelist(i) - nr = this%dis%get_nodenumber(nu, 0) - this%gwfpackages(ip)%nodelist(i) = nr - end do - ip = ip + 1 - end select - end do - else - ! - ! -- write message to indicate that flows are being reused - write (this%iout, fmtbudkstpkper) kstp, kper, this%bfr%kstp, this%bfr%kper - ! - ! -- set the flag to indicate that flows were not updated - this%iflowsupdated = 0 - end if - end subroutine advance_bfr - - subroutine finalize_bfr(this) -! ****************************************************************************** -! finalize_bfr -- finalize the budget file reader -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - class(GwtFmiType) :: this - ! -- dummy -! ------------------------------------------------------------------------------ - ! - ! -- Finalize the budget file reader - call this%bfr%finalize() - ! - end subroutine finalize_bfr - - subroutine initialize_hfr(this) -! ****************************************************************************** -! initialize_hfr -- initalize the head file reader -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - class(GwtFmiType) :: 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 - end subroutine initialize_hfr - - subroutine advance_hfr(this) -! ****************************************************************************** -! advance_hfr -- advance the head file reader -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use TdisModule, only: kstp, kper - class(GwtFmiType) :: this - integer(I4B) :: nu, nr, i, ilay - integer(I4B) :: ncpl - real(DP) :: val - logical :: readnext - logical :: success - character(len=*), parameter :: fmtkstpkper = & - &"(1x,/1x,'FMI READING HEAD FOR KSTP ', i0, ' KPER ', i0)" - character(len=*), parameter :: fmthdskstpkper = & - "(1x,/1x, 'FMI SETTING HEAD FOR KSTP ', i0, ' AND KPER ', & - &i0, ' TO BINARY FILE HEADS FROM KSTP ', i0, ' AND KPER ', i0)" -! ------------------------------------------------------------------------------ - ! - ! -- If the latest record read from the head file is from a stress - ! -- period with only one time step, reuse that record (do not read a - ! -- new record) if the GWT model is still in that same stress period, - ! -- or if that record is the last one in the head file. - readnext = .true. - if (kstp * kper > 1) then - if (this%hfr%kstp == 1) then - if (this%hfr%kpernext == kper + 1) then - readnext = .false. - else if (this%hfr%endoffile) then - readnext = .false. - end if - else if (this%hfr%endoffile) then - write (errmsg, '(a)') 'Reached end of GWF head & - &file before reading sufficient head information for this & - &GWT simulation.' - call store_error(errmsg) - call store_error_unit(this%iuhds) - end if - end if - ! - ! -- Read the next record - if (readnext) then - ! - ! -- write to list file that heads are being read - write (this%iout, fmtkstpkper) kstp, kper - ! - ! -- loop through the layered heads for this time step - do ilay = 1, this%hfr%nlay - ! - ! -- read next head chunk - call this%hfr%read_record(success, this%iout) - if (.not. success) then - write (errmsg, '(a)') 'GWF head read not successful' - call store_error(errmsg) - call store_error_unit(this%iuhds) - end if - ! - ! -- Ensure kper is same between model and head file - if (kper /= this%hfr%kper) then - write (errmsg, '(a)') 'Period number in head file & - &does not match period number in transport model. If there & - &is more than one time step in the head file for a given stress & - &period, head file time steps must match gwt model time steps & - &one-for-one in that stress period.' - call store_error(errmsg) - call store_error_unit(this%iuhds) - end if - ! - ! -- if head file kstp > 1, then kstp must match - if (this%hfr%kstp > 1 .and. (kstp /= this%hfr%kstp)) then - write (errmsg, '(a)') 'Time step number in head file & - &does not match time step number in transport model. If there & - &is more than one time step in the head file for a given stress & - &period, head file time steps must match gwt model time steps & - &one-for-one in that stress period.' - call store_error(errmsg) - call store_error_unit(this%iuhds) - end if - ! - ! -- fill the head array for this layer and - ! compress into reduced form - ncpl = size(this%hfr%head) - do i = 1, ncpl - nu = (ilay - 1) * ncpl + i - nr = this%dis%get_nodenumber(nu, 0) - val = this%hfr%head(i) - if (nr > 0) this%gwfhead(nr) = val - end do - end do - else - write (this%iout, fmthdskstpkper) kstp, kper, this%hfr%kstp, this%hfr%kper - end if - end subroutine advance_hfr - - subroutine finalize_hfr(this) -! ****************************************************************************** -! finalize_hfr -- finalize the head file reader -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - class(GwtFmiType) :: this - ! -- dummy -! ------------------------------------------------------------------------------ - ! - ! -- Finalize the head file reader - close (this%iuhds) - ! - end subroutine finalize_hfr - + !> @brief Initialize terms and count unique terms/packages in file subroutine initialize_gwfterms_from_bfr(this) -! ****************************************************************************** -! initialize_gwfterms_from_bfr -- initalize terms and figure out how many -! different terms and packages are contained within the file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use MemoryManagerModule, only: mem_allocate use SimModule, only: store_error, store_error_unit, count_errors @@ -1309,7 +810,6 @@ subroutine initialize_gwfterms_from_bfr(this) logical :: found_stoss logical :: found_stosy integer(I4B), dimension(:), allocatable :: imap -! ------------------------------------------------------------------------------ ! ! -- Calculate the number of gwf flow packages allocate (imap(this%bfr%nbudterms)) @@ -1392,14 +892,8 @@ subroutine initialize_gwfterms_from_bfr(this) return end subroutine initialize_gwfterms_from_bfr + !> @brief Initialize flow terms from a gwf-gwt exchange subroutine initialize_gwfterms_from_gwfbndlist(this) -! ****************************************************************************** -! initialize_gwfterms_from_gwfbndlist -- flows are coming from a gwf-gwt -! exchange -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use BndModule, only: BndType, GetBndFromList ! -- dummy @@ -1413,7 +907,6 @@ subroutine initialize_gwfterms_from_gwfbndlist(this) integer(I4B) :: iterm character(len=LENPACKAGENAME) :: budtxt class(BndType), pointer :: packobj => null() -! ------------------------------------------------------------------------------ ! ! -- determine size of gwf terms ngwfpack = this%gwfbndlist%Count() @@ -1466,15 +959,13 @@ subroutine initialize_gwfterms_from_gwfbndlist(this) return end subroutine initialize_gwfterms_from_gwfbndlist - subroutine allocate_gwfpackages(this, ngwfterms) -! ****************************************************************************** -! allocate_gwfpackages -- gwfpackages is an array of PackageBudget objects. -! This routine allocates gwfpackages to the proper size and initializes some -! member variables. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Allocate GWF packages + !! + !! This routine allocates gwfpackages (an array of PackageBudget + !! objects) to the proper size and initializes member variables. + !! + !< + subroutine gwtfmi_allocate_gwfpackages(this, ngwfterms) ! -- modules use ConstantsModule, only: LENMEMPATH use MemoryManagerModule, only: mem_allocate @@ -1484,7 +975,6 @@ subroutine allocate_gwfpackages(this, ngwfterms) ! -- local integer(I4B) :: n character(len=LENMEMPATH) :: memPath -! ------------------------------------------------------------------------------ ! ! -- direct allocate allocate (this%gwfpackages(ngwfterms)) @@ -1510,21 +1000,15 @@ subroutine allocate_gwfpackages(this, ngwfterms) ! ! -- return return - end subroutine allocate_gwfpackages + end subroutine gwtfmi_allocate_gwfpackages - subroutine deallocate_gwfpackages(this) -! ****************************************************************************** -! deallocate_gwfpackages -- memory in the gwfpackages array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Deallocate memory in the gwfpackages array + subroutine gwtfmi_deallocate_gwfpackages(this) ! -- modules ! -- dummy class(GwtFmiType) :: this ! -- local integer(I4B) :: n -! ------------------------------------------------------------------------------ ! ! -- initialize do n = 1, this%nflowpack @@ -1533,38 +1017,6 @@ subroutine deallocate_gwfpackages(this) ! ! -- return return - end subroutine deallocate_gwfpackages - - subroutine get_package_index(this, name, idx) -! ****************************************************************************** -! get_package_index -- find the package index for package called name -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - use BndModule, only: BndType, GetBndFromList - class(GwtFmiType) :: this - character(len=*), intent(in) :: name - integer(I4B), intent(inout) :: idx - ! -- local - integer(I4B) :: ip -! ------------------------------------------------------------------------------ - ! - ! -- Look through all the packages and return the index with name - idx = 0 - do ip = 1, size(this%flowpacknamearray) - if (this%flowpacknamearray(ip) == name) then - idx = ip - exit - end if - end do - if (idx == 0) then - call store_error('Error in get_package_index. Could not find '//name, & - terminate=.TRUE.) - end if - ! - ! -- return - return - end subroutine get_package_index + end subroutine gwtfmi_deallocate_gwfpackages end module GwtFmiModule diff --git a/src/Model/ModelUtilities/FlowModelInterface.f90 b/src/Model/ModelUtilities/FlowModelInterface.f90 new file mode 100644 index 00000000000..d5c4fd8d0af --- /dev/null +++ b/src/Model/ModelUtilities/FlowModelInterface.f90 @@ -0,0 +1,965 @@ +module FlowModelInterfaceModule + + use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: DONE, DZERO, DHALF, LINELENGTH, LENBUDTXT, & + LENPACKAGENAME + use SimModule, only: store_error, store_error_unit + use SimVariablesModule, only: errmsg + use NumericalPackageModule, only: NumericalPackageType + use BaseDisModule, only: DisBaseType + use ListModule, only: ListType + use BudgetFileReaderModule, only: BudgetFileReaderType + use HeadFileReaderModule, only: HeadFileReaderType + use PackageBudgetModule, only: PackageBudgetType + use BudgetObjectModule, only: BudgetObjectType, budgetobject_cr_bfr + + implicit none + private + public :: FlowModelInterfaceType + + type, extends(NumericalPackageType) :: FlowModelInterfaceType + + character(len=LENPACKAGENAME) :: text = '' !< text string for package + 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 + 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 + real(DP), dimension(:), pointer, contiguous :: gwfsat => null() !< pointer to the GWF saturation array + integer(I4B), dimension(:), pointer, contiguous :: ibdgwfsat0 => null() !< mark cells with saturation = 0 to exclude from dispersion + real(DP), dimension(:), pointer, contiguous :: gwfstrgss => null() !< pointer to flow model QSTOSS + real(DP), dimension(:), pointer, contiguous :: gwfstrgsy => null() !< pointer to flow model QSTOSY + integer(I4B), pointer :: igwfstrgss => null() !< indicates if gwfstrgss is available + integer(I4B), pointer :: igwfstrgsy => null() !< indicates if gwfstrgsy is available + 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 :: 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(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 + character(len=16), dimension(:), allocatable :: flowpacknamearray !< array of boundary package names (e.g. LAK-1, SFR-3, etc.) + contains + + procedure :: advance_bfr + procedure :: advance_hfr + procedure :: allocate_arrays + procedure :: allocate_gwfpackages + procedure :: allocate_scalars + procedure :: deallocate_gwfpackages + procedure :: finalize_bfr + procedure :: finalize_hfr + procedure :: fmi_ar + procedure :: fmi_da + procedure :: fmi_df + procedure :: get_package_index + procedure :: initialize_bfr + procedure :: initialize_gwfterms_from_bfr + procedure :: initialize_gwfterms_from_gwfbndlist + procedure :: initialize_hfr + procedure :: read_options + procedure :: read_packagedata + + end type FlowModelInterfaceType + +contains + + !> @brief Define the flow model interface + subroutine fmi_df(this, dis) + ! -- modules + use SimModule, only: store_error + ! -- dummy + class(FlowModelInterfaceType) :: this + class(DisBaseType), pointer, intent(in) :: dis + ! -- formats + character(len=*), parameter :: fmtfmi = & + "(1x,/1x,'FMI -- FLOW MODEL INTERFACE, VERSION 2, 8/17/2023', & + &' INPUT READ FROM UNIT ', i0, //)" + character(len=*), parameter :: fmtfmi0 = & + "(1x,/1x,'FMI -- FLOW MODEL INTERFACE,'& + &' VERSION 2, 8/17/2023')" + + ! --print a message identifying the FMI package. + if (this%inunit /= 0) then + write (this%iout, fmtfmi) this%inunit + else + write (this%iout, fmtfmi0) + if (this%flows_from_file) then + write (this%iout, '(a)') ' FLOWS ARE ASSUMED TO BE ZERO.' + else + write (this%iout, '(a)') ' FLOWS PROVIDED BY A GWF MODEL IN THIS & + &SIMULATION' + end if + end if + ! + ! -- Store pointers + this%dis => dis + ! + ! -- Read fmi options + if (this%inunit /= 0) then + call this%read_options() + end if + ! + ! -- Read packagedata options + if (this%inunit /= 0 .and. this%flows_from_file) then + call this%read_packagedata() + call this%initialize_gwfterms_from_bfr() + end if + ! + ! -- If GWF-Model exchange is active, setup flow terms + if (.not. this%flows_from_file) then + call this%initialize_gwfterms_from_gwfbndlist() + end if + ! + ! -- Return + return + end subroutine fmi_df + + !> @brief Allocate the package + subroutine fmi_ar(this, ibound) + ! -- modules + use SimModule, only: store_error + ! -- dummy + class(FlowModelInterfaceType) :: this + integer(I4B), dimension(:), pointer, contiguous :: ibound + ! + ! -- store pointers to arguments that were passed in + this%ibound => ibound + ! + ! -- Allocate arrays + call this%allocate_arrays(this%dis%nodes) + ! + ! -- Return + return + end subroutine fmi_ar + + !> @brief Deallocate variables + subroutine fmi_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + ! -- dummy + class(FlowModelInterfaceType) :: this + ! -- todo: finalize hfr and bfr either here or in a finalize routine + ! + ! -- deallocate any memory stored with gwfpackages + call this%deallocate_gwfpackages() + ! + ! -- deallocate fmi arrays + deallocate (this%gwfpackages) + deallocate (this%flowpacknamearray) + call mem_deallocate(this%igwfmvrterm) + call mem_deallocate(this%ibdgwfsat0) + ! + if (this%flows_from_file) then + call mem_deallocate(this%gwfstrgss) + call mem_deallocate(this%gwfstrgsy) + end if + ! + ! -- special treatment, these could be from mem_checkin + call mem_deallocate(this%gwfhead, 'GWFHEAD', this%memoryPath) + call mem_deallocate(this%gwfsat, 'GWFSAT', this%memoryPath) + call mem_deallocate(this%gwfspdis, 'GWFSPDIS', this%memoryPath) + call mem_deallocate(this%gwfflowja, 'GWFFLOWJA', this%memoryPath) + ! + ! -- deallocate scalars + call mem_deallocate(this%flows_from_file) + call mem_deallocate(this%iflowsupdated) + call mem_deallocate(this%igwfstrgss) + call mem_deallocate(this%igwfstrgsy) + call mem_deallocate(this%iubud) + call mem_deallocate(this%iuhds) + call mem_deallocate(this%iumvr) + call mem_deallocate(this%nflowpack) + ! + ! -- deallocate parent + call this%NumericalPackageType%da() + ! + ! -- Return + return + end subroutine fmi_da + + !> @brief Allocate scalars + subroutine allocate_scalars(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate, mem_setptr + ! -- dummy + class(FlowModelInterfaceType) :: this + ! -- local + ! + ! -- allocate scalars in NumericalPackageType + call this%NumericalPackageType%allocate_scalars() + ! + ! -- Allocate + call mem_allocate(this%flows_from_file, 'FLOWS_FROM_FILE', this%memoryPath) + call mem_allocate(this%iflowsupdated, 'IFLOWSUPDATED', this%memoryPath) + call mem_allocate(this%igwfstrgss, 'IGWFSTRGSS', this%memoryPath) + call mem_allocate(this%igwfstrgsy, 'IGWFSTRGSY', this%memoryPath) + 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%nflowpack, 'NFLOWPACK', this%memoryPath) + ! + ! ! + ! -- Initialize + this%flows_from_file = .true. + this%iflowsupdated = 1 + this%igwfstrgss = 0 + this%igwfstrgsy = 0 + this%iubud = 0 + this%iuhds = 0 + this%iumvr = 0 + this%nflowpack = 0 + ! + ! -- Return + return + end subroutine allocate_scalars + + !> @brief Allocate arrays + subroutine allocate_arrays(this, nodes) + use MemoryManagerModule, only: mem_allocate + !modules + use ConstantsModule, only: DZERO + ! -- dummy + class(FlowModelInterfaceType) :: this + integer(I4B), intent(in) :: nodes + ! -- local + integer(I4B) :: n + ! + ! -- Allocate ibdgwfsat0, which is an indicator array marking cells with + ! saturation greater than 0.0 with a value of 1 + call mem_allocate(this%ibdgwfsat0, nodes, 'IBDGWFSAT0', this%memoryPath) + do n = 1, nodes + this%ibdgwfsat0(n) = 1 + end do + ! + ! -- Allocate differently depending on whether or not flows are + ! being read from a file. + if (this%flows_from_file) then + call mem_allocate(this%gwfflowja, this%dis%con%nja, & + 'GWFFLOWJA', this%memoryPath) + call mem_allocate(this%gwfsat, nodes, 'GWFSAT', this%memoryPath) + call mem_allocate(this%gwfhead, nodes, 'GWFHEAD', this%memoryPath) + call mem_allocate(this%gwfspdis, 3, nodes, 'GWFSPDIS', this%memoryPath) + do n = 1, nodes + this%gwfsat(n) = DONE + this%gwfhead(n) = DZERO + this%gwfspdis(:, n) = DZERO + end do + do n = 1, size(this%gwfflowja) + this%gwfflowja(n) = DZERO + end do + ! + ! -- allocate and initialize storage arrays + if (this%igwfstrgss == 0) then + call mem_allocate(this%gwfstrgss, 1, 'GWFSTRGSS', this%memoryPath) + else + call mem_allocate(this%gwfstrgss, nodes, 'GWFSTRGSS', this%memoryPath) + end if + if (this%igwfstrgsy == 0) then + call mem_allocate(this%gwfstrgsy, 1, 'GWFSTRGSY', this%memoryPath) + else + call mem_allocate(this%gwfstrgsy, nodes, 'GWFSTRGSY', this%memoryPath) + end if + do n = 1, size(this%gwfstrgss) + this%gwfstrgss(n) = DZERO + end do + do n = 1, size(this%gwfstrgsy) + this%gwfstrgsy(n) = DZERO + end do + ! + ! -- If there is no fmi package, then there are no flows at all or a + ! connected GWF model, so allocate gwfpackages to zero + if (this%inunit == 0) call this%allocate_gwfpackages(this%nflowpack) + end if + ! + ! -- Return + return + end subroutine allocate_arrays + + !> @brief Read options from input file + subroutine read_options(this) + ! -- modules + use ConstantsModule, only: LINELENGTH, DEM6 + use InputOutputModule, only: getunit, openfile, urdaux + use SimModule, only: store_error, store_error_unit + ! -- dummy + class(FlowModelInterfaceType) :: this + ! -- local + character(len=LINELENGTH) :: keyword + integer(I4B) :: ierr + logical :: isfound, endOfBlock + ! + ! -- get options block + call this%parser%GetBlock('OPTIONS', isfound, ierr, blockRequired=.false., & + supportOpenClose=.true.) + ! + ! -- parse options block if detected + if (isfound) then + write (this%iout, '(1x,a)') 'PROCESSING FMI OPTIONS' + do + call this%parser%GetNextLine(endOfBlock) + if (endOfBlock) exit + call this%parser%GetStringCaps(keyword) + select case (keyword) + case ('SAVE_FLOWS') + this%ipakcb = -1 + case default + write (errmsg, '(a,3(1x,a))') & + 'UNKNOWN', trim(adjustl(this%text)), 'OPTION:', trim(keyword) + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end select + end do + write (this%iout, '(1x,a)') 'END OF FMI OPTIONS' + end if + ! + ! -- return + return + end subroutine read_options + + !> @brief Read packagedata block from input file + subroutine read_packagedata(this) + ! -- modules + use OpenSpecModule, only: ACCESS, FORM + use ConstantsModule, only: LINELENGTH, DEM6, LENPACKAGENAME + use InputOutputModule, only: getunit, openfile, urdaux + use SimModule, only: store_error, store_error_unit + ! -- dummy + class(FlowModelInterfaceType) :: this + ! -- local + character(len=LINELENGTH) :: keyword, fname + integer(I4B) :: ierr + integer(I4B) :: inunit + integer(I4B) :: iapt + logical :: isfound, endOfBlock + logical :: blockrequired + logical :: exist + ! + ! -- initialize + iapt = 0 + blockrequired = .true. + ! + ! -- get options block + call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, & + blockRequired=blockRequired, & + supportOpenClose=.true.) + ! + ! -- parse options block if detected + if (isfound) then + write (this%iout, '(1x,a)') 'PROCESSING FMI PACKAGEDATA' + do + call this%parser%GetNextLine(endOfBlock) + if (endOfBlock) exit + call this%parser%GetStringCaps(keyword) + select case (keyword) + case ('GWFBUDGET') + call this%parser%GetStringCaps(keyword) + if (keyword /= 'FILEIN') then + call store_error('GWFBUDGET KEYWORD MUST BE FOLLOWED BY '// & + '"FILEIN" then by filename.') + call this%parser%StoreErrorUnit() + end if + call this%parser%GetString(fname) + inunit = getunit() + inquire (file=trim(fname), exist=exist) + if (.not. exist) then + call store_error('Could not find file '//trim(fname)) + call this%parser%StoreErrorUnit() + end if + call openfile(inunit, this%iout, fname, 'DATA(BINARY)', FORM, & + ACCESS, 'UNKNOWN') + this%iubud = inunit + call this%initialize_bfr() + case ('GWFHEAD') + call this%parser%GetStringCaps(keyword) + if (keyword /= 'FILEIN') then + call store_error('GWFHEAD KEYWORD MUST BE FOLLOWED BY '// & + '"FILEIN" then by filename.') + call this%parser%StoreErrorUnit() + end if + call this%parser%GetString(fname) + inquire (file=trim(fname), exist=exist) + if (.not. exist) then + call store_error('Could not find file '//trim(fname)) + call this%parser%StoreErrorUnit() + end if + inunit = getunit() + call openfile(inunit, this%iout, fname, 'DATA(BINARY)', FORM, & + ACCESS, 'UNKNOWN') + this%iuhds = inunit + call this%initialize_hfr() + case ('GWFMOVER') + call this%parser%GetStringCaps(keyword) + if (keyword /= 'FILEIN') then + call store_error('GWFMOVER 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%iumvr = inunit + call budgetobject_cr_bfr(this%mvrbudobj, 'MVT', this%iumvr, & ! kluge note: MVT? + this%iout) + call this%mvrbudobj%fill_from_bfr(this%dis, this%iout) + case default + write (errmsg, '(a,3(1x,a))') & + 'UNKNOWN', trim(adjustl(this%text)), 'PACKAGEDATA:', trim(keyword) + call store_error(errmsg) + end select + end do + write (this%iout, '(1x,a)') 'END OF FMI PACKAGEDATA' + end if + ! + ! -- return + return + 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 + 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 + use TdisModule, only: kstp, kper + ! -- dummy + class(FlowModelInterfaceType) :: this + ! -- local + logical :: success + integer(I4B) :: n + integer(I4B) :: ipos + integer(I4B) :: nu, nr + integer(I4B) :: ip, i + logical :: readnext + ! -- format + character(len=*), parameter :: fmtkstpkper = & + "(1x,/1x,'FMI READING BUDGET TERMS & + &FOR KSTP ', i0, ' KPER ', i0)" + character(len=*), parameter :: fmtbudkstpkper = & + "(1x,/1x, 'FMI SETTING BUDGET TERMS & + &FOR KSTP ', i0, ' AND KPER ', & + &i0, ' TO BUDGET FILE TERMS FROM & + &KSTP ', i0, ' AND KPER ', i0)" + ! + ! -- If the latest record read from the budget file is from a stress + ! -- period with only one time step, reuse that record (do not read a + ! -- new record) if the running model is still in that same stress period, + ! -- or if that record is the last one in the budget file. + readnext = .true. + if (kstp * kper > 1) then + if (this%bfr%kstp == 1) then + if (this%bfr%kpernext == kper + 1) then + readnext = .false. + else if (this%bfr%endoffile) then + readnext = .false. + end if + else if (this%bfr%endoffile) then + write (errmsg, '(4x,a)') 'REACHED END OF GWF BUDGET & + &FILE BEFORE READING SUFFICIENT BUDGET INFORMATION FOR THIS & + &GWT SIMULATION.' + call store_error(errmsg) + call store_error_unit(this%iubud) + end if + end if + ! + ! -- Read the next record + if (readnext) then + ! + ! -- Write the current time step and stress period + write (this%iout, fmtkstpkper) kstp, kper + ! + ! -- loop through the budget terms for this stress period + ! i is the counter for gwf flow packages + ip = 1 + do n = 1, this%bfr%nbudterms + call this%bfr%read_record(success, this%iout) + if (.not. success) then + write (errmsg, '(4x,a)') 'GWF BUDGET READ NOT SUCCESSFUL' + call store_error(errmsg) + call store_error_unit(this%iubud) + end if + ! + ! -- Ensure kper is same between model and budget file + if (kper /= this%bfr%kper) then + write (errmsg, '(4x,a)') 'PERIOD NUMBER IN BUDGET FILE & + &DOES NOT MATCH PERIOD NUMBER IN TRANSPORT MODEL. IF THERE & + &IS MORE THAN ONE TIME STEP IN THE BUDGET FILE FOR A GIVEN STRESS & + &PERIOD, BUDGET FILE TIME STEPS MUST MATCH GWT MODEL TIME STEPS & + &ONE-FOR-ONE IN THAT STRESS PERIOD.' + call store_error(errmsg) + call store_error_unit(this%iubud) + end if + ! + ! -- if budget file kstp > 1, then kstp must match + if (this%bfr%kstp > 1 .and. (kstp /= this%bfr%kstp)) then + write (errmsg, '(4x,a)') 'TIME STEP NUMBER IN BUDGET FILE & + &DOES NOT MATCH TIME STEP NUMBER IN TRANSPORT MODEL. IF THERE & + &IS MORE THAN ONE TIME STEP IN THE BUDGET FILE FOR A GIVEN STRESS & + &PERIOD, BUDGET FILE TIME STEPS MUST MATCH GWT MODEL TIME STEPS & + &ONE-FOR-ONE IN THAT STRESS PERIOD.' + call store_error(errmsg) + call store_error_unit(this%iubud) + end if + ! + ! -- parse based on the type of data, and compress all user node + ! numbers into reduced node numbers + select case (trim(adjustl(this%bfr%budtxt))) + case ('FLOW-JA-FACE') + ! + ! -- bfr%flowja contains only reduced connections so there is + ! a one-to-one match with this%gwfflowja + do ipos = 1, size(this%bfr%flowja) + this%gwfflowja(ipos) = this%bfr%flowja(ipos) + end do + case ('DATA-SPDIS') + do i = 1, this%bfr%nlist + nu = this%bfr%nodesrc(i) + nr = this%dis%get_nodenumber(nu, 0) + if (nr <= 0) cycle + this%gwfspdis(1, nr) = this%bfr%auxvar(1, i) + this%gwfspdis(2, nr) = this%bfr%auxvar(2, i) + this%gwfspdis(3, nr) = this%bfr%auxvar(3, i) + end do + case ('DATA-SAT') + do i = 1, this%bfr%nlist + nu = this%bfr%nodesrc(i) + nr = this%dis%get_nodenumber(nu, 0) + if (nr <= 0) cycle + this%gwfsat(nr) = this%bfr%auxvar(1, i) + end do + case ('STO-SS') + do nu = 1, this%dis%nodesuser + nr = this%dis%get_nodenumber(nu, 0) + if (nr <= 0) cycle + this%gwfstrgss(nr) = this%bfr%flow(nu) + end do + case ('STO-SY') + do nu = 1, this%dis%nodesuser + nr = this%dis%get_nodenumber(nu, 0) + if (nr <= 0) cycle + this%gwfstrgsy(nr) = this%bfr%flow(nu) + end do + case default + call this%gwfpackages(ip)%copy_values( & + this%bfr%nlist, & + this%bfr%nodesrc, & + this%bfr%flow, & + this%bfr%auxvar) + do i = 1, this%gwfpackages(ip)%nbound + nu = this%gwfpackages(ip)%nodelist(i) + nr = this%dis%get_nodenumber(nu, 0) + this%gwfpackages(ip)%nodelist(i) = nr + end do + ip = ip + 1 + end select + end do + else + ! + ! -- write message to indicate that flows are being reused + write (this%iout, fmtbudkstpkper) kstp, kper, this%bfr%kstp, this%bfr%kper + ! + ! -- set the flag to indicate that flows were not updated + this%iflowsupdated = 0 + end if + 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 + end subroutine initialize_hfr + + !> @brief Advance the head file reader + subroutine advance_hfr(this) + ! -- modules + use TdisModule, only: kstp, kper + class(FlowModelInterfaceType) :: this + integer(I4B) :: nu, nr, i, ilay + integer(I4B) :: ncpl + real(DP) :: val + logical :: readnext + logical :: success + character(len=*), parameter :: fmtkstpkper = & + "(1x,/1x,'FMI READING HEAD FOR & + &KSTP ', i0, ' KPER ', i0)" + character(len=*), parameter :: fmthdskstpkper = & + "(1x,/1x, 'FMI SETTING HEAD FOR KSTP ', i0, ' AND KPER ', & + &i0, ' TO BINARY FILE HEADS FROM KSTP ', i0, ' AND KPER ', i0)" + ! + ! -- If the latest record read from the head file is from a stress + ! -- period with only one time step, reuse that record (do not read a + ! -- new record) if the running model is still in that same stress period, + ! -- or if that record is the last one in the head file. + readnext = .true. + if (kstp * kper > 1) then + if (this%hfr%kstp == 1) then + if (this%hfr%kpernext == kper + 1) then + readnext = .false. + else if (this%hfr%endoffile) then + readnext = .false. + end if + else if (this%hfr%endoffile) then + write (errmsg, '(4x,a)') 'REACHED END OF GWF HEAD & + &FILE BEFORE READING SUFFICIENT HEAD INFORMATION FOR THIS & + &GWT SIMULATION.' + call store_error(errmsg) + call store_error_unit(this%iuhds) + end if + end if + ! + ! -- Read the next record + if (readnext) then + ! + ! -- write to list file that heads are being read + write (this%iout, fmtkstpkper) kstp, kper + ! + ! -- loop through the layered heads for this time step + do ilay = 1, this%hfr%nlay + ! + ! -- read next head chunk + call this%hfr%read_record(success, this%iout) + if (.not. success) then + write (errmsg, '(4x,a)') 'GWF HEAD READ NOT SUCCESSFUL' + call store_error(errmsg) + call store_error_unit(this%iuhds) + end if + ! + ! -- Ensure kper is same between model and head file + if (kper /= this%hfr%kper) then + write (errmsg, '(4x,a)') 'PERIOD NUMBER IN HEAD FILE & + &DOES NOT MATCH PERIOD NUMBER IN TRANSPORT MODEL. IF THERE & + &IS MORE THAN ONE TIME STEP IN THE HEAD FILE FOR A GIVEN STRESS & + &PERIOD, HEAD FILE TIME STEPS MUST MATCH GWT MODEL TIME STEPS & + &ONE-FOR-ONE IN THAT STRESS PERIOD.' + call store_error(errmsg) + call store_error_unit(this%iuhds) + end if + ! + ! -- if head file kstp > 1, then kstp must match + if (this%hfr%kstp > 1 .and. (kstp /= this%hfr%kstp)) then + write (errmsg, '(4x,a)') 'TIME STEP NUMBER IN HEAD FILE & + &DOES NOT MATCH TIME STEP NUMBER IN TRANSPORT MODEL. IF THERE & + &IS MORE THAN ONE TIME STEP IN THE HEAD FILE FOR A GIVEN STRESS & + &PERIOD, HEAD FILE TIME STEPS MUST MATCH GWT MODEL TIME STEPS & + &ONE-FOR-ONE IN THAT STRESS PERIOD.' + call store_error(errmsg) + call store_error_unit(this%iuhds) + end if + ! + ! -- fill the head array for this layer and + ! compress into reduced form + ncpl = size(this%hfr%head) + do i = 1, ncpl + nu = (ilay - 1) * ncpl + i + nr = this%dis%get_nodenumber(nu, 0) + val = this%hfr%head(i) + if (nr > 0) this%gwfhead(nr) = val + end do + end do + else + write (this%iout, fmthdskstpkper) kstp, kper, this%hfr%kstp, this%hfr%kper + end if + 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 + !! + !! initalize terms and figure out how many + !! different terms and packages are contained within the file + !! + !< + subroutine initialize_gwfterms_from_bfr(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + use SimModule, only: store_error, store_error_unit, count_errors + ! -- dummy + class(FlowModelInterfaceType) :: this + ! -- local + integer(I4B) :: nflowpack + integer(I4B) :: i, ip + integer(I4B) :: naux + logical :: found_flowja + logical :: found_dataspdis + logical :: found_datasat + logical :: found_stoss + logical :: found_stosy + integer(I4B), dimension(:), allocatable :: imap + ! + ! -- Calculate the number of gwf flow packages + allocate (imap(this%bfr%nbudterms)) + imap(:) = 0 + nflowpack = 0 + found_flowja = .false. + found_dataspdis = .false. + found_datasat = .false. + found_stoss = .false. + found_stosy = .false. + do i = 1, this%bfr%nbudterms + select case (trim(adjustl(this%bfr%budtxtarray(i)))) + case ('FLOW-JA-FACE') + found_flowja = .true. + case ('DATA-SPDIS') + found_dataspdis = .true. + case ('DATA-SAT') + found_datasat = .true. + case ('STO-SS') + found_stoss = .true. + this%igwfstrgss = 1 + case ('STO-SY') + found_stosy = .true. + this%igwfstrgsy = 1 + case default + nflowpack = nflowpack + 1 + imap(i) = 1 + end select + end do + ! + ! -- allocate gwfpackage arrays + call this%allocate_gwfpackages(nflowpack) + ! + ! -- Copy the package name and aux names from budget file reader + ! to the gwfpackages derived-type variable + ip = 1 + do i = 1, this%bfr%nbudterms + if (imap(i) == 0) cycle + call this%gwfpackages(ip)%set_name(this%bfr%dstpackagenamearray(i), & + this%bfr%budtxtarray(i)) + naux = this%bfr%nauxarray(i) + call this%gwfpackages(ip)%set_auxname(naux, this%bfr%auxtxtarray(1:naux, i)) + ip = ip + 1 + end do + ! + ! -- Copy just the package names for the boundary packages into + ! the flowpacknamearray + ip = 1 + do i = 1, size(imap) + if (imap(i) == 1) then + this%flowpacknamearray(ip) = this%bfr%dstpackagenamearray(i) + ip = ip + 1 + end if + end do + ! + ! -- Error if specific discharge, saturation or flowja not found + if (.not. found_dataspdis) then + write (errmsg, '(4x,a)') 'SPECIFIC DISCHARGE NOT FOUND IN & + &BUDGET FILE. SAVE_SPECIFIC_DISCHARGE AND & + &SAVE_FLOWS MUST BE ACTIVATED IN THE NPF PACKAGE.' + call store_error(errmsg) + end if + if (.not. found_datasat) then + write (errmsg, '(4x,a)') 'SATURATION NOT FOUND IN & + &BUDGET FILE. SAVE_SATURATION AND & + &SAVE_FLOWS MUST BE ACTIVATED IN THE NPF PACKAGE.' + call store_error(errmsg) + end if + if (.not. found_flowja) then + write (errmsg, '(4x,a)') 'FLOWJA NOT FOUND IN & + &BUDGET FILE. SAVE_FLOWS MUST & + &BE ACTIVATED IN THE NPF PACKAGE.' + call store_error(errmsg) + end if + if (count_errors() > 0) then + call this%parser%StoreErrorUnit() + end if + ! + ! -- return + return + 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 + ! -- dummy + class(FlowModelInterfaceType) :: this + ! -- local + integer(I4B) :: ngwfpack + integer(I4B) :: ngwfterms + integer(I4B) :: ip + integer(I4B) :: imover + integer(I4B) :: ntomvr + integer(I4B) :: iterm + character(len=LENPACKAGENAME) :: budtxt + class(BndType), pointer :: packobj => null() + ! + ! -- determine size of gwf terms + ngwfpack = this%gwfbndlist%Count() + ! + ! -- Count number of to-mvr terms, but do not include advanced packages + ! as those mover terms are not losses from the cell, but rather flows + ! within the advanced package + ntomvr = 0 + do ip = 1, ngwfpack + packobj => GetBndFromList(this%gwfbndlist, ip) + imover = packobj%imover + if (packobj%isadvpak /= 0) imover = 0 + if (imover /= 0) then + ntomvr = ntomvr + 1 + end if + end do + ! + ! -- Allocate arrays in fmi of size ngwfterms, which is the number of + ! packages plus the number of packages with mover terms. + ngwfterms = ngwfpack + ntomvr + call this%allocate_gwfpackages(ngwfterms) + ! + ! -- Assign values in the fmi package + iterm = 1 + do ip = 1, ngwfpack + ! + ! -- set and store names + packobj => GetBndFromList(this%gwfbndlist, ip) + budtxt = adjustl(packobj%text) + call this%gwfpackages(iterm)%set_name(packobj%packName, budtxt) + this%flowpacknamearray(iterm) = packobj%packName + iterm = iterm + 1 + ! + ! -- if this package has a mover associated with it, then add another + ! term that corresponds to the mover flows + imover = packobj%imover + if (packobj%isadvpak /= 0) imover = 0 + if (imover /= 0) then + budtxt = trim(adjustl(packobj%text))//'-TO-MVR' + call this%gwfpackages(iterm)%set_name(packobj%packName, budtxt) + this%flowpacknamearray(iterm) = packobj%packName + this%igwfmvrterm(iterm) = 1 + iterm = iterm + 1 + end if + end do + return + end subroutine initialize_gwfterms_from_gwfbndlist + + !> @brief Allocate budget packages + !! + !! gwfpackages is an array of PackageBudget objects. + !! This routine allocates gwfpackages to the proper size and initializes some + !! member variables. + !! + !< + subroutine allocate_gwfpackages(this, ngwfterms) + ! -- modules + use ConstantsModule, only: LENMEMPATH + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(FlowModelInterfaceType) :: this + integer(I4B), intent(in) :: ngwfterms + ! -- local + integer(I4B) :: n + character(len=LENMEMPATH) :: memPath + ! + ! -- direct allocate + allocate (this%gwfpackages(ngwfterms)) + allocate (this%flowpacknamearray(ngwfterms)) + ! + ! -- mem_allocate + call mem_allocate(this%igwfmvrterm, ngwfterms, 'IGWFMVRTERM', this%memoryPath) + ! + ! -- initialize + this%nflowpack = ngwfterms + do n = 1, this%nflowpack + this%igwfmvrterm(n) = 0 + this%flowpacknamearray(n) = '' + ! + ! -- Create a mempath for each individual flow package data set + ! of the form, MODELNAME/FMI-FTn + write (memPath, '(a, i0)') trim(this%memoryPath)//'-FT', n + call this%gwfpackages(n)%initialize(memPath) + end do + ! + ! -- return + return + 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 + ! + ! -- return + return + 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 + character(len=*), intent(in) :: name + integer(I4B), intent(inout) :: idx + ! -- local + integer(I4B) :: ip + ! + ! -- Look through all the packages and return the index with name + idx = 0 + do ip = 1, size(this%flowpacknamearray) + if (this%flowpacknamearray(ip) == name) then + idx = ip + exit + end if + end do + if (idx == 0) then + call store_error('Error in get_package_index. Could not find '//name, & + terminate=.TRUE.) + end if + ! + ! -- return + return + end subroutine get_package_index + +end module FlowModelInterfaceModule diff --git a/src/meson.build b/src/meson.build index 57ce656dd37..08417b774e9 100644 --- a/src/meson.build +++ b/src/meson.build @@ -107,6 +107,7 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'Connections.f90', 'Model' / 'ModelUtilities' / 'DiscretizationBase.f90', 'Model' / 'ModelUtilities' / 'DisvGeom.f90', + 'Model' / 'ModelUtilities' / 'FlowModelInterface.f90', 'Model' / 'ModelUtilities' / 'GwfBuyInputData.f90', 'Model' / 'ModelUtilities' / 'GwfMvrPeriodData.f90', 'Model' / 'ModelUtilities' / 'GwfNpfOptions.f90',