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',