diff --git a/autotest/test_gwf_sfr_kinematic01.py b/autotest/test_gwf_sfr_kinematic01.py new file mode 100644 index 00000000000..0bfeb38af8d --- /dev/null +++ b/autotest/test_gwf_sfr_kinematic01.py @@ -0,0 +1,185 @@ +""" +This test sfr kinematic wave approximation using example +in Ponce (1989) Example 9-3 part 1. In Ponce (1989) the +Courant number was specified, was invariant, and was used +directly in the kinematic wave equation. As a result, the +of the sfr result is approximate since the sfr kinematic +wave approximation is a numerical method where the Courant +number is a function of the solution. + +The test has 1 sfr reach that is not connected to the +groundwater model. + +The test evaluates EXT-OUTFLOW for reaches 1 against +the known numerical solution for storage weight values +set to 0.5 and 1.0. + +Ponce, V. M. (1989). Engineering Hydrology, Principles and Practices. +""" + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +paktest = "sfr" +cases = ("sfr-kwe01", "sfr-kwe02") +storage_weights = (0.5, 1.0) + + +def build_models(idx, test): + name = cases[idx] + + dt = 60.0 * 60.0 + flows = np.array( + [ + 0.0, + 30.0, + 60.0, + 90.0, + 120.0, + 150.0, + 120.0, + 90.0, + 60.0, + 30.0, + 0.0, + 0.0, + 0.0, + 0.0, + ] + ) + times = np.arange(0.0, (flows.shape[0]) * dt, dt) + + # static model data + # temporal discretization + nper = times.shape[0] - 1 + nstp = 1 + perioddata = [(dt, nstp, 1.0) for idx in range(nper)] + + # spatial discretization data + nlay, nrow, ncol = 1, 1, 6 + delr, delc = 100.0, 100.0 + dx = 7200.0 + top = 100.0 + botm = 0.0 + + sim = flopy.mf6.MFSimulation( + sim_name=name, + sim_ws=test.workspace, + version="mf6", + exe_name="mf6", + ) + tdis = flopy.mf6.ModflowTdis( + sim, time_units="seconds", nper=nper, perioddata=perioddata + ) + ims = flopy.mf6.ModflowIms(sim) + + gwf = flopy.mf6.ModflowGwf(sim, modelname=name) + dis = flopy.mf6.ModflowGwfdis( + gwf, + length_units="meters", + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=dx, + delc=dx, + top=top, + botm=botm, + ) + npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1) + ic = flopy.mf6.ModflowGwfic(gwf, strt=top) + sto = flopy.mf6.ModflowGwfsto( + gwf, iconvert=1, ss=1e-6, sy=0.2, transient={0: True} + ) + ghb = flopy.mf6.ModflowGwfghb( + gwf, stress_period_data=[(0, 0, 0, top, 5.0)] + ) + oc = flopy.mf6.ModflowGwfoc(gwf, printrecord=[("budget", "all")]) + + # sfr file + nreaches = 1 + slope = 1.0 / dx + roughness = 0.03574737676661647 + + # [] [] + pak_data = [ + ( + ifno, + -1, + -1, + -1, + dx, + 10.0, + slope, + top, + 1.0, + 0.0, + roughness, + 0, + 0.0, + 0, + ) + for ifno in range(nreaches) + ] + sfr_spd = {idx: [(0, "inflow", q)] for idx, q in enumerate(flows[1:])} + sfr = flopy.mf6.ModflowGwfsfr( + gwf, + print_flows=True, + storage=True, + storage_weight=storage_weights[idx], + nreaches=nreaches, + packagedata=pak_data, + connectiondata=[(0,)], + perioddata=sfr_spd, + pname="sfr-1", + ) + fname = f"{name}.sfr.obs" + sfr_obs = { + f"{fname}.csv": [ + ("inflow", "ext-inflow", (0,)), + ("outflow", "ext-outflow", (0,)), + ] + } + sfr.obs.initialize(filename=fname, print_input=True, continuous=sfr_obs) + + return sim + + +def check_output(idx, test): + print("Checking sfr external outflow") + name = cases[idx] + obs_values = flopy.utils.Mf6Obs( + test.workspace / f"{name}.sfr.obs.csv" + ).get_data() + # fmt: off + test_values = {storage_weights[0]: np.array( + [ + 0. , -10.68920803, -58.58988171, -92.59647491, -125.23616408, + -144.76383592, -117.40352509, -89.57181956, -64.81693279, -53.46056541, + -11.22910318, -4.95517783, -2.72890377 + ] + ), + storage_weights[1]: np.array( + [ + 0. , -15.3918031 , -58.58988171, -92.59647491, -125.23616408, + -144.76383592, -117.40352509, -88.84252125, -62.90365264, -47.01932528, + -19.87588927, -9.96585666, -5.63155787 + ] + ), + } + assert np.allclose( + obs_values["OUTFLOW"], test_values[storage_weights[idx]] + ), f"failed comparison for '{name}' observation" + + +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 7ff4bd0b791..c3b6c4e497e 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -1,5 +1,5 @@ -% Use this template for starting initializing the release notes -% after a release has just been made. +% Use this template for starting initializing the release notes after a release +% has just been made. \item \currentmodflowversion @@ -7,7 +7,8 @@ \begin{itemize} \item A new adaptive time stepping (ATS) capability was added to the Advection (ADV) Package of the Groundwater Transport (GWT) Model. A new input option, called ATS\_PERCEL, specifies the fractional cell distance that a particle of water can travel within one time step. When ATS\_PERCEL is specified by the user, and the ATS utility is activated in the TDIS Package, the ADV Package will calculate the largest time step that will meet this fractional cell distance constraint, and will submit this time step to the ATS utility. This option may improve time stepping for solute transport models and for variable-density flow and transport models by allowing step lengths to be calculated as a function of the flow system rather than being specified as input by the user. \item Added the capability to write sorbate concentrations to a binary output file. A new SORBATE option is now available in the Mass Storage and Transfer (MST) Package of the GWT Model to provide the name of the binary output file for sorbate concentrations. Sorbate concentrations will be written to the binary output file whenever concentrations for the GWT model are saved, as determined by the GWT Output Control option. - % \item xxx + \item Add kinematic-wave routing option for the Streamflow Routing (SFR) Package. Prior to this change, the SFR Package could only simulate unidirectional, steady, uniform flow. With kinematic-wave routing, unidirectional waves can now propagate through the SFR network by explicitly including a storage term in the reach continuity equation. The kinematic-wave routing option is based on the ``TRANSROUTE'' option available in the SFR Package in MODFLOW-NWT. The kinematic-wave routing option is enabled by specifying ``STORAGE'' in the SFR Package OPTIONS block. + % \item xxx \end{itemize} \underline{EXAMPLES} diff --git a/doc/mf6io/gwf/sfr.tex b/doc/mf6io/gwf/sfr.tex index 7c873e6de3f..80bcdfa344e 100644 --- a/doc/mf6io/gwf/sfr.tex +++ b/doc/mf6io/gwf/sfr.tex @@ -28,6 +28,10 @@ \subsubsection{Structure of Blocks} \noindent \textit{IF ndv IS GREATER THAN ZERO FOR ANY REACH} \lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-sfr-diversions.dat} +\vspace{5mm} +\noindent \textit{INITIALSTAGES BLOCK IS OPTIONAL} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-sfr-initialstages.dat} + \vspace{5mm} \noindent \textit{FOR ANY STRESS PERIOD} \lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-sfr-period.dat} diff --git a/doc/mf6io/ims.tex b/doc/mf6io/ims.tex index 9efed16d1b9..f3cf7c565d2 100644 --- a/doc/mf6io/ims.tex +++ b/doc/mf6io/ims.tex @@ -12,7 +12,7 @@ \subsection{Explanation of Variables} \input{./mf6ivar/tex/sln-ims-desc.tex} \end{description} -\noindent \textbf{IMS Default and Specified Complexity Values} +\subsection{IMS Default and Specified Complexity Values} The values that are assigned to the \texttt{nonlinear} and \texttt{linear} variables for the the \texttt{simple}, \texttt{moderate}, and \texttt{complex} complexity options are summarized in table~\ref{table:imsopt}. The values defined for the \texttt{simple} complexity option are assigned if the \texttt{COMPLEXITY} keyword is not specified in the \texttt{OPTIONS} block. diff --git a/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn b/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn index eb3967e441d..270a66bd75e 100644 --- a/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn @@ -2,6 +2,22 @@ # flopy multi-package # package-type advanced-stress-package +block options +name storage +type keyword +reader urword +optional true +longname activate reach storage +description keyword that activates storage contributions to the stream-flow routing package continuity equation. + +block options +name storage_weight +type double precision +reader urword +optional true +longname reach storage time weighting +description real number value that defines the time weighting factor used to calculate the change in channel storage. STORAGE\_WEIGHT must have a value between 0.5 and 1. Default STORAGE\_WEIGHT value is 0.5. + block options name auxiliary type string @@ -641,6 +657,41 @@ reader urword longname iprior code description character string value that defines the the prioritization system for the diversion, such as when insufficient water is available to meet all diversion stipulations, and is used in conjunction with the value of FLOW value specified in the STRESS\_PERIOD\_DATA section. Available diversion options include: (1) CPRIOR = `FRACTION', then the amount of the diversion is computed as a fraction of the streamflow leaving reach IFNO ($Q_{DS}$); in this case, 0.0 $\le$ DIVFLOW $\le$ 1.0. (2) CPRIOR = `EXCESS', a diversion is made only if $Q_{DS}$ for reach IFNO exceeds the value of DIVFLOW. If this occurs, then the quantity of water diverted is the excess flow ($Q_{DS} -$ DIVFLOW) and $Q_{DS}$ from reach IFNO is set equal to DIVFLOW. This represents a flood-control type of diversion, as described by Danskin and Hanson (2002). (3) CPRIOR = `THRESHOLD', then if $Q_{DS}$ in reach IFNO is less than the specified diversion flow DIVFLOW, no water is diverted from reach IFNO. If $Q_{DS}$ in reach IFNO is greater than or equal to DIVFLOW, DIVFLOW is diverted and $Q_{DS}$ is set to the remainder ($Q_{DS} -$ DIVFLOW)). This approach assumes that once flow in the stream is sufficiently low, diversions from the stream cease, and is the `priority' algorithm that originally was programmed into the STR1 Package (Prudic, 1989). (4) CPRIOR = `UPTO' -- if $Q_{DS}$ in reach IFNO is greater than or equal to the specified diversion flow DIVFLOW, $Q_{DS}$ is reduced by DIVFLOW. If $Q_{DS}$ in reach IFNO is less than DIVFLOW, DIVFLOW is set to $Q_{DS}$ and there will be no flow available for reaches connected to downstream end of reach IFNO. +# --------------------- gwf initial stages --------------------- + +block initialstages +name initialstages +type recarray ifno initialstage +shape (maxbound) +valid +optional false +reader urword +longname +description + +block initialstages +name ifno +type integer +shape +tagged false +in_record true +optional false +reader urword +longname reach number for this entry +description integer value that defines the feature (reach) number associated with the specified initial stage. Initial stage data must be specified for every reach or the program will terminate with an error. The program will also terminate with a error if IFNO is less than one or greater than NREACHES. +numeric_index true + +block initialstages +name initialstage +type double precision +shape +tagged false +in_record true +optional false +reader urword +longname initial reach stage +description real value that defines the initial stage for the reach. The program will terminate with an error if INITIALSTAGE is less than the RTP value for reach IFNO defined in the PACKAGEDATA block. INITIALSTAGE data are used only if STORAGE is specified in the Options block and the first stress period is transient or for reaches defined to use the SIMPLE STATUS in the Period block. + # --------------------- gwf sfr period --------------------- diff --git a/make/makefile b/make/makefile index 29312d15bae..71bc053718f 100644 --- a/make/makefile +++ b/make/makefile @@ -505,6 +505,8 @@ $(OBJDIR)/sparsekit.o \ $(OBJDIR)/rcm.o \ $(OBJDIR)/blas1_d.o \ $(OBJDIR)/Iunit.o \ +$(OBJDIR)/GwfSfrCommon.o \ +$(OBJDIR)/gwf-sfr-transient.o \ $(OBJDIR)/gwf-sfr-steady.o \ $(OBJDIR)/gwf-sfr-constant.o \ $(OBJDIR)/RectangularGeometry.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index e2e0e874fab..cc81f11cd98 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -261,7 +261,8 @@ - + + @@ -313,6 +314,7 @@ + diff --git a/src/Model/GroundWaterFlow/gwf-sfr.f90 b/src/Model/GroundWaterFlow/gwf-sfr.f90 index e8fd8a5d7b2..7e06f2effc1 100644 --- a/src/Model/GroundWaterFlow/gwf-sfr.f90 +++ b/src/Model/GroundWaterFlow/gwf-sfr.f90 @@ -59,6 +59,7 @@ module SfrModule character(len=LENBOUNDNAME), dimension(:), pointer, & contiguous :: sfrname => null() !< internal SFR reach name ! -- integers + integer(I4B), pointer :: istorage => null() !< flag for using kinematic wave approximation integer(I4B), pointer :: iprhed => null() !< flag for printing stages to listing file integer(I4B), pointer :: istageout => null() !< flag and unit number for binary stage output integer(I4B), pointer :: ibudgetout => null() !< flag and unit number for binary sfr budget output @@ -80,6 +81,7 @@ module SfrModule real(DP), pointer :: timeconv => NULL() !< time conversion factor (SI to model units) real(DP), pointer :: dmaxchg => NULL() !< maximum depth change allowed real(DP), pointer :: deps => NULL() !< perturbation value + real(DP), pointer :: storage_weight => NULL() !< time weighting factor for kinematic wave approximation ! -- integer vectors integer(I4B), dimension(:), pointer, contiguous :: isfrorder => null() !< sfr reach order determined from DAG of upstream reaches integer(I4B), dimension(:), pointer, contiguous :: ia => null() !< CRS row pointer for SFR reaches @@ -113,13 +115,18 @@ module SfrModule integer(I4B), dimension(:), pointer, contiguous :: ndiv => null() !< number of diversions for each reach real(DP), dimension(:), pointer, contiguous :: usflow => null() !< upstream reach flow real(DP), dimension(:), pointer, contiguous :: dsflow => null() !< downstream reach flow + real(DP), dimension(:), pointer, contiguous :: dsflowold => null() !< downstream reach flow for previous time step + real(DP), dimension(:), pointer, contiguous :: usinflow => null() !< upstream reach flow for previous time step + real(DP), dimension(:), pointer, contiguous :: usinflowold => null() !< upstream reach flow for previous time step real(DP), dimension(:), pointer, contiguous :: depth => null() !< reach depth real(DP), dimension(:), pointer, contiguous :: stage => null() !< reach stage + real(DP), dimension(:), pointer, contiguous :: stageold => null() !< reach stage for last timestep real(DP), dimension(:), pointer, contiguous :: gwflow => null() !< flow from groundwater to reach real(DP), dimension(:), pointer, contiguous :: simevap => null() !< simulated reach evaporation real(DP), dimension(:), pointer, contiguous :: simrunoff => null() !< simulated reach runoff real(DP), dimension(:), pointer, contiguous :: stage0 => null() !< previous reach stage iterate real(DP), dimension(:), pointer, contiguous :: usflow0 => null() !< previous upstream reach flow iterate + real(DP), dimension(:), pointer, contiguous :: storage => null() !< previous upstream reach flow iterate ! -- cross-section data integer(I4B), pointer :: ncrossptstot => null() !< total number of cross-section points integer(I4B), dimension(:), pointer, contiguous :: ncrosspts => null() !< number of cross-section points for each reach @@ -183,6 +190,7 @@ module SfrModule procedure, private :: sfr_set_stressperiod procedure, private :: sfr_solve procedure, private :: sfr_calc_constant + procedure, private :: sfr_calc_transient procedure, private :: sfr_calc_steady procedure, private :: sfr_update_flows procedure, private :: sfr_adjust_ro_ev @@ -204,14 +212,17 @@ module SfrModule procedure, private :: sfr_read_crossection procedure, private :: sfr_read_connectiondata procedure, private :: sfr_read_diversions + procedure, private :: sfr_read_initial_stages ! -- calculations procedure, private :: sfr_calc_reach_depth procedure, private :: sfr_calc_xs_depth ! -- error checking procedure, private :: sfr_check_conversion + procedure, private :: sfr_check_storage_weight procedure, private :: sfr_check_reaches procedure, private :: sfr_check_connections procedure, private :: sfr_check_diversions + procedure, private :: sfr_check_initialstages procedure, private :: sfr_check_ustrf ! -- budget procedure, private :: sfr_setup_budobj @@ -244,6 +255,25 @@ module subroutine sfr_calc_steady(this, n, d1, hgwf, & end subroutine end interface + interface + module subroutine sfr_calc_transient(this, n, d1, hgwf, & + qu, qi, qfrommvr, qr, qe, qro, & + qgwf, qd) + class(SfrType) :: this !< SfrType object + integer(I4B), intent(in) :: n !< reach number + real(DP), intent(inout) :: d1 !< current reach depth estimate + real(DP), intent(in) :: hgwf !< head in gw cell + real(DP), intent(in) :: qu !< reach upstream flow + real(DP), intent(in) :: qi !< reach specified inflow + real(DP), intent(in) :: qfrommvr !< reach flow from mover + real(DP), intent(in) :: qr !< reach rainfall + real(DP), intent(in) :: qe !< reach evaporation + real(DP), intent(in) :: qro !< reach runoff flow + real(DP), intent(inout) :: qgwf !< reach-aquifer exchange + real(DP), intent(inout) :: qd !< reach outflow + end subroutine + end interface + interface module subroutine sfr_calc_constant(this, n, d1, hgwf, qgwf, qd) class(SfrType) :: this !< SfrType object @@ -320,6 +350,7 @@ subroutine sfr_allocate_scalars(this) call this%BndType%allocate_scalars() ! ! -- allocate the object and assign values to object variables + call mem_allocate(this%istorage, 'ISTORAGE', this%memoryPath) call mem_allocate(this%iprhed, 'IPRHED', this%memoryPath) call mem_allocate(this%istageout, 'ISTAGEOUT', this%memoryPath) call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath) @@ -335,6 +366,7 @@ subroutine sfr_allocate_scalars(this) call mem_allocate(this%timeconv, 'TIMECONV', this%memoryPath) call mem_allocate(this%dmaxchg, 'DMAXCHG', this%memoryPath) call mem_allocate(this%deps, 'DEPS', this%memoryPath) + call mem_allocate(this%storage_weight, 'STORAGE_WEIGHT', this%memoryPath) call mem_allocate(this%nconn, 'NCONN', this%memoryPath) call mem_allocate(this%icheck, 'ICHECK', this%memoryPath) call mem_allocate(this%iconvchk, 'ICONVCHK', this%memoryPath) @@ -346,6 +378,7 @@ subroutine sfr_allocate_scalars(this) call mem_setptr(this%gwfiss, 'ISS', create_mem_path(this%name_model)) ! ! -- Set values + this%istorage = 0 this%iprhed = 0 this%istageout = 0 this%ibudgetout = 0 @@ -361,6 +394,7 @@ subroutine sfr_allocate_scalars(this) this%timeconv = DNODATA this%dmaxchg = DEM5 this%deps = DP999 * this%dmaxchg + this%storage_weight = DNODATA this%nconn = 0 this%icheck = 1 this%iconvchk = 1 @@ -420,6 +454,20 @@ subroutine sfr_allocate_arrays(this) call mem_allocate(this%stage0, this%maxbound, 'STAGE0', this%memoryPath) call mem_allocate(this%usflow0, this%maxbound, 'USFLOW0', this%memoryPath) ! + ! -- stage, usflow, inflow, and dsflow for previous timestep + if (this%istorage == 1) then + call mem_allocate(this%stageold, this%maxbound, 'STAGEOLD', & + this%memoryPath) + call mem_allocate(this%usinflow, this%maxbound, 'USINFLOW', & + this%memoryPath) + call mem_allocate(this%usinflowold, this%maxbound, 'USINFLOWOLD', & + this%memoryPath) + call mem_allocate(this%dsflowold, this%maxbound, 'DSFLOWOLD', & + this%memoryPath) + call mem_allocate(this%storage, this%maxbound, 'STORAGE', & + this%memoryPath) + end if + ! ! -- reach order and connection data call mem_allocate(this%isfrorder, this%maxbound, 'ISFRORDER', & this%memoryPath) @@ -485,6 +533,15 @@ subroutine sfr_allocate_arrays(this) this%stage0(i) = DZERO this%usflow0(i) = DZERO ! + ! -- stage + if (this%istorage == 1) then + this%stageold(i) = DZERO + this%usinflow(i) = DZERO + this%usinflowold(i) = DZERO + this%dsflowold(i) = DZERO + this%storage(i) = DZERO + end if + ! ! -- boundary data this%rough(i) = DZERO this%rain(i) = DZERO @@ -641,6 +698,9 @@ subroutine sfr_read_dimensions(this) ! -- read diversion data call this%sfr_read_diversions() ! + ! -- read initial stage data + call this%sfr_read_initial_stages() + ! ! -- setup the budget object call this%sfr_setup_budobj() ! @@ -682,10 +742,27 @@ subroutine sfr_options(this, option, found) character(len=*), parameter :: fmtsfrbin = & "(4x, 'SFR ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, /4x, & &'OPENED ON UNIT: ', I0)" + character(len=*), parameter :: fmtstoweight = & + &"(4x, 'KINEMATIC STORAGE WEIGHT (',g0,') SPECIFIED.')" ! ! -- Check for SFR options found = .true. select case (option) + case ('STORAGE') + this%istorage = 1 + write (this%iout, '(4x,a)') trim(adjustl(this%text))// & + ' REACH STORAGE IS ACTIVE.' + case ('STORAGE_WEIGHT') + r = this%parser%GetDouble() + if (r < DHALF .or. r > DONE) then + write (errmsg, '(a,g0,a)') & + "STORAGE_WEIGHT SPECIFIED TO BE '", r, & + "' BUT CANNOT BE LESS THAN 0.5 OR GREATER THAN 1.0" + call store_error(errmsg) + else + this%storage_weight = r + write (this%iout, fmtstoweight) this%storage_weight + end if case ('PRINT_STAGE') this%iprhed = 1 write (this%iout, '(4x,a)') trim(adjustl(this%text))// & @@ -838,6 +915,9 @@ subroutine sfr_ar(this) ! -- check the sfr unit conversion data call this%sfr_check_conversion() ! + ! -- check the storage_weight + call this%sfr_check_storage_weight() + ! ! -- check the sfr reach data call this%sfr_check_reaches() @@ -848,6 +928,11 @@ subroutine sfr_ar(this) if (this%idiversions /= 0) then call this%sfr_check_diversions() end if + + ! -- check the diversion data + if (this%istorage == 1) then + call this%sfr_check_initialstages() + end if ! ! -- terminate if errors were detected in any of the static sfr data ierr = count_errors() @@ -1745,6 +1830,105 @@ subroutine sfr_read_diversions(this) return end subroutine sfr_read_diversions + !> @ brief Read initialstages data for the package + !! + !! Method to read initialstages data for each reach for the SFR package. + !! + !< + subroutine sfr_read_initial_stages(this) + ! -- modules + use TimeSeriesManagerModule, only: read_value_or_time_series_adv + ! -- dummy variables + class(SfrType), intent(inout) :: this !< SfrType object + ! -- local variables + integer(I4B) :: n + integer(I4B) :: ierr + logical(LGP) :: isfound + logical(LGP) :: endOfBlock + integer(I4B) :: i + real(DP) :: rval + integer, allocatable, dimension(:) :: nboundchk + ! + ! -- read data + call this%parser%GetBlock('INITIALSTAGES', isfound, ierr, & + supportOpenClose=.true., & + blockRequired=.false.) + ! + ! -- parse block if detected + if (isfound) then + write (this%iout, '(/1x,a)') & + 'PROCESSING '//trim(adjustl(this%text))//' INITIALSTAGES' + + allocate (nboundchk(this%maxbound)) + do n = 1, this%maxbound + nboundchk(n) = 0 + end do + + do + call this%parser%GetNextLine(endOfBlock) + if (endOfBlock) exit + + ! -- read reach number + n = this%parser%GetInteger() + + if (n < 1 .or. n > this%maxbound) then + write (errmsg, '(a,i0,a,1x,i0,a)') & + 'Reach number (', n, ') must be greater than 0 and less & + &than or equal to', this%maxbound, '.' + call store_error(errmsg) + cycle + end if + + ! -- increment nboundchk + nboundchk(n) = nboundchk(n) + 1 + + rval = this%parser%GetDouble() + this%stage(n) = rval + this%depth(n) = rval - this%strtop(n) + + if (rval < this%strtop(n)) then + write (errmsg, '(a,g0,a,1x,i0,1x,a,g0,a)') & + 'Initial stage (', rval, ') for reach', n, & + 'is less than the reach top (', this%strtop(n), ').' + call store_error(errmsg) + end if + end do + + write (this%iout, '(1x,a)') & + 'END OF '//trim(adjustl(this%text))//' INITIALSTAGES' + + ! + ! -- Check to make sure that every reach is specified and that no reach + ! is specified more than once. + do i = 1, this%maxbound + if (nboundchk(i) == 0) then + write (errmsg, '(a,i0,1x,a)') & + 'Information for reach ', i, 'not specified in initialstages block.' + call store_error(errmsg) + else if (nboundchk(i) > 1) then + write (errmsg, '(a,1x,i0,1x,a,1x,i0)') & + 'Initial stage information specified', & + nboundchk(i), 'times for reach', i + call store_error(errmsg) + end if + end do + deallocate (nboundchk) + else + ! -- set default initial stage based on a zero depth + if (this%istorage == 1) then + do n = 1, this%maxbound + rval = this%strtop(n) + this%stage(n) = rval + end do + end if + end if + ! + ! -- terminate if errors encountered in reach block + if (count_errors() > 0) then + call this%parser%StoreErrorUnit() + end if + end subroutine sfr_read_initial_stages + !> @ brief Read and prepare period data for package !! !! Method to read and prepare period data for the SFR package. @@ -1944,6 +2128,15 @@ subroutine sfr_ad(this) ! -- local variables integer(I4B) :: n integer(I4B) :: iaux + + ! -- update previous values + if (this%istorage == 1) then + do n = 1, this%maxbound + this%stageold(n) = this%stage(n) + this%usinflowold(n) = this%usinflow(n) + this%dsflowold(n) = this%dsflow(n) + end do + end if ! ! -- Most advanced package AD routines have to restore state if ! the solution failed and the time step is being retried with a smaller @@ -2707,6 +2900,15 @@ subroutine sfr_da(this) call mem_deallocate(this%denseterms) call mem_deallocate(this%viscratios) ! + ! -- stage, usflow, and dsflow for previous timestep + if (this%istorage == 1) then + call mem_deallocate(this%stageold) + call mem_deallocate(this%dsflowold) + call mem_deallocate(this%storage) + call mem_deallocate(this%usinflow) + call mem_deallocate(this%usinflowold) + end if + ! ! -- deallocate reach order and connection data call mem_deallocate(this%isfrorder) call mem_deallocate(this%ia) @@ -2765,6 +2967,8 @@ subroutine sfr_da(this) end if ! ! -- deallocate scalars + call mem_deallocate(this%istorage) + call mem_deallocate(this%storage_weight) call mem_deallocate(this%iprhed) call mem_deallocate(this%istageout) call mem_deallocate(this%ibudgetout) @@ -3509,9 +3713,15 @@ subroutine sfr_solve(this, n, h, hcof, rhs, update) if (this%iboundpak(n) < 0) then call this%sfr_calc_constant(n, d1, hgwf, qgwf, qd) else - call this%sfr_calc_steady(n, d1, hgwf, qu, qi, & - qfrommvr, qr, qe, qro, & - qgwf, qd) + if (this%gwfiss == 0 .and. this%istorage == 1) then + call this%sfr_calc_transient(n, d1, hgwf, qu, qi, & + qfrommvr, qr, qe, qro, & + qgwf, qd) + else + call this%sfr_calc_steady(n, d1, hgwf, qu, qi, & + qfrommvr, qr, qe, qro, & + qgwf, qd) + end if end if ! -- update sfr stage @@ -4183,6 +4393,31 @@ subroutine sfr_check_conversion(this) return end subroutine sfr_check_conversion + !> @brief Check storage weight + !! + !! Method to check the kinematic storage weight for a SFR package. + !! If the kinematic storage weight has not been set it is set to + !! the default value. + !! + !< + subroutine sfr_check_storage_weight(this) + ! -- dummy variables + class(SfrType) :: this !< SfrType object + ! -- formats + character(len=*), parameter :: fmtweight = & + &"(1x,'SFR PACKAGE (',a,') SETTING DEFAULT',& + &/4x,'STORAGE_WEIGHT VALUE (',g0,').',/)" + ! + ! -- set storage weight if it has not been defined yet + if (this%istorage == 1) then + if (this%storage_weight == DNODATA) then + this%storage_weight = DHALF + write (this%iout, fmtweight) & + trim(adjustl(this%packName)), this%storage_weight + end if + end if + end subroutine sfr_check_storage_weight + !> @brief Check reach data !! !! Method to check specified data for a SFR package. This method @@ -4696,6 +4931,63 @@ subroutine sfr_check_diversions(this) return end subroutine sfr_check_diversions + !> @brief Check initial stage data + !! + !! Method to check initial data for a SFR package and calculates + !! the initial upstream and downstream flows for the reach based + !! on the initial staalso creates the tables used to print input + !! data, if this option in enabled in the SFR package. + !! + !< + subroutine sfr_check_initialstages(this) + class(SfrType) :: this !< SfrType object + + character(len=LINELENGTH) :: title + character(len=LINELENGTH) :: text + character(len=5) :: crch + integer(I4B) :: n + real(DP) :: qman + + ! skip check if storage is not activated + if (this%istorage == 0) return + + ! write header + if (this%iprpak /= 0) then + ! + ! -- reset the input table object + title = trim(adjustl(this%text))//' PACKAGE ('// & + trim(adjustl(this%packName))//') REACH INITIAL STAGE DATA' + call table_cr(this%inputtab, this%packName, title) + call this%inputtab%table_df(this%maxbound, 4, this%iout) + text = 'REACH' + call this%inputtab%initialize_column(text, 10, alignment=TABCENTER) + text = 'INITIAL STAGE' + call this%inputtab%initialize_column(text, 10, alignment=TABCENTER) + text = 'INITIAL DEPTH' + call this%inputtab%initialize_column(text, 10, alignment=TABCENTER) + text = 'INITIAL FLOW' + call this%inputtab%initialize_column(text, 10, alignment=TABCENTER) + end if + ! + ! -- check that data are correct + do n = 1, this%maxbound + write (crch, '(i5)') n + + ! calculate the initial flows + call this%sfr_calc_qman(n, this%depth(n), qman) + this%usinflow(n) = qman + this%dsflow(n) = qman + + ! add terms to the table + if (this%iprpak /= 0) then + call this%inputtab%add_term(n) + call this%inputtab%add_term(this%stage(n)) + call this%inputtab%add_term(this%depth(n)) + call this%inputtab%add_term(qman) + end if + end do + end subroutine sfr_check_initialstages + !> @brief Check upstream fraction data !! !! Method to check upstream fraction data for a SFR package. @@ -5300,6 +5592,9 @@ subroutine sfr_fill_budobj(this) d = this%depth(n) a = this%calc_surface_area_wet(n, d) this%qauxcbc(1) = a * d + if (this%gwfiss == 0 .and. this%istorage == 1) then + q = this%storage(n) + end if else q = DZERO this%qauxcbc(1) = DZERO diff --git a/src/Model/GroundWaterFlow/submodules/gwf-sfr-transient.f90 b/src/Model/GroundWaterFlow/submodules/gwf-sfr-transient.f90 new file mode 100644 index 00000000000..921f2e9d7b3 --- /dev/null +++ b/src/Model/GroundWaterFlow/submodules/gwf-sfr-transient.f90 @@ -0,0 +1,190 @@ +submodule(SfrModule) SfrModuleTransient +contains + + !> @brief Method for solving for transient reaches + !! + !! Method to solve the continuity equation for a SFR package + !! reach using the kinematic wave approximation. + !! + !< + module procedure sfr_calc_transient + use TdisModule, only: delt + + integer(I4B) :: igwfconn + integer(I4B) :: number_picard + integer(I4B) :: i + integer(I4B) :: j + real(DP) :: kinematic_residual + real(DP) :: kinematic_storage + real(DP) :: weight + real(DP) :: celerity + real(DP) :: courant + real(DP) :: dq + real(DP) :: qtol + real(DP) :: qsrc + real(DP) :: qlat + real(DP) :: da + real(DP) :: db + real(DP) :: dc + real(DP) :: dd + real(DP) :: qa + real(DP) :: qb + real(DP) :: qc + real(DP) :: xsa_a + real(DP) :: xsa_b + real(DP) :: xsa_c + real(DP) :: q + real(DP) :: q2 + real(DP) :: d + real(DP) :: d2 + real(DP) :: a + real(DP) :: a2 + real(DP) :: d1old + real(DP) :: qd2 + real(DP) :: ad + real(DP) :: ad2 + real(DP) :: residual + real(DP) :: residual2 + real(DP) :: residual_final + real(DP) :: qderv + real(DP) :: delq + real(DP) :: delh + real(DP) :: dd2 + + weight = this%storage_weight + dq = this%deps + qtol = dq * DTWO + + celerity = DZERO + qgwf = DZERO + + qlat = qr + qro - qe + + this%usinflow(n) = qu + qi + qfrommvr + + qa = this%usinflowold(n) + qb = this%dsflowold(n) + call this%sfr_calc_reach_depth(n, qa, da) + call this%sfr_calc_reach_depth(n, qb, db) + + qc = this%usinflow(n) + call this%sfr_calc_reach_depth(n, qc, dc) + + xsa_a = this%calc_area_wet(n, da) + xsa_b = this%calc_area_wet(n, db) + xsa_c = this%calc_area_wet(n, dc) + + ! estimate qd + qd = this%dsflow(n) + if (qd == DZERO) then + qd = (qc + qb) * DHALF + end if + call this%sfr_calc_reach_depth(n, qd, dd) + ad = this%calc_area_wet(n, dd) + + ! estimate the depth at the midpoint + d1 = (dc + dd) * DHALF + d1old = d1 + + ! estimate qgwf + igwfconn = this%sfr_gwf_conn(n) + if (igwfconn == 1) then + q = qu + qi + qr - qe + qro + qfrommvr + call this%sfr_calc_qgwf(n, d1, hgwf, qgwf) + qgwf = -qgwf + if (qgwf > q) then + qgwf = q + end if + end if + + ! calculate maximum wave speed and courant number + q = qc + qlat - qgwf + call this%sfr_calc_reach_depth(n, q, d) + a = this%calc_area_wet(n, d) + if (d > DZERO) then + q2 = q + dq + call this%sfr_calc_reach_depth(n, q2, d2) + a2 = this%calc_area_wet(n, d2) + celerity = (q2 - q) / (a2 - a) + courant = celerity * delt / this%length(n) + ! write(*,*) this%length(n), courant * delt + end if + + qlat = qlat / this%length(n) + + number_picard = this%maxsfrpicard + if (igwfconn == 1) then + number_picard = this%maxsfrpicard + else + number_picard = 1 + end if + + kinematicpicard: do i = 1, number_picard + if (igwfconn == 1) then + q = qu + qi + qr - qe + qro + qfrommvr + call this%sfr_calc_qgwf(n, d1, hgwf, qgwf) + qgwf = -qgwf + if (qgwf > q) then + qgwf = q + end if + end if + + qsrc = qlat - qgwf / this%length(n) + + newton: do j = 1, this%maxsfrit + qd2 = qd + dq + call this%sfr_calc_reach_depth(n, qd2, dd2) + ad2 = this%calc_area_wet(n, dd2) + + residual = kinematic_residual(qa, qb, qc, qd, & + xsa_a, xsa_b, xsa_c, ad, & + qsrc, this%length(n), weight, delt, & + courant) + + residual2 = kinematic_residual(qa, qb, qc, qd2, & + xsa_a, xsa_b, xsa_c, ad2, & + qsrc, this%length(n), weight, delt, & + courant) + qderv = (residual2 - residual) / dq + if (qderv > DZERO) then + delq = -residual / qderv + else + delq = DZERO + end if + + if (qd + delq < DEM30) then + delq = -qd + end if + + qd = qd + delq + + call this%sfr_calc_reach_depth(n, qd, dd) + ad = this%calc_area_wet(n, dd) + residual_final = kinematic_residual(qa, qb, qc, qd, & + xsa_a, xsa_b, xsa_c, ad, & + qsrc, this%length(n), weight, delt, & + courant) + + if (abs(delq) < qtol .and. abs(residual_final) < qtol) then + exit newton + end if + + end do newton + + qd = max(qd, DZERO) + d1 = (dc + dd) * DHALF + delh = (d1 - d1old) + + if (i > 1 .and. abs(delh) < this%dmaxchg) then + exit kinematicpicard + end if + + end do kinematicpicard + + this%storage(n) = kinematic_storage(xsa_a, xsa_b, xsa_c, ad, & + this%length(n), delt, & + courant) + + end procedure sfr_calc_transient + +end submodule diff --git a/src/Model/ModelUtilities/GwfSfrCommon.f90 b/src/Model/ModelUtilities/GwfSfrCommon.f90 new file mode 100644 index 00000000000..32a2957233f --- /dev/null +++ b/src/Model/ModelUtilities/GwfSfrCommon.f90 @@ -0,0 +1,70 @@ + !> @brief Kinematic routing equation residual + !! + !! Method to calculate the kinematic-wave routing + !! residual. + !! + !< + function kinematic_residual(qa, qb, qc, qd, & + aa, ab, ac, ad, & + qsrc, length, weight, delt, & + courant) + use KindModule, only: DP + use ConstantsModule, only: DZERO, DHALF, DONE + ! -- return variable + real(DP) :: kinematic_residual !< kinematic-wave residual + ! -- dummy variables + real(DP), intent(in) :: qa !< upstream flow at previous time + real(DP), intent(in) :: qb !< downstream flow at previous time + real(DP), intent(in) :: qc !< upstream flow + real(DP), intent(in) :: qd !< downstream flow + real(DP), intent(in) :: aa !< upstream area at previous time + real(DP), intent(in) :: ab !< downstream area at previous time + real(DP), intent(in) :: ac !< upstream area + real(DP), intent(in) :: ad !< downstream area + real(DP), intent(in) :: qsrc !< lateral flow term (L3T-1L-1) + real(DP), intent(in) :: length !< reach length + real(DP), intent(in) :: weight !< temporal weight + real(DP), intent(in) :: delt !< time step length + real(DP), intent(in) :: courant !< courant number + ! --local variables + real(DP) :: f11 + real(DP) :: f12 + + if (courant > DONE) then + f11 = (qd - qc) / length + f12 = (ac - aa) / delt + else + f11 = (weight * (qd - qc) + (DONE - weight) * (qb - qa)) / length + f12 = DHALF * ((ad - ab) + (ac - aa)) / delt + end if + kinematic_residual = f11 + f12 - qsrc + + end function kinematic_residual + + !> @brief Kinematic routing equation storage term + !! + !! Method to calculate the kinematic-wave routing + !! storage term. + !! + !< + function kinematic_storage(aa, ab, ac, ad, length, delt, courant) + use KindModule, only: DP + use ConstantsModule, only: DHALF + + real(DP) :: kinematic_storage !< kinematic-wave storage change + + real(DP), intent(in) :: aa !< upstream area at previous time + real(DP), intent(in) :: ab !< downstream area at previous time + real(DP), intent(in) :: ac !< upstream area + real(DP), intent(in) :: ad !< downstream area + real(DP), intent(in) :: length !< reach length + real(DP), intent(in) :: delt !< time step length + real(DP), intent(in) :: courant !< courant number + + if (courant > DONE) then + kinematic_storage = (aa - ac) * length / delt + else + kinematic_storage = DHALF * ((ac + ad) - (aa + ab)) * length / delt + end if + + end function kinematic_storage diff --git a/src/meson.build b/src/meson.build index eb6c23ee1e3..90e47691f2f 100644 --- a/src/meson.build +++ b/src/meson.build @@ -194,6 +194,7 @@ modflow_sources = files( 'Model' / 'GroundWaterFlow' / 'gwf-vsc.f90', 'Model' / 'GroundWaterFlow' / 'gwf-wel.f90', 'Model' / 'GroundWaterFlow' / 'submodules' / 'gwf-sfr-steady.f90', + 'Model' / 'GroundWaterFlow' / 'submodules' / 'gwf-sfr-transient.f90', 'Model' / 'GroundWaterFlow' / 'submodules' / 'gwf-sfr-constant.f90', 'Model' / 'GroundWaterTransport' / 'gwt.f90', 'Model' / 'GroundWaterTransport' / 'gwt-cnc.f90', @@ -234,6 +235,7 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'GwfConductanceUtils.f90', 'Model' / 'ModelUtilities' / 'GwfMvrPeriodData.f90', 'Model' / 'ModelUtilities' / 'GwfNpfOptions.f90', + 'Model' / 'ModelUtilities' / 'GwfSfrCommon.f90', 'Model' / 'ModelUtilities' / 'GwfStorageUtils.f90', 'Model' / 'ModelUtilities' / 'GwfVscInputData.f90', 'Model' / 'ModelUtilities' / 'GwtDspOptions.f90', diff --git a/utils/zonebudget/make/makefile b/utils/zonebudget/make/makefile index 7cc155da31c..4cb973e9759 100644 --- a/utils/zonebudget/make/makefile +++ b/utils/zonebudget/make/makefile @@ -1,4 +1,4 @@ -# makefile created by pymake (version 1.2.11.dev0) for the 'zbud6' executable. +# makefile created by pymake (version 1.2.10.dev0) for the 'zbud6' executable. include ./makedefaults