From f7b6f76f3f8b4770f0dad5486f91a5c6a98b5777 Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Fri, 13 Oct 2023 15:26:06 -0700 Subject: [PATCH] refactor(tsp): Elevate OBS to generalized transport class (#1400) --- make/makefile | 6 +- msvs/mf6core.vfproj | 2 +- src/Model/Connection/GwtInterfaceModel.f90 | 4 +- src/Model/GroundWaterTransport/gwt1.f90 | 242 +--------------- src/Model/TransportModel/tsp1.f90 | 261 +++++++++++++++++- .../tsp1obs1.f90} | 173 +++++------- src/meson.build | 2 +- 7 files changed, 351 insertions(+), 339 deletions(-) rename src/Model/{GroundWaterTransport/gwt1obs1.f90 => TransportModel/tsp1obs1.f90} (56%) diff --git a/make/makefile b/make/makefile index 137466623fd..6037f3af64d 100644 --- a/make/makefile +++ b/make/makefile @@ -194,6 +194,7 @@ $(OBJDIR)/SimStages.o \ $(OBJDIR)/NumericalModel.o \ $(OBJDIR)/FlowModelInterface.o \ $(OBJDIR)/OutputControlData.o \ +$(OBJDIR)/gwf3ic8.o \ $(OBJDIR)/Xt3dAlgorithm.o \ $(OBJDIR)/gwf3tvbase8.o \ $(OBJDIR)/gwf3sfr8.o \ @@ -210,7 +211,7 @@ $(OBJDIR)/BaseExchange.o \ $(OBJDIR)/tsp1fmi1.o \ $(OBJDIR)/GwtSpc.o \ $(OBJDIR)/OutputControl.o \ -$(OBJDIR)/gwf3ic8.o \ +$(OBJDIR)/tsp1ic1.o \ $(OBJDIR)/TspAdvOptions.o \ $(OBJDIR)/UzfCellGroup.o \ $(OBJDIR)/Xt3dInterface.o \ @@ -225,8 +226,8 @@ $(OBJDIR)/CellWithNbrs.o \ $(OBJDIR)/NumericalExchange.o \ $(OBJDIR)/tsp1ssm1.o \ $(OBJDIR)/tsp1oc1.o \ +$(OBJDIR)/tsp1obs1.o \ $(OBJDIR)/tsp1mvt1.o \ -$(OBJDIR)/tsp1ic1.o \ $(OBJDIR)/tsp1adv1.o \ $(OBJDIR)/gwf3disv8.o \ $(OBJDIR)/gwf3disu8.o \ @@ -254,7 +255,6 @@ $(OBJDIR)/tsp1.o \ $(OBJDIR)/gwt1uzt1.o \ $(OBJDIR)/gwt1src1.o \ $(OBJDIR)/gwt1sft1.o \ -$(OBJDIR)/gwt1obs1.o \ $(OBJDIR)/gwt1mwt1.o \ $(OBJDIR)/gwt1lkt1.o \ $(OBJDIR)/gwt1ist1.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index a70e28355f0..9c91d969729 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -175,7 +175,6 @@ - @@ -207,6 +206,7 @@ + diff --git a/src/Model/Connection/GwtInterfaceModel.f90 b/src/Model/Connection/GwtInterfaceModel.f90 index 435410a9712..4adfda4c36d 100644 --- a/src/Model/Connection/GwtInterfaceModel.f90 +++ b/src/Model/Connection/GwtInterfaceModel.f90 @@ -12,7 +12,7 @@ module GwtInterfaceModelModule use GwtDspModule, only: dsp_cr, GwtDspType use GwtDspOptionsModule, only: GwtDspOptionsType use GwtMstModule, only: mst_cr - use GwtObsModule, only: gwt_obs_cr + use TspObsModule, only: tsp_obs_cr use GridConnectionModule implicit none @@ -87,7 +87,7 @@ subroutine gwtifmod_cr(this, name, iout, gridConn) call adv_cr(this%adv, this%name, adv_unit, this%iout, this%fmi, & this%ieqnsclfac) call dsp_cr(this%dsp, this%name, '', -dsp_unit, this%iout, this%fmi) - call gwt_obs_cr(this%obs, inobs) + call tsp_obs_cr(this%obs, inobs) end subroutine gwtifmod_cr diff --git a/src/Model/GroundWaterTransport/gwt1.f90 b/src/Model/GroundWaterTransport/gwt1.f90 index b1cc864d250..3d12b5db786 100644 --- a/src/Model/GroundWaterTransport/gwt1.f90 +++ b/src/Model/GroundWaterTransport/gwt1.f90 @@ -12,13 +12,13 @@ module GwtModule LENPAKLOC, LENVARNAME use VersionModule, only: write_listfile_header use NumericalModelModule, only: NumericalModelType - use TransportModelModule, only: TransportModelType + use BaseModelModule, only: BaseModelType use BndModule, only: BndType, AddBndToList, GetBndFromList use GwtDspModule, only: GwtDspType use GwtMstModule, only: GwtMstType - use GwtObsModule, only: GwtObsType use BudgetModule, only: BudgetType + use TransportModelModule use MatrixBaseModule implicit none @@ -35,10 +35,8 @@ module GwtModule type(GwtMstType), pointer :: mst => null() ! mass storage and transfer package type(GwtDspType), pointer :: dsp => null() ! dispersion package - type(GwtObsType), pointer :: obs => null() ! observation package integer(I4B), pointer :: inmst => null() ! unit number MST integer(I4B), pointer :: indsp => null() ! DSP enabled flag - integer(I4B), pointer :: inobs => null() ! unit number OBS contains @@ -58,11 +56,6 @@ module GwtModule procedure :: model_bdentry => gwt_bdentry procedure :: allocate_scalars procedure :: get_iasym => gwt_get_iasym - procedure, private :: gwt_ot_flow - procedure, private :: gwt_ot_flowja - procedure, private :: gwt_ot_dv - procedure, private :: gwt_ot_bdsummary - procedure, private :: gwt_ot_obs procedure :: create_packages => create_gwt_packages procedure, private :: create_bndpkgs procedure, private :: package_create @@ -77,7 +70,7 @@ subroutine gwt_cr(filename, id, modelname) ! -- modules use ListsModule, only: basemodellist use BaseModelModule, only: AddBaseModelToList - use ConstantsModule, only: LINELENGTH + use ConstantsModule, only: LINELENGTH, LENPACKAGENAME use MemoryHelperModule, only: create_mem_path use MemoryManagerExtModule, only: mem_set_value use SimVariablesModule, only: idm_context @@ -281,7 +274,7 @@ subroutine gwt_ar(this) if (this%inadv > 0) call this%adv%adv_ar(this%dis, this%ibound) if (this%indsp > 0) call this%dsp%dsp_ar(this%ibound, this%mst%thetam) if (this%inssm > 0) call this%ssm%ssm_ar(this%dis, this%ibound, this%x) - if (this%inobs > 0) call this%obs%gwt_obs_ar(this%ic, this%x, this%flowja) + if (this%inobs > 0) call this%obs%tsp_obs_ar(this%ic, this%x, this%flowja) ! ! -- Set governing equation scale factor. Note that this scale factor ! -- cannot be set arbitrarily. For solute transport, it must be set @@ -595,229 +588,27 @@ end subroutine gwt_bd !! Call the parent class output routine !< subroutine gwt_ot(this) - ! -- modules - use TdisModule, only: kstp, kper, tdis_ot, endofperiod ! -- dummy class(GwtModelType) :: this ! -- local - integer(I4B) :: idvsave - integer(I4B) :: idvprint integer(I4B) :: icbcfl integer(I4B) :: icbcun - integer(I4B) :: ibudfl - integer(I4B) :: ipflag - ! -- formats - 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 - idvprint = 0 + ! + ! -- Initialize icbcfl = 0 - ibudfl = 0 - if (this%oc%oc_save('CONCENTRATION')) idvsave = 1 - if (this%oc%oc_print('CONCENTRATION')) idvprint = 1 + ! + ! -- Because mst belongs to gwt, call mst_ot_flow directly (and not from parent) if (this%oc%oc_save('BUDGET')) icbcfl = 1 - if (this%oc%oc_print('BUDGET')) ibudfl = 1 icbcun = this%oc%oc_save_unit('BUDGET') - ! - ! -- Override ibudfl and idvprint flags for nonconvergence - ! and end of period - ibudfl = this%oc%set_print_flag('BUDGET', this%icnvg, endofperiod) - idvprint = this%oc%set_print_flag('CONCENTRATION', this%icnvg, endofperiod) - ! - ! Calculate and save observations - call this%gwt_ot_obs() - ! - ! Save and print flows - call this%gwt_ot_flow(icbcfl, ibudfl, icbcun) - ! - ! Save and print dependent variables - call this%gwt_ot_dv(idvsave, idvprint, ipflag) - ! - ! Print budget summaries - call this%gwt_ot_bdsummary(ibudfl, ipflag) - ! - ! -- Timing Output; if any dependendent variables or budgets - ! are printed, then ipflag is set to 1. - if (ipflag == 1) call tdis_ot(this%iout) - ! - ! -- Write non-convergence message - if (this%icnvg == 0) then - write (this%iout, fmtnocnvg) kstp, kper - end if - ! - ! -- Return - return - end subroutine gwt_ot - - !> @brief Calculate and save observations - !< - subroutine gwt_ot_obs(this) - class(GwtModelType) :: this - class(BndType), pointer :: packobj - integer(I4B) :: ip - - ! -- Calculate and save observations - call this%obs%obs_bd() - call this%obs%obs_ot() - - ! -- Calculate and save package obserations - do ip = 1, this%bndlist%Count() - packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_bd_obs() - call packobj%bnd_ot_obs() - end do - ! - ! -- Return - return - end subroutine gwt_ot_obs - - !> @brief Save flows - !< - subroutine gwt_ot_flow(this, icbcfl, ibudfl, icbcun) - class(GwtModelType) :: this - integer(I4B), intent(in) :: icbcfl - integer(I4B), intent(in) :: ibudfl - integer(I4B), intent(in) :: icbcun - class(BndType), pointer :: packobj - integer(I4B) :: ip - - ! -- Save GWT flows - call this%gwt_ot_flowja(this%nja, this%flowja, icbcfl, icbcun) if (this%inmst > 0) call this%mst%mst_ot_flow(icbcfl, icbcun) - if (this%infmi > 0) call this%fmi%fmi_ot_flow(icbcfl, icbcun) - if (this%inssm > 0) then - call this%ssm%ssm_ot_flow(icbcfl=icbcfl, ibudfl=0, icbcun=icbcun) - end if - do ip = 1, this%bndlist%Count() - packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=0, icbcun=icbcun) - end do - - ! -- Save advanced package flows - do ip = 1, this%bndlist%Count() - packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_ot_package_flows(icbcfl=icbcfl, ibudfl=0) - end do - if (this%inmvt > 0) then - call this%mvt%mvt_ot_saveflow(icbcfl, ibudfl) - end if - ! - ! -- Print GWF flows - ! no need to print flowja - ! no need to print mst - ! no need to print fmi - if (this%inssm > 0) then - call this%ssm%ssm_ot_flow(icbcfl=icbcfl, ibudfl=ibudfl, icbcun=0) - end if - do ip = 1, this%bndlist%Count() - packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=ibudfl, icbcun=0) - end do ! - ! -- Print advanced package flows - do ip = 1, this%bndlist%Count() - packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_ot_package_flows(icbcfl=0, ibudfl=ibudfl) - end do - if (this%inmvt > 0) then - call this%mvt%mvt_ot_printflow(icbcfl, ibudfl) - end if + ! -- Call parent class _ot routines. + call this%tsp_ot(this%inmst) ! ! -- Return return - end subroutine gwt_ot_flow - - !> @brief Write intercell flows - !< - subroutine gwt_ot_flowja(this, nja, flowja, icbcfl, icbcun) - ! -- dummy - class(GwtModelType) :: this - integer(I4B), intent(in) :: nja - real(DP), dimension(nja), intent(in) :: flowja - integer(I4B), intent(in) :: icbcfl - integer(I4B), intent(in) :: icbcun - ! -- local - integer(I4B) :: ibinun - ! -- formats - ! - ! -- Set unit number for binary output - if (this%ipakcb < 0) then - ibinun = icbcun - elseif (this%ipakcb == 0) then - ibinun = 0 - else - ibinun = this%ipakcb - end if - if (icbcfl == 0) ibinun = 0 - ! - ! -- Write the face flows if requested - if (ibinun /= 0) then - call this%dis%record_connection_array(flowja, ibinun, this%iout) - end if - ! - ! -- Return - 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 - integer(I4B), intent(in) :: idvprint - integer(I4B), intent(inout) :: ipflag - class(BndType), pointer :: packobj - integer(I4B) :: ip - ! - ! -- Print advanced package dependent variables - do ip = 1, this%bndlist%Count() - packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_ot_dv(idvsave, idvprint) - end do - ! - ! -- save head and print head - call this%oc%oc_ot(ipflag) - ! - ! -- Return - return - 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 - integer(I4B), intent(in) :: ibudfl - integer(I4B), intent(inout) :: ipflag - class(BndType), pointer :: packobj - integer(I4B) :: ip - ! - ! -- Package budget summary - do ip = 1, this%bndlist%Count() - packobj => GetBndFromList(this%bndlist, ip) - call packobj%bnd_ot_bdsummary(kstp, kper, this%iout, ibudfl) - end do - ! - ! -- mover budget summary - if (this%inmvt > 0) then - call this%mvt%mvt_ot_bdsummary(ibudfl) - end if - ! - ! -- model budget summary - if (ibudfl /= 0) then - ipflag = 1 - call this%budget%budget_ot(kstp, kper, this%iout) - end if - ! - ! -- Write to budget csv - call this%budget%writecsv(totim) - ! - ! -- Return - return - end subroutine gwt_ot_bdsummary + end subroutine gwt_ot !> @brief Deallocate !! @@ -873,7 +664,6 @@ subroutine gwt_da(this) ! -- Scalars call mem_deallocate(this%indsp) call mem_deallocate(this%inmst) - call mem_deallocate(this%inobs) ! ! -- Parent class members call this%TransportModelType%tsp_da() @@ -956,14 +746,12 @@ subroutine allocate_scalars(this, modelname) ! -- allocate parent class scalars call this%allocate_tsp_scalars(modelname) ! - ! -- allocate members that are part of model class + ! -- allocate additional members specific to GWT model type call mem_allocate(this%inmst, 'INMST', this%memoryPath) call mem_allocate(this%indsp, 'INDSP', this%memoryPath) - call mem_allocate(this%inobs, 'INOBS', this%memoryPath) ! this%inmst = 0 this%indsp = 0 - this%inobs = 0 ! ! -- Return return @@ -1131,7 +919,6 @@ subroutine create_gwt_packages(this, indis) use SimVariablesModule, only: idm_context use GwtMstModule, only: mst_cr use GwtDspModule, only: dsp_cr - use GwtObsModule, only: gwt_obs_cr ! -- dummy class(GwtModelType) :: this integer(I4B), intent(in) :: indis @@ -1177,10 +964,6 @@ subroutine create_gwt_packages(this, indis) case ('DSP6') this%indsp = 1 mempathdsp = mempath - case ('SSM6') - this%inssm = inunit - case ('OBS6') - this%inobs = inunit case ('CNC6', 'SRC6', 'LKT6', 'SFT6', & 'MWT6', 'UZT6', 'IST6', 'API6') call expandarray(bndpkgs) @@ -1194,7 +977,6 @@ subroutine create_gwt_packages(this, indis) call mst_cr(this%mst, this%name, this%inmst, this%iout, this%fmi) call dsp_cr(this%dsp, this%name, mempathdsp, this%indsp, this%iout, & this%fmi) - call gwt_obs_cr(this%obs, this%inobs) ! ! -- Check to make sure that required ftype's have been specified call this%ftype_check(indis, this%inmst) diff --git a/src/Model/TransportModel/tsp1.f90 b/src/Model/TransportModel/tsp1.f90 index 50547148205..1ed3e16b313 100644 --- a/src/Model/TransportModel/tsp1.f90 +++ b/src/Model/TransportModel/tsp1.f90 @@ -11,12 +11,14 @@ module TransportModelModule LENMEMPATH, LENVARNAME use SimVariablesModule, only: errmsg use NumericalModelModule, only: NumericalModelType + use BndModule, only: BndType, GetBndFromList use TspIcModule, only: TspIcType use TspFmiModule, only: TspFmiType use TspAdvModule, only: TspAdvType use TspSsmModule, only: TspSsmType use TspMvtModule, only: TspMvtType use TspOcModule, only: TspOcType + use TspObsModule, only: TspObsType use BudgetModule, only: BudgetType use MatrixBaseModule @@ -28,10 +30,12 @@ module TransportModelModule type, extends(NumericalModelType) :: TransportModelType - type(TspFmiType), pointer :: fmi => null() ! flow model interface + ! Generalized transport package types common to either GWT or GWE type(TspAdvType), pointer :: adv => null() !< advection package + type(TspFmiType), pointer :: fmi => null() !< flow model interface type(TspIcType), pointer :: ic => null() !< initial conditions package type(TspMvtType), pointer :: mvt => null() !< mover transport package + type(TspObsType), pointer :: obs => null() !< observation package type(TspOcType), pointer :: oc => null() !< output control package type(TspSsmType), pointer :: ssm => null() !< source sink mixing package type(BudgetType), pointer :: budget => null() !< budget object @@ -40,6 +44,7 @@ module TransportModelModule integer(I4B), pointer :: inic => null() !< unit number IC integer(I4B), pointer :: inmvt => null() !< unit number MVT integer(I4B), pointer :: inoc => null() !< unit number OC + integer(I4B), pointer :: inobs => null() !< unit number OBS integer(I4B), pointer :: inssm => null() !< unit number SSM real(DP), pointer :: eqnsclfac => null() !< constant factor by which all terms in the model's governing equation are scaled (divided) for formulation and solution @@ -64,10 +69,16 @@ module TransportModelModule procedure, public :: tsp_cc procedure, public :: tsp_cq procedure, public :: tsp_bd + procedure, public :: tsp_ot procedure, public :: allocate_tsp_scalars procedure, public :: set_tsp_labels procedure, public :: ftype_check ! -- private + procedure, private :: tsp_ot_obs + procedure, private :: tsp_ot_flow + procedure, private :: tsp_ot_flowja + procedure, private :: tsp_ot_dv + procedure, private :: tsp_ot_bdsummary procedure, private :: create_lstfile procedure, private :: create_tsp_packages procedure, private :: log_namfile_options @@ -267,6 +278,244 @@ subroutine tsp_bd(this, icnvg, isuppress_output) return end subroutine tsp_bd + !> @brief Generalized transport model output routine + !! + !! Generalized transport model output + !< + subroutine tsp_ot(this, inmst) + ! -- modules + use TdisModule, only: kstp, kper, tdis_ot, endofperiod + ! -- dummy + class(TransportModelType) :: this + integer(I4B), intent(in) :: inmst + ! -- local + integer(I4B) :: idvsave + integer(I4B) :: idvprint + integer(I4B) :: icbcfl + integer(I4B) :: icbcun + integer(I4B) :: ibudfl + integer(I4B) :: ipflag + ! -- formats + 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 + idvprint = 0 + icbcfl = 0 + ibudfl = 0 + if (this%oc%oc_save(trim(this%depvartype))) idvsave = 1 + if (this%oc%oc_print(trim(this%depvartype))) idvprint = 1 + if (this%oc%oc_save('BUDGET')) icbcfl = 1 + if (this%oc%oc_print('BUDGET')) ibudfl = 1 + icbcun = this%oc%oc_save_unit('BUDGET') + ! + ! -- Override ibudfl and idvprint flags for nonconvergence + ! and end of period + ibudfl = this%oc%set_print_flag('BUDGET', this%icnvg, endofperiod) + idvprint = this%oc%set_print_flag(trim(this%depvartype), & + this%icnvg, endofperiod) + ! + ! Calculate and save observations + call this%tsp_ot_obs() + ! + ! Save and print flows + call this%tsp_ot_flow(icbcfl, ibudfl, icbcun, inmst) + ! + ! Save and print dependent variables + call this%tsp_ot_dv(idvsave, idvprint, ipflag) + ! + ! Print budget summaries + call this%tsp_ot_bdsummary(ibudfl, ipflag) + ! + ! -- Timing Output; if any dependendent variables or budgets + ! are printed, then ipflag is set to 1. + if (ipflag == 1) call tdis_ot(this%iout) + ! + ! -- Write non-convergence message + if (this%icnvg == 0) then + write (this%iout, fmtnocnvg) kstp, kper + end if + ! + ! -- Return + return + end subroutine tsp_ot + + !> @brief Generalized transport model output routine + !! + !! Calculate and save observations + !< + subroutine tsp_ot_obs(this) + class(TransportModelType) :: this + class(BndType), pointer :: packobj + integer(I4B) :: ip + ! -- Calculate and save observations + call this%obs%obs_bd() + call this%obs%obs_ot() + + ! -- Calculate and save package obserations + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_bd_obs() + call packobj%bnd_ot_obs() + end do + + end subroutine tsp_ot_obs + + !> @brief Generalized transport model output routine + !! + !! Save and print flows + !< + subroutine tsp_ot_flow(this, icbcfl, ibudfl, icbcun, inmst) + ! -- dummy + class(TransportModelType) :: this + integer(I4B), intent(in) :: icbcfl + integer(I4B), intent(in) :: ibudfl + integer(I4B), intent(in) :: icbcun + integer(I4B), intent(in) :: inmst + ! -- local + class(BndType), pointer :: packobj + integer(I4B) :: ip + ! + ! -- Save TSP flows + call this%tsp_ot_flowja(this%nja, this%flowja, icbcfl, icbcun) + if (this%infmi > 0) call this%fmi%fmi_ot_flow(icbcfl, icbcun) + if (this%inssm > 0) then + call this%ssm%ssm_ot_flow(icbcfl=icbcfl, ibudfl=0, icbcun=icbcun) + end if + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=0, icbcun=icbcun) + end do + + ! -- Save advanced package flows + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_package_flows(icbcfl=icbcfl, ibudfl=0) + end do + if (this%inmvt > 0) then + call this%mvt%mvt_ot_saveflow(icbcfl, ibudfl) + end if + + ! -- Print Model (GWT or GWE) flows + ! no need to print flowja + ! no need to print mst + ! no need to print fmi + if (this%inssm > 0) then + call this%ssm%ssm_ot_flow(icbcfl=icbcfl, ibudfl=ibudfl, icbcun=0) + end if + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=ibudfl, icbcun=0) + end do + + ! -- Print advanced package flows + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_package_flows(icbcfl=0, ibudfl=ibudfl) + end do + if (this%inmvt > 0) then + call this%mvt%mvt_ot_printflow(icbcfl, ibudfl) + end if + + end subroutine tsp_ot_flow + + !> @brief Generalized transport model output routine + !! + !! Write intercell flows for the transport model + !< + subroutine tsp_ot_flowja(this, nja, flowja, icbcfl, icbcun) + ! -- dummy + class(TransportModelType) :: this + integer(I4B), intent(in) :: nja + real(DP), dimension(nja), intent(in) :: flowja + integer(I4B), intent(in) :: icbcfl + integer(I4B), intent(in) :: icbcun + ! -- local + integer(I4B) :: ibinun + ! -- formats + ! + ! -- Set unit number for binary output + if (this%ipakcb < 0) then + ibinun = icbcun + elseif (this%ipakcb == 0) then + ibinun = 0 + else + ibinun = this%ipakcb + end if + if (icbcfl == 0) ibinun = 0 + ! + ! -- Write the face flows if requested + if (ibinun /= 0) then + call this%dis%record_connection_array(flowja, ibinun, this%iout) + end if + ! + ! -- Return + return + end subroutine tsp_ot_flowja + + !> @brief Generalized tranpsort model output routine + !! + !! Loop through attached packages saving and printing dependent variables + !< + subroutine tsp_ot_dv(this, idvsave, idvprint, ipflag) + class(TransportModelType) :: this + integer(I4B), intent(in) :: idvsave + integer(I4B), intent(in) :: idvprint + integer(I4B), intent(inout) :: ipflag + class(BndType), pointer :: packobj + integer(I4B) :: ip + ! + ! -- Print advanced package dependent variables + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_dv(idvsave, idvprint) + end do + + ! -- save head and print head + call this%oc%oc_ot(ipflag) + ! + ! -- Return + return + end subroutine tsp_ot_dv + + !> @brief Generalized tranpsort model output budget summary + !! + !! Loop through attached packages and write budget summaries + !< + subroutine tsp_ot_bdsummary(this, ibudfl, ipflag) + use TdisModule, only: kstp, kper, totim + class(TransportModelType) :: this + integer(I4B), intent(in) :: ibudfl + integer(I4B), intent(inout) :: ipflag + class(BndType), pointer :: packobj + integer(I4B) :: ip + ! + ! -- Package budget summary + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_bdsummary(kstp, kper, this%iout, ibudfl) + end do + + ! -- mover budget summary + if (this%inmvt > 0) then + call this%mvt%mvt_ot_bdsummary(ibudfl) + end if + + ! -- model budget summary + if (ibudfl /= 0) then + ipflag = 1 + call this%budget%budget_ot(kstp, kper, this%iout) + end if + + ! -- Write to budget csv + call this%budget%writecsv(totim) + ! + ! -- Return + return + end subroutine tsp_ot_bdsummary + !> @brief Allocate scalar variables for transport model !! !! Method to allocate memory for non-allocatable members. @@ -288,6 +537,7 @@ subroutine allocate_tsp_scalars(this, modelname) call mem_allocate(this%inadv, 'INADV', this%memoryPath) call mem_allocate(this%inssm, 'INSSM', this%memoryPath) call mem_allocate(this%inoc, 'INOC ', this%memoryPath) + call mem_allocate(this%inobs, 'INOBS', this%memoryPath) call mem_allocate(this%eqnsclfac, 'EQNSCLFAC', this%memoryPath) ! this%inic = 0 @@ -296,6 +546,7 @@ subroutine allocate_tsp_scalars(this, modelname) this%inadv = 0 this%inssm = 0 this%inoc = 0 + this%inobs = 0 this%eqnsclfac = DZERO ! ! -- Return @@ -349,6 +600,7 @@ subroutine tsp_da(this) call mem_deallocate(this%inssm) call mem_deallocate(this%inmvt) call mem_deallocate(this%inoc) + call mem_deallocate(this%inobs) call mem_deallocate(this%eqnsclfac) ! ! -- Return @@ -507,6 +759,7 @@ subroutine create_tsp_packages(this, indis) use TspSsmModule, only: ssm_cr use TspMvtModule, only: mvt_cr use TspOcModule, only: oc_cr + use TspObsModule, only: tsp_obs_cr ! -- dummy class(TransportModelType) :: this integer(I4B), intent(inout) :: indis ! DIS enabled flag @@ -569,6 +822,10 @@ subroutine create_tsp_packages(this, indis) this%inssm = inunit case ('OC6') this%inoc = inunit + case ('OBS6') + this%inobs = inunit + !case default + ! TODO end select end do ! @@ -584,7 +841,7 @@ subroutine create_tsp_packages(this, indis) call mvt_cr(this%mvt, this%name, this%inmvt, this%iout, this%fmi, & this%eqnsclfac) call oc_cr(this%oc, this%name, this%inoc, this%iout) - + call tsp_obs_cr(this%obs, this%inobs) ! ! -- Return return diff --git a/src/Model/GroundWaterTransport/gwt1obs1.f90 b/src/Model/TransportModel/tsp1obs1.f90 similarity index 56% rename from src/Model/GroundWaterTransport/gwt1obs1.f90 rename to src/Model/TransportModel/tsp1obs1.f90 index fc81782f491..55e4114bf0a 100644 --- a/src/Model/GroundWaterTransport/gwt1obs1.f90 +++ b/src/Model/TransportModel/tsp1obs1.f90 @@ -1,4 +1,4 @@ -module GwtObsModule +module TspObsModule use KindModule, only: DP, I4B use ConstantsModule, only: LINELENGTH, MAXOBSTYPES @@ -11,40 +11,37 @@ module GwtObsModule implicit none private - public :: GwtObsType, gwt_obs_cr + public :: TspObsType, tsp_obs_cr - type, extends(ObsType) :: GwtObsType + type, extends(ObsType) :: TspObsType ! -- Private members type(TspIcType), pointer, private :: ic => null() ! initial conditions real(DP), dimension(:), pointer, contiguous, private :: x => null() ! concentration real(DP), dimension(:), pointer, contiguous, private :: flowja => null() ! intercell flows contains ! -- Public procedures - procedure, public :: gwt_obs_ar - procedure, public :: obs_bd => gwt_obs_bd - procedure, public :: obs_df => gwt_obs_df - procedure, public :: obs_rp => gwt_obs_rp - procedure, public :: obs_da => gwt_obs_da + procedure, public :: tsp_obs_ar + procedure, public :: obs_bd => tsp_obs_bd + procedure, public :: obs_df => tsp_obs_df + procedure, public :: obs_rp => tsp_obs_rp + procedure, public :: obs_da => tsp_obs_da ! -- Private procedures procedure, private :: set_pointers - end type GwtObsType + end type TspObsType contains - subroutine gwt_obs_cr(obs, inobs) -! ****************************************************************************** -! gwt_obs_cr -- Create a new GwtObsType object -! Subroutine: (1) creates object -! (2) allocates pointers -! (3) initializes values -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Create a new TspObsType object + !! + !! This routine: + !! - creates an observation object + !! - allocates pointers + !! - initializes values + !< + subroutine tsp_obs_cr(obs, inobs) ! -- dummy - type(GwtObsType), pointer, intent(out) :: obs + type(TspObsType), pointer, intent(out) :: obs integer(I4B), pointer, intent(in) :: inobs -! ------------------------------------------------------------------------------ ! allocate (obs) call obs%allocate_scalars() @@ -52,22 +49,20 @@ subroutine gwt_obs_cr(obs, inobs) obs%inputFilename = '' obs%inUnitObs => inobs ! + ! -- Return return - end subroutine gwt_obs_cr + end subroutine tsp_obs_cr - subroutine gwt_obs_ar(this, ic, x, flowja) -! ****************************************************************************** -! gwt_obs_ar -- allocate and read -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Allocate and read method for package + !! + !! Method to allocate and read static data for the package. + !< + subroutine tsp_obs_ar(this, ic, x, flowja) ! -- dummy - class(GwtObsType), intent(inout) :: this + class(TspObsType), intent(inout) :: this type(TspIcType), pointer, intent(in) :: ic real(DP), dimension(:), pointer, contiguous, intent(in) :: x real(DP), dimension(:), pointer, contiguous, intent(in) :: flowja -! ------------------------------------------------------------------------------ ! ! Call ar method of parent class call this%obs_ar() @@ -75,25 +70,21 @@ subroutine gwt_obs_ar(this, ic, x, flowja) ! set pointers call this%set_pointers(ic, x, flowja) ! + ! -- Return return - end subroutine gwt_obs_ar + end subroutine tsp_obs_ar - subroutine gwt_obs_df(this, iout, pkgname, filtyp, dis) -! ****************************************************************************** -! gwt_obs_df -- define -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Define observation object + !< + subroutine tsp_obs_df(this, iout, pkgname, filtyp, dis) ! -- dummy - class(GwtObsType), intent(inout) :: this + class(TspObsType), intent(inout) :: this integer(I4B), intent(in) :: iout character(len=*), intent(in) :: pkgname character(len=*), intent(in) :: filtyp class(DisBaseType), pointer :: dis ! -- local integer(I4B) :: indx -! ------------------------------------------------------------------------------ ! ! Call overridden method of parent class call this%ObsType%obs_df(iout, pkgname, filtyp, dis) @@ -107,25 +98,21 @@ subroutine gwt_obs_df(this, iout, pkgname, filtyp, dis) ! ! -- Store obs type and assign procedure pointer for flow-ja-face observation type call this%StoreObsType('flow-ja-face', .true., indx) - this%obsData(indx)%ProcessIdPtr => gwt_process_intercell_obs_id + this%obsData(indx)%ProcessIdPtr => tsp_process_intercell_obs_id ! + ! -- Return return - end subroutine gwt_obs_df + end subroutine tsp_obs_df - subroutine gwt_obs_bd(this) -! ****************************************************************************** -! gwt_obs_bd -- save obs -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Save observations + !< + subroutine tsp_obs_bd(this) ! -- dummy - class(GwtObsType), intent(inout) :: this + class(TspObsType), intent(inout) :: this ! -- local integer(I4B) :: i, jaindex, nodenumber character(len=100) :: msg class(ObserveType), pointer :: obsrv => null() -! ------------------------------------------------------------------------------ ! call this%obs_bd_clear() ! @@ -148,72 +135,60 @@ subroutine gwt_obs_bd(this) end do end if ! + ! -- Return return - end subroutine gwt_obs_bd + end subroutine tsp_obs_bd - subroutine gwt_obs_rp(this) -! ****************************************************************************** -! gwt_obs_rp -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - class(GwtObsType), intent(inout) :: this -! ------------------------------------------------------------------------------ + !> @brief If transport model observations need checks, add them here + !< + subroutine tsp_obs_rp(this) + ! -- dummy + class(TspObsType), intent(inout) :: this ! ! Do GWT observations need any checking? If so, add checks here + ! + ! -- Return return - end subroutine gwt_obs_rp + end subroutine tsp_obs_rp - subroutine gwt_obs_da(this) -! ****************************************************************************** -! gwt_obs_da -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> Deallocate memory + !! + !! Deallocate memory associated with transport model + subroutine tsp_obs_da(this) ! -- dummy - class(GwtObsType), intent(inout) :: this -! ------------------------------------------------------------------------------ + class(TspObsType), intent(inout) :: this ! nullify (this%ic) nullify (this%x) nullify (this%flowja) call this%ObsType%obs_da() ! + ! -- Return return - end subroutine gwt_obs_da + end subroutine tsp_obs_da + !> @brief Set pointers needed by the transport OBS package + !< subroutine set_pointers(this, ic, x, flowja) -! ****************************************************************************** -! set_pointers -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy - class(GwtObsType), intent(inout) :: this + class(TspObsType), intent(inout) :: this type(TspIcType), pointer, intent(in) :: ic real(DP), dimension(:), pointer, contiguous, intent(in) :: x real(DP), dimension(:), pointer, contiguous, intent(in) :: flowja -! ------------------------------------------------------------------------------ ! this%ic => ic this%x => x this%flowja => flowja ! + ! -- Return return end subroutine set_pointers - ! -- Procedures related to GWF observations (NOT type-bound) - + !> @brief Procedure related to Tsp observations (NOT type-bound) + !! + !! Process a specific observation ID + !< subroutine gwt_process_concentration_obs_id(obsrv, dis, inunitobs, iout) -! ****************************************************************************** -! gwt_process_concentration_obs_id -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy type(ObserveType), intent(inout) :: obsrv class(DisBaseType), intent(in) :: dis @@ -223,7 +198,6 @@ subroutine gwt_process_concentration_obs_id(obsrv, dis, inunitobs, iout) integer(I4B) :: nn1 integer(I4B) :: icol, istart, istop character(len=LINELENGTH) :: ermsg, strng -! ------------------------------------------------------------------------------ ! ! -- Initialize variables strng = obsrv%IDstring @@ -242,16 +216,15 @@ subroutine gwt_process_concentration_obs_id(obsrv, dis, inunitobs, iout) call store_error_unit(inunitobs) end if ! + ! -- Return return end subroutine gwt_process_concentration_obs_id - subroutine gwt_process_intercell_obs_id(obsrv, dis, inunitobs, iout) -! ****************************************************************************** -! gwt_process_intercell_obs_id -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Procedure related to Tsp observations (NOT type-bound) + !! + !! Process an intercell observation requested by the user + !< + subroutine tsp_process_intercell_obs_id(obsrv, dis, inunitobs, iout) ! -- dummy type(ObserveType), intent(inout) :: obsrv class(DisBaseType), intent(in) :: dis @@ -263,7 +236,6 @@ subroutine gwt_process_intercell_obs_id(obsrv, dis, inunitobs, iout) character(len=LINELENGTH) :: ermsg, strng ! formats 70 format('Error: No connection exists between cells identified in text: ', a) -! ------------------------------------------------------------------------------ ! ! -- Initialize variables strng = obsrv%IDstring @@ -304,7 +276,8 @@ subroutine gwt_process_intercell_obs_id(obsrv, dis, inunitobs, iout) call store_error_unit(inunitobs) end if ! + ! -- Return return - end subroutine gwt_process_intercell_obs_id + end subroutine tsp_process_intercell_obs_id -end module GwtObsModule +end module TspObsModule diff --git a/src/meson.build b/src/meson.build index 4f1192993d7..147a1ed8cef 100644 --- a/src/meson.build +++ b/src/meson.build @@ -102,7 +102,6 @@ modflow_sources = files( 'Model' / 'GroundWaterTransport' / 'gwt1lkt1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1mst1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1mwt1.f90', - 'Model' / 'GroundWaterTransport' / 'gwt1obs1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1sft1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1src1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1uzt1.f90', @@ -131,6 +130,7 @@ modflow_sources = files( 'Model' / 'TransportModel' / 'tsp1adv1.f90', 'Model' / 'TransportModel' / 'tsp1fmi1.f90', 'Model' / 'TransportModel' / 'tsp1ic1.f90', + 'Model' / 'TransportModel' / 'tsp1obs1.f90', 'Model' / 'TransportModel' / 'tsp1oc1.f90', 'Model' / 'TransportModel' / 'tsp1mvt1.f90', 'Model' / 'TransportModel' / 'tsp1ssm1.f90',