From 24a56f2f2399e6eb0f83825ae57a5fe403f158e7 Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D." Date: Mon, 27 Jan 2025 13:40:09 -0600 Subject: [PATCH 1/4] feat(swf-evp): add new evaporation package for surface water models --- doc/mf6io/chf/chf.tex | 4 + doc/mf6io/chf/evp.tex | 47 ++ doc/mf6io/chf/namefile.tex | 2 + doc/mf6io/mf6ivar/dfn/chf-evp.dfn | 211 ++++++ doc/mf6io/mf6ivar/dfn/olf-evp.dfn | 211 ++++++ doc/mf6io/mf6ivar/dfn/swf-evp.dfn | 211 ++++++ .../mf6ivar/examples/chf-evp-example.dat | 13 + .../mf6ivar/examples/olf-evp-example.dat | 13 + doc/mf6io/olf/evp.tex | 47 ++ doc/mf6io/olf/namefile.tex | 2 + doc/mf6io/olf/olf.tex | 6 +- msvs/mf6core.vfproj | 4 + src/Idm/chf-evpidm.f90 | 437 +++++++++++++ src/Idm/olf-evpidm.f90 | 437 +++++++++++++ src/Idm/selector/IdmChfDfnSelector.f90 | 13 + src/Idm/selector/IdmOlfDfnSelector.f90 | 13 + src/Idm/selector/IdmSwfDfnSelector.f90 | 13 + src/Idm/swf-evpidm.f90 | 437 +++++++++++++ src/Model/ChannelFlow/chf.f90 | 3 +- src/Model/ModelUtilities/SwfCxsUtils.f90 | 12 +- src/Model/OverlandFlow/olf.f90 | 3 +- src/Model/SurfaceWaterFlow/swf-cxs.f90 | 2 +- src/Model/SurfaceWaterFlow/swf-evp.f90 | 600 ++++++++++++++++++ src/Model/SurfaceWaterFlow/swf-pcp.f90 | 3 + src/Model/SurfaceWaterFlow/swf.f90 | 6 +- src/meson.build | 4 + utils/idmloader/dfns.txt | 3 + 27 files changed, 2750 insertions(+), 7 deletions(-) create mode 100644 doc/mf6io/chf/evp.tex create mode 100644 doc/mf6io/mf6ivar/dfn/chf-evp.dfn create mode 100644 doc/mf6io/mf6ivar/dfn/olf-evp.dfn create mode 100644 doc/mf6io/mf6ivar/dfn/swf-evp.dfn create mode 100644 doc/mf6io/mf6ivar/examples/chf-evp-example.dat create mode 100644 doc/mf6io/mf6ivar/examples/olf-evp-example.dat create mode 100644 doc/mf6io/olf/evp.tex create mode 100644 src/Idm/chf-evpidm.f90 create mode 100644 src/Idm/olf-evpidm.f90 create mode 100644 src/Idm/swf-evpidm.f90 create mode 100644 src/Model/SurfaceWaterFlow/swf-evp.f90 diff --git a/doc/mf6io/chf/chf.tex b/doc/mf6io/chf/chf.tex index 5279aeeb9c6..c9b50a82ab1 100644 --- a/doc/mf6io/chf/chf.tex +++ b/doc/mf6io/chf/chf.tex @@ -50,6 +50,10 @@ \subsection{Inflow (FLW) Package} \subsection{Precipitation (PCP) Package} \input{chf/pcp} +\newpage +\subsection{Evaporation (EVP) Package} +\input{chf/evp} + \newpage \subsection{Zero-Depth Gradient (ZDG) Package} \input{chf/zdg} diff --git a/doc/mf6io/chf/evp.tex b/doc/mf6io/chf/evp.tex new file mode 100644 index 00000000000..6be1cdb883b --- /dev/null +++ b/doc/mf6io/chf/evp.tex @@ -0,0 +1,47 @@ +Input to the Evaporation (EVP) Package is read from the file that has type ``EVP6'' in the Name File. Any number of EVP Packages can be specified for a single surface water flow model. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\vspace{5mm} + +\noindent \textit{FOR EACH SIMULATION} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/chf-evp-options.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/chf-evp-dimensions.dat} +\vspace{5mm} +\noindent \textit{FOR ANY STRESS PERIOD} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/chf-evp-period.dat} +\packageperioddescription + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/chf-evp-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/chf-evp-example.dat} + +%\vspace{5mm} +%\subsubsection{Available observation types} +%Well Package observations include the simulated well rates (\texttt{wel}), the well discharge that is available for the MVR package (\texttt{to-mvr}), and the reduction in the specified \texttt{q} when the \texttt{AUTO\_FLOW\_REDUCE} option is enabled. The data required for each WEL Package observation type is defined in table~\ref{table:gwf-welobstype}. The sum of \texttt{wel} and \texttt{to-mvr} is equal to the simulated well discharge rate, which may be less than the specified \texttt{q} if the \texttt{AUTO\_FLOW\_REDUCE} option is enabled. The \texttt{DNODATA} value is returned if the \texttt{wel-reduction} observation is specified but the \texttt{AUTO\_FLOW\_REDUCE} option is not enabled. Negative and positive values for an observation represent a loss from and gain to the GWF model, respectively. + +%\begin{longtable}{p{2cm} p{2.75cm} p{2cm} p{1.25cm} p{7cm}} +%\caption{Available WEL Package observation types} \tabularnewline + +%\hline +%\hline +%\textbf{Stress Package} & \textbf{Observation type} & \textbf{ID} & \textbf{ID2} & \textbf{Description} \\ +%\hline +%\endhead + +%\hline +%\endfoot + +%\input{../Common/gwf-welobs.tex} +%\label{table:gwf-welobstype} +%\end{longtable} + +%\vspace{5mm} +%\subsubsection{Example Observation Input File} +%\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwf-wel-example-obs.dat} diff --git a/doc/mf6io/chf/namefile.tex b/doc/mf6io/chf/namefile.tex index bda985534de..65e1df2168c 100644 --- a/doc/mf6io/chf/namefile.tex +++ b/doc/mf6io/chf/namefile.tex @@ -30,6 +30,8 @@ \subsubsection{Explanation of Variables} STO6 & Storage Package \\ CHD6 & Constant Head Package & * \\ FLW6 & Inflow Package & * \\ +PCP6 & Precipitation Package & * \\ +EVP6 & Evaporation Package & * \\ ZDG6 & Zero-Depth Gradient Package & * \\ CDB6 & Critical Depth Boundary Package & * \\ OBS6 & Observations Option \\ diff --git a/doc/mf6io/mf6ivar/dfn/chf-evp.dfn b/doc/mf6io/mf6ivar/dfn/chf-evp.dfn new file mode 100644 index 00000000000..e2c7321ac41 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/chf-evp.dfn @@ -0,0 +1,211 @@ +# --------------------- swf evp options --------------------- +# flopy multi-package + +block options +name auxiliary +type string +shape (naux) +reader urword +optional true +longname keyword to specify aux variables +description REPLACE auxnames {'{#1}': 'Surface Water Flow'} + +block options +name auxmultname +type string +shape +reader urword +optional true +longname name of auxiliary variable for multiplier +description REPLACE auxmultname {'{#1}': 'evaporation'} + +block options +name boundnames +type keyword +shape +reader urword +optional true +longname +description REPLACE boundnames {'{#1}': 'evaporation'} + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'evaporation'} +mf6internal iprpak + +block options +name print_flows +type keyword +reader urword +optional true +longname print evaporation rates to listing file +description REPLACE print_flows {'{#1}': 'evaporation'} +mf6internal iprflow + +block options +name save_flows +type keyword +reader urword +optional true +longname save evaporation to budget file +description REPLACE save_flows {'{#1}': 'evaporation'} +mf6internal ipakcb + +block options +name ts_filerecord +type record ts6 filein ts6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name ts6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname head keyword +description keyword to specify that record corresponds to a time-series file. + +block options +name filein +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname file keyword +description keyword to specify that an input filename is expected next. + +block options +name ts6_filename +type string +preserve_case true +in_record true +reader urword +optional false +tagged false +longname file name of time series information +description REPLACE timeseriesfile {} + +block options +name obs_filerecord +type record obs6 filein obs6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name obs6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname obs keyword +description keyword to specify that record corresponds to an observations file. + +block options +name obs6_filename +type string +preserve_case true +in_record true +tagged false +reader urword +optional false +longname obs6 input filename +description REPLACE obs6_filename {'{#1}': 'Evaporation'} + +# --------------------- swf evp dimensions --------------------- + +block dimensions +name maxbound +type integer +reader urword +optional false +longname maximum number of evaporation cells +description REPLACE maxbound {'{#1}': 'evaporation cells'} + + +# --------------------- swf evp period --------------------- + +block period +name iper +type integer +block_variable True +in_record true +tagged false +shape +valid +reader urword +optional false +longname stress period number +description REPLACE iper {} + +block period +name stress_period_data +type recarray cellid evaporation aux boundname +shape (maxbound) +reader urword +longname +description +mf6internal spd + +block period +name cellid +type integer +shape (ncelldim) +tagged false +in_record true +reader urword +longname cell identifier +description REPLACE cellid {} + +block period +name evaporation +type double precision +shape +tagged false +in_record true +reader urword +time_series true +longname evaporation rate +description is the evaporation flux rate ($LT^{-1}$). This rate is multiplied inside the program by the water surface area of the cell to calculate the volumetric evaporation rate. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name aux +type double precision +in_record true +tagged false +shape (naux) +reader urword +optional true +time_series true +longname auxiliary variables +description REPLACE aux {'{#1}': 'evaporation'} +mf6internal auxvar + +block period +name boundname +type string +shape +tagged false +in_record true +reader urword +optional true +longname evaporation name +description REPLACE boundname {'{#1}': 'evaporation'} diff --git a/doc/mf6io/mf6ivar/dfn/olf-evp.dfn b/doc/mf6io/mf6ivar/dfn/olf-evp.dfn new file mode 100644 index 00000000000..e2c7321ac41 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/olf-evp.dfn @@ -0,0 +1,211 @@ +# --------------------- swf evp options --------------------- +# flopy multi-package + +block options +name auxiliary +type string +shape (naux) +reader urword +optional true +longname keyword to specify aux variables +description REPLACE auxnames {'{#1}': 'Surface Water Flow'} + +block options +name auxmultname +type string +shape +reader urword +optional true +longname name of auxiliary variable for multiplier +description REPLACE auxmultname {'{#1}': 'evaporation'} + +block options +name boundnames +type keyword +shape +reader urword +optional true +longname +description REPLACE boundnames {'{#1}': 'evaporation'} + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'evaporation'} +mf6internal iprpak + +block options +name print_flows +type keyword +reader urword +optional true +longname print evaporation rates to listing file +description REPLACE print_flows {'{#1}': 'evaporation'} +mf6internal iprflow + +block options +name save_flows +type keyword +reader urword +optional true +longname save evaporation to budget file +description REPLACE save_flows {'{#1}': 'evaporation'} +mf6internal ipakcb + +block options +name ts_filerecord +type record ts6 filein ts6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name ts6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname head keyword +description keyword to specify that record corresponds to a time-series file. + +block options +name filein +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname file keyword +description keyword to specify that an input filename is expected next. + +block options +name ts6_filename +type string +preserve_case true +in_record true +reader urword +optional false +tagged false +longname file name of time series information +description REPLACE timeseriesfile {} + +block options +name obs_filerecord +type record obs6 filein obs6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name obs6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname obs keyword +description keyword to specify that record corresponds to an observations file. + +block options +name obs6_filename +type string +preserve_case true +in_record true +tagged false +reader urword +optional false +longname obs6 input filename +description REPLACE obs6_filename {'{#1}': 'Evaporation'} + +# --------------------- swf evp dimensions --------------------- + +block dimensions +name maxbound +type integer +reader urword +optional false +longname maximum number of evaporation cells +description REPLACE maxbound {'{#1}': 'evaporation cells'} + + +# --------------------- swf evp period --------------------- + +block period +name iper +type integer +block_variable True +in_record true +tagged false +shape +valid +reader urword +optional false +longname stress period number +description REPLACE iper {} + +block period +name stress_period_data +type recarray cellid evaporation aux boundname +shape (maxbound) +reader urword +longname +description +mf6internal spd + +block period +name cellid +type integer +shape (ncelldim) +tagged false +in_record true +reader urword +longname cell identifier +description REPLACE cellid {} + +block period +name evaporation +type double precision +shape +tagged false +in_record true +reader urword +time_series true +longname evaporation rate +description is the evaporation flux rate ($LT^{-1}$). This rate is multiplied inside the program by the water surface area of the cell to calculate the volumetric evaporation rate. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name aux +type double precision +in_record true +tagged false +shape (naux) +reader urword +optional true +time_series true +longname auxiliary variables +description REPLACE aux {'{#1}': 'evaporation'} +mf6internal auxvar + +block period +name boundname +type string +shape +tagged false +in_record true +reader urword +optional true +longname evaporation name +description REPLACE boundname {'{#1}': 'evaporation'} diff --git a/doc/mf6io/mf6ivar/dfn/swf-evp.dfn b/doc/mf6io/mf6ivar/dfn/swf-evp.dfn new file mode 100644 index 00000000000..e2c7321ac41 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/swf-evp.dfn @@ -0,0 +1,211 @@ +# --------------------- swf evp options --------------------- +# flopy multi-package + +block options +name auxiliary +type string +shape (naux) +reader urword +optional true +longname keyword to specify aux variables +description REPLACE auxnames {'{#1}': 'Surface Water Flow'} + +block options +name auxmultname +type string +shape +reader urword +optional true +longname name of auxiliary variable for multiplier +description REPLACE auxmultname {'{#1}': 'evaporation'} + +block options +name boundnames +type keyword +shape +reader urword +optional true +longname +description REPLACE boundnames {'{#1}': 'evaporation'} + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'evaporation'} +mf6internal iprpak + +block options +name print_flows +type keyword +reader urword +optional true +longname print evaporation rates to listing file +description REPLACE print_flows {'{#1}': 'evaporation'} +mf6internal iprflow + +block options +name save_flows +type keyword +reader urword +optional true +longname save evaporation to budget file +description REPLACE save_flows {'{#1}': 'evaporation'} +mf6internal ipakcb + +block options +name ts_filerecord +type record ts6 filein ts6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name ts6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname head keyword +description keyword to specify that record corresponds to a time-series file. + +block options +name filein +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname file keyword +description keyword to specify that an input filename is expected next. + +block options +name ts6_filename +type string +preserve_case true +in_record true +reader urword +optional false +tagged false +longname file name of time series information +description REPLACE timeseriesfile {} + +block options +name obs_filerecord +type record obs6 filein obs6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name obs6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname obs keyword +description keyword to specify that record corresponds to an observations file. + +block options +name obs6_filename +type string +preserve_case true +in_record true +tagged false +reader urword +optional false +longname obs6 input filename +description REPLACE obs6_filename {'{#1}': 'Evaporation'} + +# --------------------- swf evp dimensions --------------------- + +block dimensions +name maxbound +type integer +reader urword +optional false +longname maximum number of evaporation cells +description REPLACE maxbound {'{#1}': 'evaporation cells'} + + +# --------------------- swf evp period --------------------- + +block period +name iper +type integer +block_variable True +in_record true +tagged false +shape +valid +reader urword +optional false +longname stress period number +description REPLACE iper {} + +block period +name stress_period_data +type recarray cellid evaporation aux boundname +shape (maxbound) +reader urword +longname +description +mf6internal spd + +block period +name cellid +type integer +shape (ncelldim) +tagged false +in_record true +reader urword +longname cell identifier +description REPLACE cellid {} + +block period +name evaporation +type double precision +shape +tagged false +in_record true +reader urword +time_series true +longname evaporation rate +description is the evaporation flux rate ($LT^{-1}$). This rate is multiplied inside the program by the water surface area of the cell to calculate the volumetric evaporation rate. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name aux +type double precision +in_record true +tagged false +shape (naux) +reader urword +optional true +time_series true +longname auxiliary variables +description REPLACE aux {'{#1}': 'evaporation'} +mf6internal auxvar + +block period +name boundname +type string +shape +tagged false +in_record true +reader urword +optional true +longname evaporation name +description REPLACE boundname {'{#1}': 'evaporation'} diff --git a/doc/mf6io/mf6ivar/examples/chf-evp-example.dat b/doc/mf6io/mf6ivar/examples/chf-evp-example.dat new file mode 100644 index 00000000000..d6c462d1c26 --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/chf-evp-example.dat @@ -0,0 +1,13 @@ +BEGIN OPTIONS + PRINT_INPUT + PRINT_FLOWS +END OPTIONS + +BEGIN dimensions + MAXBOUND 2 +END dimensions + +BEGIN PERIOD 1 + 1 0.005 + 2 0.007 +END PERIOD 1 diff --git a/doc/mf6io/mf6ivar/examples/olf-evp-example.dat b/doc/mf6io/mf6ivar/examples/olf-evp-example.dat new file mode 100644 index 00000000000..d6c462d1c26 --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/olf-evp-example.dat @@ -0,0 +1,13 @@ +BEGIN OPTIONS + PRINT_INPUT + PRINT_FLOWS +END OPTIONS + +BEGIN dimensions + MAXBOUND 2 +END dimensions + +BEGIN PERIOD 1 + 1 0.005 + 2 0.007 +END PERIOD 1 diff --git a/doc/mf6io/olf/evp.tex b/doc/mf6io/olf/evp.tex new file mode 100644 index 00000000000..515bae0425a --- /dev/null +++ b/doc/mf6io/olf/evp.tex @@ -0,0 +1,47 @@ +Input to the Evaporation (EVP) Package is read from the file that has type ``EVP6'' in the Name File. Any number of EVP Packages can be specified for a single surface water flow model. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\vspace{5mm} + +\noindent \textit{FOR EACH SIMULATION} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/olf-evp-options.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/olf-evp-dimensions.dat} +\vspace{5mm} +\noindent \textit{FOR ANY STRESS PERIOD} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/olf-evp-period.dat} +\packageperioddescription + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/olf-evp-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/olf-evp-example.dat} + +%\vspace{5mm} +%\subsubsection{Available observation types} +%Well Package observations include the simulated well rates (\texttt{wel}), the well discharge that is available for the MVR package (\texttt{to-mvr}), and the reduction in the specified \texttt{q} when the \texttt{AUTO\_FLOW\_REDUCE} option is enabled. The data required for each WEL Package observation type is defined in table~\ref{table:gwf-welobstype}. The sum of \texttt{wel} and \texttt{to-mvr} is equal to the simulated well discharge rate, which may be less than the specified \texttt{q} if the \texttt{AUTO\_FLOW\_REDUCE} option is enabled. The \texttt{DNODATA} value is returned if the \texttt{wel-reduction} observation is specified but the \texttt{AUTO\_FLOW\_REDUCE} option is not enabled. Negative and positive values for an observation represent a loss from and gain to the GWF model, respectively. + +%\begin{longtable}{p{2cm} p{2.75cm} p{2cm} p{1.25cm} p{7cm}} +%\caption{Available WEL Package observation types} \tabularnewline + +%\hline +%\hline +%\textbf{Stress Package} & \textbf{Observation type} & \textbf{ID} & \textbf{ID2} & \textbf{Description} \\ +%\hline +%\endhead + +%\hline +%\endfoot + +%\input{../Common/gwf-welobs.tex} +%\label{table:gwf-welobstype} +%\end{longtable} + +%\vspace{5mm} +%\subsubsection{Example Observation Input File} +%\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwf-wel-example-obs.dat} diff --git a/doc/mf6io/olf/namefile.tex b/doc/mf6io/olf/namefile.tex index fe1e94e7707..b83d68c2232 100644 --- a/doc/mf6io/olf/namefile.tex +++ b/doc/mf6io/olf/namefile.tex @@ -30,6 +30,8 @@ \subsubsection{Explanation of Variables} STO6 & Storage Package \\ CHD6 & Constant Head Package & * \\ FLW6 & Inflow Package & * \\ +PCP6 & Precipitation Package & * \\ +EVP6 & Evaporation Package & * \\ ZDG6 & Zero-Depth Gradient Package & * \\ CDB6 & Critical Depth Boundary Package & * \\ OBS6 & Observations Option \\ diff --git a/doc/mf6io/olf/olf.tex b/doc/mf6io/olf/olf.tex index 74e75f288cc..a164a937794 100644 --- a/doc/mf6io/olf/olf.tex +++ b/doc/mf6io/olf/olf.tex @@ -47,9 +47,13 @@ \subsection{Inflow (FLW) Package} \input{olf/flw} \newpage -\subsection{Inflow (PCP) Package} +\subsection{Precipitation (PCP) Package} \input{olf/pcp} +\newpage +\subsection{Evaporation (EVP) Package} +\input{olf/evp} + \newpage \subsection{Zero-Depth Gradient (ZDG) Package} \input{olf/zdg} diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index c8b4dfcc745..173962bf69b 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -145,6 +145,7 @@ + @@ -197,6 +198,7 @@ + @@ -216,6 +218,7 @@ + @@ -354,6 +357,7 @@ + diff --git a/src/Idm/chf-evpidm.f90 b/src/Idm/chf-evpidm.f90 new file mode 100644 index 00000000000..ea7f5836092 --- /dev/null +++ b/src/Idm/chf-evpidm.f90 @@ -0,0 +1,437 @@ +! ** Do Not Modify! MODFLOW 6 system generated file. ** +module ChfEvpInputModule + use ConstantsModule, only: LENVARNAME + use InputDefinitionModule, only: InputParamDefinitionType, & + InputBlockDefinitionType + private + public chf_evp_param_definitions + public chf_evp_aggregate_definitions + public chf_evp_block_definitions + public ChfEvpParamFoundType + public chf_evp_multi_package + public chf_evp_subpackages + + type ChfEvpParamFoundType + logical :: auxiliary = .false. + logical :: auxmultname = .false. + logical :: boundnames = .false. + logical :: iprpak = .false. + logical :: iprflow = .false. + logical :: ipakcb = .false. + logical :: ts_filerecord = .false. + logical :: ts6 = .false. + logical :: filein = .false. + logical :: ts6_filename = .false. + logical :: obs_filerecord = .false. + logical :: obs6 = .false. + logical :: obs6_filename = .false. + logical :: maxbound = .false. + logical :: cellid = .false. + logical :: evaporation = .false. + logical :: auxvar = .false. + logical :: boundname = .false. + end type ChfEvpParamFoundType + + logical :: chf_evp_multi_package = .true. + + character(len=16), parameter :: & + chf_evp_subpackages(*) = & + [ & + ' ' & + ] + + type(InputParamDefinitionType), parameter :: & + chfevp_auxiliary = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'AUXILIARY', & ! tag name + 'AUXILIARY', & ! fortran variable + 'STRING', & ! type + 'NAUX', & ! shape + 'keyword to specify aux variables', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_auxmultname = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'AUXMULTNAME', & ! tag name + 'AUXMULTNAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'name of auxiliary variable for multiplier', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_boundnames = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'BOUNDNAMES', & ! tag name + 'BOUNDNAMES', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_iprpak = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'PRINT_INPUT', & ! tag name + 'IPRPAK', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'print input to listing file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_iprflow = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'PRINT_FLOWS', & ! tag name + 'IPRFLOW', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'print evaporation rates to listing file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_ipakcb = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'SAVE_FLOWS', & ! tag name + 'IPAKCB', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'save evaporation to budget file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_ts_filerecord = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS_FILERECORD', & ! tag name + 'TS_FILERECORD', & ! fortran variable + 'RECORD TS6 FILEIN TS6_FILENAME', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_ts6 = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS6', & ! tag name + 'TS6', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'head keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_filein = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'FILEIN', & ! tag name + 'FILEIN', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'file keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_ts6_filename = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS6_FILENAME', & ! tag name + 'TS6_FILENAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'file name of time series information', & ! longname + .true., & ! required + .true., & ! multi-record + .true., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_obs_filerecord = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS_FILERECORD', & ! tag name + 'OBS_FILERECORD', & ! fortran variable + 'RECORD OBS6 FILEIN OBS6_FILENAME', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_obs6 = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6', & ! tag name + 'OBS6', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'obs keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_obs6_filename = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6_FILENAME', & ! tag name + 'OBS6_FILENAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'obs6 input filename', & ! longname + .true., & ! required + .true., & ! multi-record + .true., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_maxbound = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'DIMENSIONS', & ! block + 'MAXBOUND', & ! tag name + 'MAXBOUND', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + 'maximum number of evaporation cells', & ! longname + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_cellid = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'CELLID', & ! tag name + 'CELLID', & ! fortran variable + 'INTEGER1D', & ! type + 'NCELLDIM', & ! shape + 'cell identifier', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_evaporation = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'EVAPORATION', & ! tag name + 'EVAPORATION', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + 'evaporation rate', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .true. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_auxvar = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'AUX', & ! tag name + 'AUXVAR', & ! fortran variable + 'DOUBLE1D', & ! type + 'NAUX', & ! shape + 'auxiliary variables', & ! longname + .false., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .true. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chfevp_boundname = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'BOUNDNAME', & ! tag name + 'BOUNDNAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'evaporation name', & ! longname + .false., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chf_evp_param_definitions(*) = & + [ & + chfevp_auxiliary, & + chfevp_auxmultname, & + chfevp_boundnames, & + chfevp_iprpak, & + chfevp_iprflow, & + chfevp_ipakcb, & + chfevp_ts_filerecord, & + chfevp_ts6, & + chfevp_filein, & + chfevp_ts6_filename, & + chfevp_obs_filerecord, & + chfevp_obs6, & + chfevp_obs6_filename, & + chfevp_maxbound, & + chfevp_cellid, & + chfevp_evaporation, & + chfevp_auxvar, & + chfevp_boundname & + ] + + type(InputParamDefinitionType), parameter :: & + chfevp_spd = InputParamDefinitionType & + ( & + 'CHF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'STRESS_PERIOD_DATA', & ! tag name + 'SPD', & ! fortran variable + 'RECARRAY CELLID EVAPORATION AUX BOUNDNAME', & ! type + 'MAXBOUND', & ! shape + '', & ! longname + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + chf_evp_aggregate_definitions(*) = & + [ & + chfevp_spd & + ] + + type(InputBlockDefinitionType), parameter :: & + chf_evp_block_definitions(*) = & + [ & + InputBlockDefinitionType( & + 'OPTIONS', & ! blockname + .false., & ! required + .false., & ! aggregate + .false. & ! block_variable + ), & + InputBlockDefinitionType( & + 'DIMENSIONS', & ! blockname + .true., & ! required + .false., & ! aggregate + .false. & ! block_variable + ), & + InputBlockDefinitionType( & + 'PERIOD', & ! blockname + .true., & ! required + .true., & ! aggregate + .true. & ! block_variable + ) & + ] + +end module ChfEvpInputModule diff --git a/src/Idm/olf-evpidm.f90 b/src/Idm/olf-evpidm.f90 new file mode 100644 index 00000000000..e6d8369522d --- /dev/null +++ b/src/Idm/olf-evpidm.f90 @@ -0,0 +1,437 @@ +! ** Do Not Modify! MODFLOW 6 system generated file. ** +module OlfEvpInputModule + use ConstantsModule, only: LENVARNAME + use InputDefinitionModule, only: InputParamDefinitionType, & + InputBlockDefinitionType + private + public olf_evp_param_definitions + public olf_evp_aggregate_definitions + public olf_evp_block_definitions + public OlfEvpParamFoundType + public olf_evp_multi_package + public olf_evp_subpackages + + type OlfEvpParamFoundType + logical :: auxiliary = .false. + logical :: auxmultname = .false. + logical :: boundnames = .false. + logical :: iprpak = .false. + logical :: iprflow = .false. + logical :: ipakcb = .false. + logical :: ts_filerecord = .false. + logical :: ts6 = .false. + logical :: filein = .false. + logical :: ts6_filename = .false. + logical :: obs_filerecord = .false. + logical :: obs6 = .false. + logical :: obs6_filename = .false. + logical :: maxbound = .false. + logical :: cellid = .false. + logical :: evaporation = .false. + logical :: auxvar = .false. + logical :: boundname = .false. + end type OlfEvpParamFoundType + + logical :: olf_evp_multi_package = .true. + + character(len=16), parameter :: & + olf_evp_subpackages(*) = & + [ & + ' ' & + ] + + type(InputParamDefinitionType), parameter :: & + olfevp_auxiliary = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'AUXILIARY', & ! tag name + 'AUXILIARY', & ! fortran variable + 'STRING', & ! type + 'NAUX', & ! shape + 'keyword to specify aux variables', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_auxmultname = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'AUXMULTNAME', & ! tag name + 'AUXMULTNAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'name of auxiliary variable for multiplier', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_boundnames = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'BOUNDNAMES', & ! tag name + 'BOUNDNAMES', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_iprpak = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'PRINT_INPUT', & ! tag name + 'IPRPAK', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'print input to listing file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_iprflow = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'PRINT_FLOWS', & ! tag name + 'IPRFLOW', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'print evaporation rates to listing file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_ipakcb = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'SAVE_FLOWS', & ! tag name + 'IPAKCB', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'save evaporation to budget file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_ts_filerecord = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS_FILERECORD', & ! tag name + 'TS_FILERECORD', & ! fortran variable + 'RECORD TS6 FILEIN TS6_FILENAME', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_ts6 = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS6', & ! tag name + 'TS6', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'head keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_filein = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'FILEIN', & ! tag name + 'FILEIN', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'file keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_ts6_filename = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS6_FILENAME', & ! tag name + 'TS6_FILENAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'file name of time series information', & ! longname + .true., & ! required + .true., & ! multi-record + .true., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_obs_filerecord = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS_FILERECORD', & ! tag name + 'OBS_FILERECORD', & ! fortran variable + 'RECORD OBS6 FILEIN OBS6_FILENAME', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_obs6 = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6', & ! tag name + 'OBS6', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'obs keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_obs6_filename = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6_FILENAME', & ! tag name + 'OBS6_FILENAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'obs6 input filename', & ! longname + .true., & ! required + .true., & ! multi-record + .true., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_maxbound = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'DIMENSIONS', & ! block + 'MAXBOUND', & ! tag name + 'MAXBOUND', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + 'maximum number of evaporation cells', & ! longname + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_cellid = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'CELLID', & ! tag name + 'CELLID', & ! fortran variable + 'INTEGER1D', & ! type + 'NCELLDIM', & ! shape + 'cell identifier', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_evaporation = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'EVAPORATION', & ! tag name + 'EVAPORATION', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + 'evaporation rate', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .true. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_auxvar = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'AUX', & ! tag name + 'AUXVAR', & ! fortran variable + 'DOUBLE1D', & ! type + 'NAUX', & ! shape + 'auxiliary variables', & ! longname + .false., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .true. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olfevp_boundname = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'BOUNDNAME', & ! tag name + 'BOUNDNAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'evaporation name', & ! longname + .false., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olf_evp_param_definitions(*) = & + [ & + olfevp_auxiliary, & + olfevp_auxmultname, & + olfevp_boundnames, & + olfevp_iprpak, & + olfevp_iprflow, & + olfevp_ipakcb, & + olfevp_ts_filerecord, & + olfevp_ts6, & + olfevp_filein, & + olfevp_ts6_filename, & + olfevp_obs_filerecord, & + olfevp_obs6, & + olfevp_obs6_filename, & + olfevp_maxbound, & + olfevp_cellid, & + olfevp_evaporation, & + olfevp_auxvar, & + olfevp_boundname & + ] + + type(InputParamDefinitionType), parameter :: & + olfevp_spd = InputParamDefinitionType & + ( & + 'OLF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'STRESS_PERIOD_DATA', & ! tag name + 'SPD', & ! fortran variable + 'RECARRAY CELLID EVAPORATION AUX BOUNDNAME', & ! type + 'MAXBOUND', & ! shape + '', & ! longname + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + olf_evp_aggregate_definitions(*) = & + [ & + olfevp_spd & + ] + + type(InputBlockDefinitionType), parameter :: & + olf_evp_block_definitions(*) = & + [ & + InputBlockDefinitionType( & + 'OPTIONS', & ! blockname + .false., & ! required + .false., & ! aggregate + .false. & ! block_variable + ), & + InputBlockDefinitionType( & + 'DIMENSIONS', & ! blockname + .true., & ! required + .false., & ! aggregate + .false. & ! block_variable + ), & + InputBlockDefinitionType( & + 'PERIOD', & ! blockname + .true., & ! required + .true., & ! aggregate + .true. & ! block_variable + ) & + ] + +end module OlfEvpInputModule diff --git a/src/Idm/selector/IdmChfDfnSelector.f90 b/src/Idm/selector/IdmChfDfnSelector.f90 index e44bd63aa0e..a55f9e3f2bc 100644 --- a/src/Idm/selector/IdmChfDfnSelector.f90 +++ b/src/Idm/selector/IdmChfDfnSelector.f90 @@ -14,6 +14,7 @@ module IdmChfDfnSelectorModule use ChfChdInputModule use ChfFlwInputModule use ChfPcpInputModule + use ChfEvpInputModule use ChfStoInputModule use ChfZdgInputModule @@ -69,6 +70,8 @@ function chf_param_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, chf_flw_param_definitions) case ('PCP') call set_param_pointer(input_definition, chf_pcp_param_definitions) + case ('EVP') + call set_param_pointer(input_definition, chf_evp_param_definitions) case ('STO') call set_param_pointer(input_definition, chf_sto_param_definitions) case ('ZDG') @@ -101,6 +104,8 @@ function chf_aggregate_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, chf_flw_aggregate_definitions) case ('PCP') call set_param_pointer(input_definition, chf_pcp_aggregate_definitions) + case ('EVP') + call set_param_pointer(input_definition, chf_evp_aggregate_definitions) case ('STO') call set_param_pointer(input_definition, chf_sto_aggregate_definitions) case ('ZDG') @@ -133,6 +138,8 @@ function chf_block_definitions(subcomponent) result(input_definition) call set_block_pointer(input_definition, chf_flw_block_definitions) case ('PCP') call set_block_pointer(input_definition, chf_pcp_block_definitions) + case ('EVP') + call set_block_pointer(input_definition, chf_evp_block_definitions) case ('STO') call set_block_pointer(input_definition, chf_sto_block_definitions) case ('ZDG') @@ -164,6 +171,8 @@ function chf_idm_multi_package(subcomponent) result(multi_package) multi_package = chf_flw_multi_package case ('PCP') multi_package = chf_pcp_multi_package + case ('EVP') + multi_package = chf_evp_multi_package case ('STO') multi_package = chf_sto_multi_package case ('ZDG') @@ -198,6 +207,8 @@ function chf_idm_subpackages(subcomponent) result(subpackages) call set_subpkg_pointer(subpackages, chf_flw_subpackages) case ('PCP') call set_subpkg_pointer(subpackages, chf_pcp_subpackages) + case ('EVP') + call set_subpkg_pointer(subpackages, chf_evp_subpackages) case ('STO') call set_subpkg_pointer(subpackages, chf_sto_subpackages) case ('ZDG') @@ -230,6 +241,8 @@ function chf_idm_integrated(subcomponent) result(integrated) integrated = .true. case ('PCP') integrated = .true. + case ('EVP') + integrated = .true. case ('STO') integrated = .true. case ('ZDG') diff --git a/src/Idm/selector/IdmOlfDfnSelector.f90 b/src/Idm/selector/IdmOlfDfnSelector.f90 index 33587530db8..7dee5ec1afe 100644 --- a/src/Idm/selector/IdmOlfDfnSelector.f90 +++ b/src/Idm/selector/IdmOlfDfnSelector.f90 @@ -14,6 +14,7 @@ module IdmOlfDfnSelectorModule use OlfChdInputModule use OlfFlwInputModule use OlfPcpInputModule + use OlfEvpInputModule use OlfStoInputModule use OlfZdgInputModule @@ -69,6 +70,8 @@ function olf_param_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, olf_flw_param_definitions) case ('PCP') call set_param_pointer(input_definition, olf_pcp_param_definitions) + case ('EVP') + call set_param_pointer(input_definition, olf_evp_param_definitions) case ('STO') call set_param_pointer(input_definition, olf_sto_param_definitions) case ('ZDG') @@ -101,6 +104,8 @@ function olf_aggregate_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, olf_flw_aggregate_definitions) case ('PCP') call set_param_pointer(input_definition, olf_pcp_aggregate_definitions) + case ('EVP') + call set_param_pointer(input_definition, olf_evp_aggregate_definitions) case ('STO') call set_param_pointer(input_definition, olf_sto_aggregate_definitions) case ('ZDG') @@ -133,6 +138,8 @@ function olf_block_definitions(subcomponent) result(input_definition) call set_block_pointer(input_definition, olf_flw_block_definitions) case ('PCP') call set_block_pointer(input_definition, olf_pcp_block_definitions) + case ('EVP') + call set_block_pointer(input_definition, olf_evp_block_definitions) case ('STO') call set_block_pointer(input_definition, olf_sto_block_definitions) case ('ZDG') @@ -164,6 +171,8 @@ function olf_idm_multi_package(subcomponent) result(multi_package) multi_package = olf_flw_multi_package case ('PCP') multi_package = olf_pcp_multi_package + case ('EVP') + multi_package = olf_evp_multi_package case ('STO') multi_package = olf_sto_multi_package case ('ZDG') @@ -198,6 +207,8 @@ function olf_idm_subpackages(subcomponent) result(subpackages) call set_subpkg_pointer(subpackages, olf_flw_subpackages) case ('PCP') call set_subpkg_pointer(subpackages, olf_pcp_subpackages) + case ('EVP') + call set_subpkg_pointer(subpackages, olf_evp_subpackages) case ('STO') call set_subpkg_pointer(subpackages, olf_sto_subpackages) case ('ZDG') @@ -230,6 +241,8 @@ function olf_idm_integrated(subcomponent) result(integrated) integrated = .true. case ('PCP') integrated = .true. + case ('EVP') + integrated = .true. case ('STO') integrated = .true. case ('ZDG') diff --git a/src/Idm/selector/IdmSwfDfnSelector.f90 b/src/Idm/selector/IdmSwfDfnSelector.f90 index c09d0cc160e..73c2ff1ffa2 100644 --- a/src/Idm/selector/IdmSwfDfnSelector.f90 +++ b/src/Idm/selector/IdmSwfDfnSelector.f90 @@ -16,6 +16,7 @@ module IdmSwfDfnSelectorModule use SwfChdInputModule use SwfFlwInputModule use SwfPcpInputModule + use SwfEvpInputModule use SwfStoInputModule use SwfZdgInputModule @@ -75,6 +76,8 @@ function swf_param_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, swf_flw_param_definitions) case ('PCP') call set_param_pointer(input_definition, swf_pcp_param_definitions) + case ('EVP') + call set_param_pointer(input_definition, swf_evp_param_definitions) case ('STO') call set_param_pointer(input_definition, swf_sto_param_definitions) case ('ZDG') @@ -111,6 +114,8 @@ function swf_aggregate_definitions(subcomponent) result(input_definition) call set_param_pointer(input_definition, swf_flw_aggregate_definitions) case ('PCP') call set_param_pointer(input_definition, swf_pcp_aggregate_definitions) + case ('EVP') + call set_param_pointer(input_definition, swf_evp_aggregate_definitions) case ('STO') call set_param_pointer(input_definition, swf_sto_aggregate_definitions) case ('ZDG') @@ -147,6 +152,8 @@ function swf_block_definitions(subcomponent) result(input_definition) call set_block_pointer(input_definition, swf_flw_block_definitions) case ('PCP') call set_block_pointer(input_definition, swf_pcp_block_definitions) + case ('EVP') + call set_block_pointer(input_definition, swf_evp_block_definitions) case ('STO') call set_block_pointer(input_definition, swf_sto_block_definitions) case ('ZDG') @@ -182,6 +189,8 @@ function swf_idm_multi_package(subcomponent) result(multi_package) multi_package = swf_flw_multi_package case ('PCP') multi_package = swf_pcp_multi_package + case ('EVP') + multi_package = swf_evp_multi_package case ('STO') multi_package = swf_sto_multi_package case ('ZDG') @@ -220,6 +229,8 @@ function swf_idm_subpackages(subcomponent) result(subpackages) call set_subpkg_pointer(subpackages, swf_flw_subpackages) case ('PCP') call set_subpkg_pointer(subpackages, swf_pcp_subpackages) + case ('EVP') + call set_subpkg_pointer(subpackages, swf_evp_subpackages) case ('STO') call set_subpkg_pointer(subpackages, swf_sto_subpackages) case ('ZDG') @@ -256,6 +267,8 @@ function swf_idm_integrated(subcomponent) result(integrated) integrated = .true. case ('PCP') integrated = .true. + case ('EVP') + integrated = .true. case ('STO') integrated = .true. case ('ZDG') diff --git a/src/Idm/swf-evpidm.f90 b/src/Idm/swf-evpidm.f90 new file mode 100644 index 00000000000..9ced28097e0 --- /dev/null +++ b/src/Idm/swf-evpidm.f90 @@ -0,0 +1,437 @@ +! ** Do Not Modify! MODFLOW 6 system generated file. ** +module SwfEvpInputModule + use ConstantsModule, only: LENVARNAME + use InputDefinitionModule, only: InputParamDefinitionType, & + InputBlockDefinitionType + private + public swf_evp_param_definitions + public swf_evp_aggregate_definitions + public swf_evp_block_definitions + public SwfEvpParamFoundType + public swf_evp_multi_package + public swf_evp_subpackages + + type SwfEvpParamFoundType + logical :: auxiliary = .false. + logical :: auxmultname = .false. + logical :: boundnames = .false. + logical :: iprpak = .false. + logical :: iprflow = .false. + logical :: ipakcb = .false. + logical :: ts_filerecord = .false. + logical :: ts6 = .false. + logical :: filein = .false. + logical :: ts6_filename = .false. + logical :: obs_filerecord = .false. + logical :: obs6 = .false. + logical :: obs6_filename = .false. + logical :: maxbound = .false. + logical :: cellid = .false. + logical :: evaporation = .false. + logical :: auxvar = .false. + logical :: boundname = .false. + end type SwfEvpParamFoundType + + logical :: swf_evp_multi_package = .true. + + character(len=16), parameter :: & + swf_evp_subpackages(*) = & + [ & + ' ' & + ] + + type(InputParamDefinitionType), parameter :: & + swfevp_auxiliary = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'AUXILIARY', & ! tag name + 'AUXILIARY', & ! fortran variable + 'STRING', & ! type + 'NAUX', & ! shape + 'keyword to specify aux variables', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_auxmultname = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'AUXMULTNAME', & ! tag name + 'AUXMULTNAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'name of auxiliary variable for multiplier', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_boundnames = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'BOUNDNAMES', & ! tag name + 'BOUNDNAMES', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_iprpak = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'PRINT_INPUT', & ! tag name + 'IPRPAK', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'print input to listing file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_iprflow = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'PRINT_FLOWS', & ! tag name + 'IPRFLOW', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'print evaporation rates to listing file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_ipakcb = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'SAVE_FLOWS', & ! tag name + 'IPAKCB', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'save evaporation to budget file', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_ts_filerecord = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS_FILERECORD', & ! tag name + 'TS_FILERECORD', & ! fortran variable + 'RECORD TS6 FILEIN TS6_FILENAME', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_ts6 = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS6', & ! tag name + 'TS6', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'head keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_filein = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'FILEIN', & ! tag name + 'FILEIN', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'file keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_ts6_filename = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'TS6_FILENAME', & ! tag name + 'TS6_FILENAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'file name of time series information', & ! longname + .true., & ! required + .true., & ! multi-record + .true., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_obs_filerecord = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS_FILERECORD', & ! tag name + 'OBS_FILERECORD', & ! fortran variable + 'RECORD OBS6 FILEIN OBS6_FILENAME', & ! type + '', & ! shape + '', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_obs6 = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6', & ! tag name + 'OBS6', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'obs keyword', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_obs6_filename = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6_FILENAME', & ! tag name + 'OBS6_FILENAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'obs6 input filename', & ! longname + .true., & ! required + .true., & ! multi-record + .true., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_maxbound = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'DIMENSIONS', & ! block + 'MAXBOUND', & ! tag name + 'MAXBOUND', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + 'maximum number of evaporation cells', & ! longname + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_cellid = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'CELLID', & ! tag name + 'CELLID', & ! fortran variable + 'INTEGER1D', & ! type + 'NCELLDIM', & ! shape + 'cell identifier', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_evaporation = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'EVAPORATION', & ! tag name + 'EVAPORATION', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + 'evaporation rate', & ! longname + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .true. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_auxvar = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'AUX', & ! tag name + 'AUXVAR', & ! fortran variable + 'DOUBLE1D', & ! type + 'NAUX', & ! shape + 'auxiliary variables', & ! longname + .false., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .true. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swfevp_boundname = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'BOUNDNAME', & ! tag name + 'BOUNDNAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + 'evaporation name', & ! longname + .false., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swf_evp_param_definitions(*) = & + [ & + swfevp_auxiliary, & + swfevp_auxmultname, & + swfevp_boundnames, & + swfevp_iprpak, & + swfevp_iprflow, & + swfevp_ipakcb, & + swfevp_ts_filerecord, & + swfevp_ts6, & + swfevp_filein, & + swfevp_ts6_filename, & + swfevp_obs_filerecord, & + swfevp_obs6, & + swfevp_obs6_filename, & + swfevp_maxbound, & + swfevp_cellid, & + swfevp_evaporation, & + swfevp_auxvar, & + swfevp_boundname & + ] + + type(InputParamDefinitionType), parameter :: & + swfevp_spd = InputParamDefinitionType & + ( & + 'SWF', & ! component + 'EVP', & ! subcomponent + 'PERIOD', & ! block + 'STRESS_PERIOD_DATA', & ! tag name + 'SPD', & ! fortran variable + 'RECARRAY CELLID EVAPORATION AUX BOUNDNAME', & ! type + 'MAXBOUND', & ! shape + '', & ! longname + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + + type(InputParamDefinitionType), parameter :: & + swf_evp_aggregate_definitions(*) = & + [ & + swfevp_spd & + ] + + type(InputBlockDefinitionType), parameter :: & + swf_evp_block_definitions(*) = & + [ & + InputBlockDefinitionType( & + 'OPTIONS', & ! blockname + .false., & ! required + .false., & ! aggregate + .false. & ! block_variable + ), & + InputBlockDefinitionType( & + 'DIMENSIONS', & ! blockname + .true., & ! required + .false., & ! aggregate + .false. & ! block_variable + ), & + InputBlockDefinitionType( & + 'PERIOD', & ! blockname + .true., & ! required + .true., & ! aggregate + .true. & ! block_variable + ) & + ] + +end module SwfEvpInputModule diff --git a/src/Model/ChannelFlow/chf.f90 b/src/Model/ChannelFlow/chf.f90 index 5a05cd10075..91cc0a29e3a 100644 --- a/src/Model/ChannelFlow/chf.f90 +++ b/src/Model/ChannelFlow/chf.f90 @@ -42,7 +42,8 @@ module ChfModule integer(I4B), parameter :: CHF_NMULTIPKG = 50 character(len=LENPACKAGETYPE), dimension(CHF_NMULTIPKG) :: CHF_MULTIPKG data CHF_MULTIPKG/'FLW6 ', 'CHD6 ', 'CDB6 ', 'ZDG6 ', 'PCP6 ', & ! 5 - &45*' '/ ! 50 + 'EVP6 ', ' ', ' ', ' ', ' ', & ! 10 + &40*' '/ ! 50 ! size of supported model package arrays integer(I4B), parameter :: NIUNIT_CHF = CHF_NBASEPKG + CHF_NMULTIPKG diff --git a/src/Model/ModelUtilities/SwfCxsUtils.f90 b/src/Model/ModelUtilities/SwfCxsUtils.f90 index 1770d907a05..c05feff3923 100644 --- a/src/Model/ModelUtilities/SwfCxsUtils.f90 +++ b/src/Model/ModelUtilities/SwfCxsUtils.f90 @@ -432,12 +432,20 @@ end function calc_qman_by_section !! !! @return w saturated top width !< - function get_saturated_topwidth(npts, stations) result(w) + function get_saturated_topwidth(npts, xfraction, width) result(w) ! -- dummy variables integer(I4B), intent(in) :: npts !< number of station-height data for a reach - real(DP), dimension(npts), intent(in) :: stations !< cross-section station distances (x-distance) + real(DP), dimension(npts), intent(in) :: xfraction !< cross-section station fractional distances (x-distance) + real(DP), intent(in) :: width ! -- local variables + integer(I4B) :: n real(DP) :: w + real(DP), dimension(npts) :: stations + ! + ! -- calculate station from xfractions and width + do n = 1, npts + stations(n) = xfraction(n) * width + end do ! ! -- calculate the saturated top width if (npts > 1) then diff --git a/src/Model/OverlandFlow/olf.f90 b/src/Model/OverlandFlow/olf.f90 index 3a2dc7fa4e7..ff3fe5e1ce0 100644 --- a/src/Model/OverlandFlow/olf.f90 +++ b/src/Model/OverlandFlow/olf.f90 @@ -42,7 +42,8 @@ module OlfModule integer(I4B), parameter :: OLF_NMULTIPKG = 50 character(len=LENPACKAGETYPE), dimension(OLF_NMULTIPKG) :: OLF_MULTIPKG data OLF_MULTIPKG/'FLW6 ', 'CHD6 ', 'CDB6 ', 'ZDG6 ', 'PCP6 ', & ! 5 - &45*' '/ ! 50 + 'EVP6 ', ' ', ' ', ' ', ' ', & ! 10 + &40*' '/ ! 50 ! size of supported model package arrays integer(I4B), parameter :: NIUNIT_OLF = OLF_NBASEPKG + OLF_NMULTIPKG diff --git a/src/Model/SurfaceWaterFlow/swf-cxs.f90 b/src/Model/SurfaceWaterFlow/swf-cxs.f90 index 3e3a6cd1b10..9be22f14e32 100644 --- a/src/Model/SurfaceWaterFlow/swf-cxs.f90 +++ b/src/Model/SurfaceWaterFlow/swf-cxs.f90 @@ -803,7 +803,7 @@ function get_maximum_top_width(this, idcxs, width) result(r) if (npts == 0) then r = width else - r = get_saturated_topwidth(npts, this%xfraction(i0:i1)) + r = get_saturated_topwidth(npts, this%xfraction(i0:i1), width) end if end function get_maximum_top_width diff --git a/src/Model/SurfaceWaterFlow/swf-evp.f90 b/src/Model/SurfaceWaterFlow/swf-evp.f90 new file mode 100644 index 00000000000..d16b29f909e --- /dev/null +++ b/src/Model/SurfaceWaterFlow/swf-evp.f90 @@ -0,0 +1,600 @@ +!> @brief This module contains the evaporation (EVP) package methods +!! +!! This module can be used to represent evaporation onto streams and +!! overland flow cells. +!< +module SwfEvpModule + + use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: DZERO, LENFTYPE, LENPACKAGENAME, MAXCHARLEN, & + LINELENGTH, DONE + use MemoryHelperModule, only: create_mem_path + use BndModule, only: BndType + use BndExtModule, only: BndExtType + use SimModule, only: store_error, store_error_filename, count_errors + use SimVariablesModule, only: errmsg + use ObsModule, only: DefaultObsIdProcessor + use TimeArraySeriesLinkModule, only: TimeArraySeriesLinkType + use BlockParserModule, only: BlockParserType + use CharacterStringModule, only: CharacterStringType + use SmoothingModule, only: sQSaturation + use MatrixBaseModule + use GeomUtilModule, only: get_node + use BaseDisModule, only: DisBaseType + use Disv1dModule, only: Disv1dType + use SwfDfwModule, only: SwfDfwType + use SwfCxsModule, only: SwfCxsType + + implicit none + + private + public :: evp_create + + character(len=LENFTYPE) :: ftype = 'EVP' + character(len=LENPACKAGENAME) :: text = ' EVP' + ! character(len=LENPACKAGENAME) :: texta = ' EVPA' + + type, extends(BndExtType) :: SwfEvpType + real(DP), dimension(:), pointer, contiguous :: evaporation => null() !< boundary evaporation array + logical, pointer, private :: read_as_arrays + + ! pointers to other objects + type(SwfDfwType), pointer :: dfw + type(SwfCxsType), pointer :: cxs + + contains + + procedure :: evp_allocate_scalars + procedure :: allocate_arrays => evp_allocate_arrays + procedure :: source_options => evp_source_options + procedure :: source_dimensions => evp_source_dimensions + procedure :: log_evp_options + procedure :: read_initial_attr => evp_read_initial_attr + procedure :: bnd_rp => evp_rp + procedure :: bnd_ck => evp_ck + procedure :: bnd_cf => evp_cf + procedure :: bnd_fc => evp_fc + procedure :: bnd_da => evp_da + procedure :: define_listlabel => evp_define_listlabel + procedure :: bound_value => evp_bound_value + procedure, private :: default_nodelist + procedure, private :: reach_length_pointer + procedure, private :: get_qevp + procedure, private :: get_evap_reduce_mult + ! for observations + procedure, public :: bnd_obs_supported => evp_obs_supported + procedure, public :: bnd_df_obs => evp_df_obs + + end type SwfEvpType + +contains + + !> @brief Create a Evaporation Package + !< + subroutine evp_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, & + mempath, dis, dfw, cxs) + ! dummy + class(BndType), pointer :: packobj !< pointer to default package type + integer(I4B), intent(in) :: id !< package id + integer(I4B), intent(in) :: ibcnum !< boundary condition number + integer(I4B), intent(in) :: inunit !< unit number of CDB package input file + integer(I4B), intent(in) :: iout !< unit number of model listing file + character(len=*), intent(in) :: namemodel !< model name + character(len=*), intent(in) :: pakname !< package name + character(len=*), intent(in) :: mempath !< input mempath + class(DisBaseType), pointer, intent(inout) :: dis !< the pointer to the discretization + type(SwfDfwType), pointer, intent(in) :: dfw !< the pointer to the dfw package + type(SwfCxsType), pointer, intent(in) :: cxs !< the pointer to the cxs package + ! local + type(SwfEvpType), pointer :: evpobj + + ! allocate evaporation object and scalar variables + allocate (evpobj) + packobj => evpobj + + ! create name and memory path + call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath) + packobj%text = text + + ! allocate scalars + call evpobj%evp_allocate_scalars() + + ! initialize package + call packobj%pack_initialize() + + packobj%inunit = inunit + packobj%iout = iout + packobj%id = id + packobj%ibcnum = ibcnum + packobj%ictMemPath = create_mem_path(namemodel, 'DFW') + + ! store pointer to dis + evpobj%dis => dis + + ! store pointer to dfw + evpobj%dfw => dfw + + ! store pointer to cxs + evpobj%cxs => cxs + end subroutine evp_create + + !> @brief Allocate scalar members + !< + subroutine evp_allocate_scalars(this) + ! dummy + class(SwfEvpType), intent(inout) :: this + + ! allocate base scalars + call this%BndExtType%allocate_scalars() + + ! allocate internal members + allocate (this%read_as_arrays) + + ! Set values + this%read_as_arrays = .false. + end subroutine evp_allocate_scalars + + !> @brief Allocate package arrays + !< + subroutine evp_allocate_arrays(this, nodelist, auxvar) + ! modules + use MemoryManagerModule, only: mem_setptr, mem_checkin + ! dummy + class(SwfEvpType) :: this + integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist + real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar + + ! allocate base arrays + call this%BndExtType%allocate_arrays(nodelist, auxvar) + + ! set input context pointers + call mem_setptr(this%evaporation, 'EVAPORATION', this%input_mempath) + + ! checkin input context pointers + call mem_checkin(this%evaporation, 'EVAPORATION', this%memoryPath, & + 'EVAPORATION', this%input_mempath) + end subroutine evp_allocate_arrays + + !> @brief Source options specific to EVPType + !< + subroutine evp_source_options(this) + ! modules + use MemoryManagerExtModule, only: mem_set_value + implicit none + ! dummy + class(SwfEvpType), intent(inout) :: this + ! local + logical(LGP) :: found_readasarrays = .false. + + ! source common bound options + call this%BndExtType%source_options() + + ! update defaults with idm sourced values + call mem_set_value(this%read_as_arrays, 'READASARRAYS', this%input_mempath, & + found_readasarrays) + + ! log evp params + call this%log_evp_options(found_readasarrays) + end subroutine evp_source_options + + !> @brief Log options specific to SwfEvpType + !< + subroutine log_evp_options(this, found_readasarrays) + implicit none + ! dummy + class(SwfEvpType), intent(inout) :: this + logical(LGP), intent(in) :: found_readasarrays + ! formats + character(len=*), parameter :: fmtreadasarrays = & + &"(4x, 'EVAPORATION INPUT WILL BE READ AS ARRAY(S).')" + + ! log found options + write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text)) & + //' OPTIONS' + + if (found_readasarrays) then + write (this%iout, fmtreadasarrays) + end if + + ! close logging block + write (this%iout, '(1x,a)') & + 'END OF '//trim(adjustl(this%text))//' OPTIONS' + end subroutine log_evp_options + + !> @brief Source the dimensions for this package + !< + subroutine evp_source_dimensions(this) + ! dummy + class(SwfEvpType), intent(inout) :: this + + if (this%read_as_arrays) then + + ! Set maxbound to the number of cells per layer, which is simply + ! nrow * ncol for a dis2d grid, and nodesuser for disv2d and disv1d + this%maxbound = this%dis%get_ncpl() + + ! verify dimensions were set + if (this%maxbound <= 0) then + write (errmsg, '(a)') & + 'MAXBOUND must be an integer greater than zero.' + call store_error(errmsg) + call store_error_filename(this%input_fname) + end if + + else + + ! source maxbound + call this%BndExtType%source_dimensions() + + end if + + ! Call define_listlabel to construct the list label that is written + ! when PRINT_INPUT option is used. + call this%define_listlabel() + end subroutine evp_source_dimensions + + !> @brief Part of allocate and read + !< + subroutine evp_read_initial_attr(this) + ! dummy + class(SwfEvpType), intent(inout) :: this + + if (this%read_as_arrays) then + call this%default_nodelist() + end if + end subroutine evp_read_initial_attr + + !> @brief Read and Prepare + !! + !! Read itmp and read new boundaries if itmp > 0 + !< + subroutine evp_rp(this) + ! modules + use TdisModule, only: kper + implicit none + ! dummy + class(SwfEvpType), intent(inout) :: this + + if (this%iper /= kper) return + + if (this%read_as_arrays) then + ! no need to do anything because this%evaporation points directly to + ! the input context evaporation, which is automatically updated by idm + else + call this%BndExtType%bnd_rp() + end if + + ! Write the list to iout if requested + if (this%iprpak /= 0) then + call this%write_list() + end if + end subroutine evp_rp + + !> @brief Ensure evaporation is positive + !< + subroutine evp_ck(this) + ! dummy + class(SwfEvpType), intent(inout) :: this + ! local + character(len=30) :: nodestr + integer(I4B) :: i, nr + character(len=*), parameter :: fmterr = & + &"('Specified stress ',i0, & + &' evaporation (',g0,') is less than zero for cell', a)" + + ! Ensure evaporation rates are positive + do i = 1, this%nbound + nr = this%nodelist(i) + if (nr <= 0) cycle + if (this%evaporation(i) < DZERO) then + call this%dis%noder_to_string(nr, nodestr) + write (errmsg, fmt=fmterr) i, this%evaporation(i), trim(nodestr) + call store_error(errmsg) + end if + end do + + ! write summary of package error messages + if (count_errors() > 0) then + call store_error_filename(this%input_fname) + end if + end subroutine evp_ck + + !> @brief Formulate the HCOF and RHS terms + !! + !! Skip if no evaporation. Otherwise, calculate hcof and rhs + !< + subroutine evp_cf(this) + ! dummy + class(SwfEvpType) :: this + ! local + integer(I4B) :: i + integer(I4B) :: node + real(DP) :: qevp + real(DP) :: evap + real(DP) :: rlen + real(DP), dimension(:), pointer :: reach_length + + ! Return if no evaporation + if (this%nbound == 0) return + + ! Set pointer to reach_length for 1d + reach_length => this%reach_length_pointer() + rlen = DZERO + + ! Calculate hcof and rhs for each evaporation entry + do i = 1, this%nbound + + ! Find the node number + node = this%nodelist(i) + + ! cycle if nonexistent bound + if (node <= 0) then + this%hcof(i) = DZERO + this%rhs(i) = DZERO + cycle + end if + + ! Initialize hcof + this%hcof(i) = DZERO + + ! assign evap rate in length per time and multiply by auxvar + evap = this%evaporation(i) + if (this%iauxmultcol > 0) then + evap = evap * this%auxvar(this%iauxmultcol, i) + end if + + ! get reach length for 1d channel + if (this%dis%is_1d()) then + rlen = reach_length(node) + end if + + ! calculate volumetric evaporation flow in L^3/T + qevp = this%get_qevp(node, rlen, this%xnew(node), evap) + + ! rhs contribution (positive value means extraction) + this%rhs(i) = qevp + + ! zero out contribution if cell is inactive or constant head + if (this%ibound(node) <= 0) then + this%rhs(i) = DZERO + cycle + end if + + end do + end subroutine evp_cf + + !> @brief Calculate qevp + !< + function get_qevp(this, node, rlen, stage, evaporation) result(qevp) + ! dummy + class(SwfEvpType) :: this !< this instance + integer(I4B), intent(in) :: node !< reduced node number + real(DP), intent(in) :: rlen !< length of reach + real(DP), intent(in) :: stage !< stage in reach + real(DP), intent(in) :: evaporation !< evaporation rate in length per time + ! return + real(DP) :: qevp + ! local + integer(I4B) :: idcxs + real(DP) :: depth + real(DP) :: bt + real(DP) :: area + real(DP) :: top_width + real(DP) :: width_channel + real(DP) :: dummy + real(DP) :: qmult + + ! calculate depth of water + bt = this%dis%bot(node) + depth = stage - bt + + ! Determine the water surface area + if (this%dis%is_2d()) then + ! overland flow case + area = this%dis%get_area(node) + else if (this%dis%is_1d()) then + ! channel case + idcxs = this%dfw%idcxs(node) + call this%dis%get_flow_width(node, node, 0, width_channel, dummy) + top_width = this%cxs%get_wetted_top_width(idcxs, width_channel, depth) + area = rlen * top_width + end if + + ! Reduce the evap rate as cell goes dry + qmult = this%get_evap_reduce_mult(stage, bt) + + ! calculate volumetric evaporation flow in L^3/T + qevp = evaporation * area * qmult + + end function get_qevp + + function get_evap_reduce_mult(this, stage, bottom) result(qmult) + ! dummy + class(SwfEvpType) :: this + real(DP), intent(in) :: stage + real(DP), intent(in) :: bottom + ! return + real(DP) :: qmult + ! local + integer(I4B) :: iflowred + real(DP) :: evap_depth + real(DP) :: tp + + iflowred = 1 + qmult = DONE + if (iflowred == 1) then + evap_depth = 0.1D0 + tp = bottom + evap_depth + qmult = sQSaturation(tp, bottom, stage) + end if + + end function get_evap_reduce_mult + + !> @brief Copy rhs and hcof into solution rhs and amat + !< + subroutine evp_fc(this, rhs, ia, idxglo, matrix_sln) + ! dummy + class(SwfEvpType) :: this + real(DP), dimension(:), intent(inout) :: rhs + integer(I4B), dimension(:), intent(in) :: ia + integer(I4B), dimension(:), intent(in) :: idxglo + class(MatrixBaseType), pointer :: matrix_sln + ! local + integer(I4B) :: i, n, ipos + + ! Copy package rhs and hcof into solution rhs and amat + do i = 1, this%nbound + n = this%nodelist(i) + if (n <= 0) cycle + rhs(n) = rhs(n) + this%rhs(i) + ipos = ia(n) + call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i)) + end do + end subroutine evp_fc + + !> @brief Deallocate memory + !< + subroutine evp_da(this) + ! modules + use MemoryManagerModule, only: mem_deallocate + ! dummy + class(SwfEvpType) :: this + + ! Deallocate parent package + call this%BndExtType%bnd_da() + + ! scalars + deallocate (this%read_as_arrays) + + ! arrays + call mem_deallocate(this%evaporation, 'EVAPORATION', this%memoryPath) + + ! pointers + nullify (this%dis) + nullify (this%dfw) + nullify (this%cxs) + end subroutine evp_da + + !> @brief Define the list heading that is written to iout when PRINT_INPUT + !! option is used. + !< + subroutine evp_define_listlabel(this) + ! dummy + class(SwfEvpType), intent(inout) :: this + ! + ! create the header list label + this%listlabel = trim(this%filtyp)//' NO.' + if (this%dis%ndim == 3) then + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER' + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW' + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL' + elseif (this%dis%ndim == 2) then + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER' + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D' + else + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE' + end if + write (this%listlabel, '(a, a16)') trim(this%listlabel), 'EVAPORATION' +! if(this%multindex > 0) & +! write(this%listlabel, '(a, a16)') trim(this%listlabel), 'MULTIPLIER' + if (this%inamedbound == 1) then + write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME' + end if + end subroutine evp_define_listlabel + + !> @brief Assign default nodelist when READASARRAYS is specified. + !< + subroutine default_nodelist(this) + ! dummy + class(SwfEvpType) :: this + ! local + integer(I4B) :: nodeu, noder + + ! This is only called for readasarrays, so nodelist will be the size of + ! the user grid, and will have a value of 0 for any entries where idomain + ! is not 1 + do nodeu = 1, this%maxbound + noder = this%dis%get_nodenumber(nodeu, 0) + this%nodelist(nodeu) = noder + end do + + ! Assign nbound + this%nbound = this%maxbound + + end subroutine default_nodelist + + ! Procedures related to observations + + !> @brief + !! + !! Overrides BndType%bnd_obs_supported() + !< + logical function evp_obs_supported(this) + implicit none + ! dummy + class(SwfEvpType) :: this + evp_obs_supported = .true. + end function evp_obs_supported + + !> @brief Implements bnd_df_obs + !! + !! Store observation type supported by EVP package. Overrides + !! BndType%bnd_df_obs + !< + subroutine evp_df_obs(this) + implicit none + ! dummy + class(SwfEvpType) :: this + ! local + integer(I4B) :: indx + + call this%obs%StoreObsType('evp', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => DefaultObsIdProcessor + end subroutine evp_df_obs + + !> @brief Return requested boundary value + !< + function evp_bound_value(this, col, row) result(bndval) + ! modules + use ConstantsModule, only: DZERO + ! dummy + class(SwfEvpType), intent(inout) :: this !< BndExtType object + integer(I4B), intent(in) :: col + integer(I4B), intent(in) :: row + ! result + real(DP) :: bndval + + select case (col) + case (1) + if (this%iauxmultcol > 0) then + bndval = this%evaporation(row) * this%auxvar(this%iauxmultcol, row) + else + bndval = this%evaporation(row) + end if + case default + errmsg = 'Programming error. EVP bound value requested column '& + &'outside range of ncolbnd (1).' + call store_error(errmsg) + call store_error_filename(this%input_fname) + end select + end function evp_bound_value + + function reach_length_pointer(this) result(ptr) + ! dummy + class(SwfEvpType) :: this !< this instance + ! return + real(DP), dimension(:), pointer :: ptr + ! local + class(DisBaseType), pointer :: dis + + ptr => null() + dis => this%dis + select type (dis) + type is (Disv1dType) + ptr => dis%length + end select + + end function reach_length_pointer + +end module SwfEvpModule + diff --git a/src/Model/SurfaceWaterFlow/swf-pcp.f90 b/src/Model/SurfaceWaterFlow/swf-pcp.f90 index 28c85938b96..5b37a6406c0 100644 --- a/src/Model/SurfaceWaterFlow/swf-pcp.f90 +++ b/src/Model/SurfaceWaterFlow/swf-pcp.f90 @@ -311,6 +311,7 @@ subroutine pcp_cf(this) real(DP) :: area real(DP) :: width_channel real(DP) :: top_width + real(DP) :: dummy real(DP), dimension(:), pointer :: reach_length ! Return if no precipitation @@ -342,6 +343,8 @@ subroutine pcp_cf(this) else if (this%dis%is_1d()) then ! this is for channel case idcxs = this%dfw%idcxs(node) + call this%dis%get_flow_width(node, node, 0, width_channel, & + dummy) top_width = this%cxs%get_maximum_top_width(idcxs, width_channel) area = reach_length(node) * top_width end if diff --git a/src/Model/SurfaceWaterFlow/swf.f90 b/src/Model/SurfaceWaterFlow/swf.f90 index 70223ee3f4e..ea7c741e733 100644 --- a/src/Model/SurfaceWaterFlow/swf.f90 +++ b/src/Model/SurfaceWaterFlow/swf.f90 @@ -837,6 +837,7 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, & use SwfCdbModule, only: cdb_create use SwfZdgModule, only: zdg_create use SwfPcpModule, only: pcp_create + use SwfEvpModule, only: evp_create ! dummy class(SwfModelType) :: this character(len=*), intent(in) :: filtyp @@ -871,6 +872,9 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, & case ('PCP6') call pcp_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & pakname, mempath, this%dis, this%dfw, this%cxs) + case ('EVP6') + call evp_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & + pakname, mempath, this%dis, this%dfw, this%cxs) case default write (errmsg, *) 'Invalid package type: ', filtyp call store_error(errmsg) @@ -1075,7 +1079,7 @@ subroutine create_packages(this) this%inoc = inunit case ('OBS6') this%inobs = inunit - case ('CHD6', 'FLW6', 'CDB6', 'ZDG6', 'PCP6') + case ('CHD6', 'FLW6', 'CDB6', 'ZDG6', 'PCP6', 'EVP6') call expandarray(bndpkgs) bndpkgs(size(bndpkgs)) = n case default diff --git a/src/meson.build b/src/meson.build index ff840c7d8fc..beb5a53dd0f 100644 --- a/src/meson.build +++ b/src/meson.build @@ -90,6 +90,7 @@ modflow_sources = files( 'Idm' / 'swf-flwidm.f90', 'Idm' / 'swf-chdidm.f90', 'Idm' / 'swf-cdbidm.f90', + 'Idm' / 'swf-evpidm.f90', 'Idm' / 'swf-pcpidm.f90', 'Idm' / 'swf-stoidm.f90', 'Idm' / 'swf-zdgidm.f90', @@ -103,6 +104,7 @@ modflow_sources = files( 'Idm' / 'chf-flwidm.f90', 'Idm' / 'chf-chdidm.f90', 'Idm' / 'chf-cdbidm.f90', + 'Idm' / 'chf-evpidm.f90', 'Idm' / 'chf-pcpidm.f90', 'Idm' / 'chf-stoidm.f90', 'Idm' / 'chf-zdgidm.f90', @@ -116,6 +118,7 @@ modflow_sources = files( 'Idm' / 'olf-flwidm.f90', 'Idm' / 'olf-chdidm.f90', 'Idm' / 'olf-cdbidm.f90', + 'Idm' / 'olf-evpidm.f90', 'Idm' / 'olf-pcpidm.f90', 'Idm' / 'olf-stoidm.f90', 'Idm' / 'olf-zdgidm.f90', @@ -217,6 +220,7 @@ modflow_sources = files( 'Model' / 'SurfaceWaterFlow' / 'swf-obs.f90', 'Model' / 'SurfaceWaterFlow' / 'swf-oc.f90', 'Model' / 'SurfaceWaterFlow' / 'swf-flw.f90', + 'Model' / 'SurfaceWaterFlow' / 'swf-evp.f90', 'Model' / 'SurfaceWaterFlow' / 'swf-pcp.f90', 'Model' / 'SurfaceWaterFlow' / 'swf-sto.f90', 'Model' / 'SurfaceWaterFlow' / 'swf-cdb.f90', diff --git a/utils/idmloader/dfns.txt b/utils/idmloader/dfns.txt index de3f3604364..6a35fec43c2 100644 --- a/utils/idmloader/dfns.txt +++ b/utils/idmloader/dfns.txt @@ -54,6 +54,7 @@ swf-cdb.dfn swf-chd.dfn swf-flw.dfn swf-pcp.dfn +swf-evp.dfn swf-sto.dfn swf-zdg.dfn @@ -67,6 +68,7 @@ chf-cdb.dfn chf-chd.dfn chf-flw.dfn chf-pcp.dfn +chf-evp.dfn chf-sto.dfn chf-zdg.dfn @@ -80,6 +82,7 @@ olf-cdb.dfn olf-chd.dfn olf-flw.dfn olf-pcp.dfn +olf-evp.dfn olf-sto.dfn olf-zdg.dfn From 9bb19f348512e896011e19b7872bea5a9be669fa Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D." Date: Tue, 28 Jan 2025 14:59:31 -0600 Subject: [PATCH 2/4] better evap add test --- autotest/test_chf_evp.py | 222 +++++++++++++++++++++++++ src/Model/SurfaceWaterFlow/swf-evp.f90 | 79 +++++++-- 2 files changed, 283 insertions(+), 18 deletions(-) create mode 100644 autotest/test_chf_evp.py diff --git a/autotest/test_chf_evp.py b/autotest/test_chf_evp.py new file mode 100644 index 00000000000..217b4efe6f1 --- /dev/null +++ b/autotest/test_chf_evp.py @@ -0,0 +1,222 @@ +""" + +Test the evaporation package for the channel model. +This example is tough to converge, because evaporation +removes all of the water in the cell. The model did +not converge well until newton unrelaxation was added +and delta-bar-delta underrelaxation was also activated. + +""" + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +cases = [ + "chf-evp01", +] + +# grid size +dx = 1.0 +nreach = 2 +total_length = dx * nreach +vertices = [] +vertices = [[j, j * dx, 0.0] for j in range(nreach + 1)] +cell1d = [] +for j in range(nreach): + cell1d.append([j, 0.5, 2, j, j + 1]) +nodes = len(cell1d) +nvert = len(vertices) +xorigin = 0.0 +yorigin = 0.0 +angrot = 0.0 + +idomain = np.ones((nodes,), dtype=int) + + +def build_models(idx, test): + perlen = [1] # 1 second + nstp = [10] + tsmult = [1.0] + nper = len(perlen) + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = "chf" + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=f"{name}_sim", version="mf6", exe_name="mf6", sim_ws=ws + ) + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, time_units="SECONDS", nper=nper, perioddata=tdis_rc + ) + + nouter, ninner = 100, 50 + hclose, rclose, relax = 1e-6, 1e-8, 1.0 + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + under_relaxation_theta=0.9, + under_relaxation_kappa=0.0001, + under_relaxation_gamma=0.0, + inner_maximum=ninner, + inner_dvclose=hclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + # backtracking_number=5, + # backtracking_tolerance=1.0, + # backtracking_reduction_factor=0.3, + # backtracking_residual_limit=100.0, + ) + + add_chf_model_disv1d(sim) + + return sim, None + + +def add_chf_model_disv1d(sim): + name = "channel" + chf = flopy.mf6.ModflowChf( + sim, modelname=name, save_flows=True, newtonoptions="UNDER_RELAXATION" + ) + + disv1d = flopy.mf6.ModflowChfdisv1D( + chf, + nodes=nodes, + nvert=nvert, + width=1.0, + bottom=0.0, + idomain=idomain, + vertices=vertices, + cell1d=cell1d, + xorigin=xorigin, + yorigin=yorigin, + angrot=angrot, + ) + + dfw = flopy.mf6.ModflowChfdfw( + chf, + print_flows=True, + save_flows=True, + length_conversion=1.0, + time_conversion=86400.0, + manningsn=0.035, + idcxs=0, + ) + + xfraction = np.array([0.0, 1.0, 2.0]) + height = [1.0, 0.0, 1.0] + npts = len(height) + mannfraction = npts * [1.0] + cxsdata = list(zip(xfraction, height, mannfraction)) + cxs = flopy.mf6.ModflowChfcxs( + chf, + nsections=1, + npoints=npts, + packagedata=[(0, npts)], + crosssectiondata=cxsdata, + ) + + ic = flopy.mf6.ModflowChfic(chf, strt=1.0) + + sto = flopy.mf6.ModflowChfsto( + chf, + steady_state={0: False}, + transient={0: True}, + ) + + evp_rate = 1.0 + evp_spd = [ + (0, evp_rate), + (1, evp_rate), + ] + + pcp = flopy.mf6.ModflowChfevp( + chf, maxbound=1, print_input=True, print_flows=True, stress_period_data=evp_spd + ) + + # output control + oc = flopy.mf6.ModflowChfoc( + chf, + budget_filerecord=f"{name}.bud", + stage_filerecord=f"{name}.stage", + saverecord=[ + ("STAGE", "ALL"), + ("BUDGET", "ALL"), + ], + printrecord=[ + ("STAGE", "ALL"), + ("BUDGET", "ALL"), + ], + ) + + return + + +def plot_output(idx, test): + import matplotlib.pyplot as plt + + mfsim = test.sims[0] + chf = mfsim.chf[0] + + sobj = chf.output.stage() + times = np.array(sobj.times) + ntimes = times.shape[0] + nreach = 2 + stage = sobj.get_alldata().reshape((ntimes, nreach)) + + bobj = chf.output.budget() + evap_list = bobj.get_data(text="EVP") + evap = [] + for rec in evap_list: + evap.append(rec[0]["q"]) + evap = np.array(evap) + + fig, ax = plt.subplots(1, 1) + ax.plot(times, stage[:, 0], label="stage") + ax.plot(times, -evap, label="Qevap") + plt.legend() + + fname = test.workspace / "results.png" + plt.savefig(fname) + + +def check_output(idx, test): + print(f"evaluating model for case {idx}...") + + mfsim = test.sims[0] + chf = mfsim.chf[0] + + # read binary stage file + sobj = chf.output.stage() + times = np.array(sobj.times) + ntimes = times.shape[0] + nreach = 2 + stage = sobj.get_alldata().reshape((ntimes, nreach)) + answer = np.linspace(0.9, 0, ntimes) + assert np.allclose(stage[:, 0], answer, atol=1.0e-5), "stage is not correct" + + +@pytest.mark.developmode +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets, plot): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + plot=lambda t: plot_output(idx, t) if plot else None, + targets=targets, + ) + test.run() diff --git a/src/Model/SurfaceWaterFlow/swf-evp.f90 b/src/Model/SurfaceWaterFlow/swf-evp.f90 index d16b29f909e..6fca101e6f5 100644 --- a/src/Model/SurfaceWaterFlow/swf-evp.f90 +++ b/src/Model/SurfaceWaterFlow/swf-evp.f90 @@ -7,7 +7,7 @@ module SwfEvpModule use KindModule, only: DP, I4B, LGP use ConstantsModule, only: DZERO, LENFTYPE, LENPACKAGENAME, MAXCHARLEN, & - LINELENGTH, DONE + LINELENGTH, DONE, DHALF, DEM6, DPREC use MemoryHelperModule, only: create_mem_path use BndModule, only: BndType use BndExtModule, only: BndExtType @@ -18,6 +18,7 @@ module SwfEvpModule use BlockParserModule, only: BlockParserType use CharacterStringModule, only: CharacterStringType use SmoothingModule, only: sQSaturation + use MathUtilModule, only: get_perturbation use MatrixBaseModule use GeomUtilModule, only: get_node use BaseDisModule, only: DisBaseType @@ -309,7 +310,11 @@ subroutine evp_cf(this) ! local integer(I4B) :: i integer(I4B) :: node - real(DP) :: qevp + integer(I4B) :: inwt + real(DP) :: q + real(DP) :: qeps + real(DP) :: eps + real(DP) :: derv real(DP) :: evap real(DP) :: rlen real(DP), dimension(:), pointer :: reach_length @@ -334,6 +339,13 @@ subroutine evp_cf(this) cycle end if + ! cycle if dry or constant head + if (this%ibound(node) <= 0) then + this%hcof(i) = DZERO + this%rhs(i) = DZERO + cycle + end if + ! Initialize hcof this%hcof(i) = DZERO @@ -348,29 +360,46 @@ subroutine evp_cf(this) rlen = reach_length(node) end if - ! calculate volumetric evaporation flow in L^3/T - qevp = this%get_qevp(node, rlen, this%xnew(node), evap) + ! Calculate volumetric evaporation flow in L^3/Td and add to rhs + q = -this%get_qevp(node, rlen, this%xnew(node), this%xold(node), evap) + this%rhs(i) = -q - ! rhs contribution (positive value means extraction) - this%rhs(i) = qevp + ! Code for adding newton terms + inwt = 1 + if (inwt == 1) then - ! zero out contribution if cell is inactive or constant head - if (this%ibound(node) <= 0) then - this%rhs(i) = DZERO - cycle + ! calculate perturbed q + eps = get_perturbation(this%xnew(node)) + qeps = -this%get_qevp(node, rlen, this%xnew(node) + eps, & + this%xold(node), evap) + + ! calculate derivative + derv = (qeps - q) / eps + + ! add derivative to hcof and update rhs with derivate contribution + this%hcof(i) = derv + this%rhs(i) = this%rhs(i) + derv * this%xnew(node) end if end do end subroutine evp_cf !> @brief Calculate qevp + !! + !! Calculate qevp for both channel and overland flow grids. + !! Approximate the average water surface width of the channel + !! as wavg = delta A over delta h, and then multiply wavg + !! by reach length to come up with surface water area for the + !! channel. Reduce evaporation when depths are small and shut + !! it off when there is no water in the cell. !< - function get_qevp(this, node, rlen, stage, evaporation) result(qevp) + function get_qevp(this, node, rlen, snew, sold, evaporation) result(qevp) ! dummy class(SwfEvpType) :: this !< this instance integer(I4B), intent(in) :: node !< reduced node number real(DP), intent(in) :: rlen !< length of reach - real(DP), intent(in) :: stage !< stage in reach + real(DP), intent(in) :: snew !< current stage in reach + real(DP), intent(in) :: sold !< previous stage in reach real(DP), intent(in) :: evaporation !< evaporation rate in length per time ! return real(DP) :: qevp @@ -379,14 +408,16 @@ function get_qevp(this, node, rlen, stage, evaporation) result(qevp) real(DP) :: depth real(DP) :: bt real(DP) :: area - real(DP) :: top_width + real(DP) :: anew + real(DP) :: aold + real(DP) :: denom + real(DP) :: wavg real(DP) :: width_channel real(DP) :: dummy real(DP) :: qmult ! calculate depth of water bt = this%dis%bot(node) - depth = stage - bt ! Determine the water surface area if (this%dis%is_2d()) then @@ -396,18 +427,30 @@ function get_qevp(this, node, rlen, stage, evaporation) result(qevp) ! channel case idcxs = this%dfw%idcxs(node) call this%dis%get_flow_width(node, node, 0, width_channel, dummy) - top_width = this%cxs%get_wetted_top_width(idcxs, width_channel, depth) - area = rlen * top_width + + depth = snew - bt + anew = this%cxs%get_area(idcxs, width_channel, depth) + depth = sold - bt + aold = this%cxs%get_area(idcxs, width_channel, depth) + wavg = this%cxs%get_wetted_top_width(idcxs, width_channel, depth) + denom = snew - sold + if (abs(denom) > DPREC) then + wavg = (anew - aold) / (snew - sold) + end if + area = rlen * wavg + end if ! Reduce the evap rate as cell goes dry - qmult = this%get_evap_reduce_mult(stage, bt) + qmult = this%get_evap_reduce_mult(snew, bt) ! calculate volumetric evaporation flow in L^3/T qevp = evaporation * area * qmult end function get_qevp + !> @brief Calculate multiplier to reduce evap as depth goes to zero + !< function get_evap_reduce_mult(this, stage, bottom) result(qmult) ! dummy class(SwfEvpType) :: this @@ -423,7 +466,7 @@ function get_evap_reduce_mult(this, stage, bottom) result(qmult) iflowred = 1 qmult = DONE if (iflowred == 1) then - evap_depth = 0.1D0 + evap_depth = DEM6 tp = bottom + evap_depth qmult = sQSaturation(tp, bottom, stage) end if From 5b2d4b27bca53a7779e35c63542ac16c790893e4 Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D." Date: Wed, 29 Jan 2025 08:01:37 -0600 Subject: [PATCH 3/4] update makefile --- make/makefile | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/make/makefile b/make/makefile index cadd1133491..905a319d484 100644 --- a/make/makefile +++ b/make/makefile @@ -131,6 +131,7 @@ $(OBJDIR)/swf-pcpidm.o \ $(OBJDIR)/swf-namidm.o \ $(OBJDIR)/swf-icidm.o \ $(OBJDIR)/swf-flwidm.o \ +$(OBJDIR)/swf-evpidm.o \ $(OBJDIR)/swf-disv2didm.o \ $(OBJDIR)/swf-disv1didm.o \ $(OBJDIR)/swf-dis2didm.o \ @@ -150,6 +151,7 @@ $(OBJDIR)/olf-pcpidm.o \ $(OBJDIR)/olf-namidm.o \ $(OBJDIR)/olf-icidm.o \ $(OBJDIR)/olf-flwidm.o \ +$(OBJDIR)/olf-evpidm.o \ $(OBJDIR)/olf-disv2didm.o \ $(OBJDIR)/olf-dis2didm.o \ $(OBJDIR)/olf-dfwidm.o \ @@ -199,6 +201,7 @@ $(OBJDIR)/chf-pcpidm.o \ $(OBJDIR)/chf-namidm.o \ $(OBJDIR)/chf-icidm.o \ $(OBJDIR)/chf-flwidm.o \ +$(OBJDIR)/chf-evpidm.o \ $(OBJDIR)/chf-disv1didm.o \ $(OBJDIR)/chf-dfwidm.o \ $(OBJDIR)/chf-cxsidm.o \ @@ -432,6 +435,7 @@ $(OBJDIR)/swf-pcp.o \ $(OBJDIR)/swf-oc.o \ $(OBJDIR)/swf-obs.o \ $(OBJDIR)/swf-flw.o \ +$(OBJDIR)/swf-evp.o \ $(OBJDIR)/swf-cdb.o \ $(OBJDIR)/GridConnection.o \ $(OBJDIR)/DistributedVariable.o \ From d80e441b1c41709785424c270e878cc36ebccc27 Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D." Date: Thu, 30 Jan 2025 14:23:59 -0600 Subject: [PATCH 4/4] add flow reduction variables --- src/Model/SurfaceWaterFlow/swf-evp.f90 | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/Model/SurfaceWaterFlow/swf-evp.f90 b/src/Model/SurfaceWaterFlow/swf-evp.f90 index 6fca101e6f5..5e974249577 100644 --- a/src/Model/SurfaceWaterFlow/swf-evp.f90 +++ b/src/Model/SurfaceWaterFlow/swf-evp.f90 @@ -9,6 +9,7 @@ module SwfEvpModule use ConstantsModule, only: DZERO, LENFTYPE, LENPACKAGENAME, MAXCHARLEN, & LINELENGTH, DONE, DHALF, DEM6, DPREC use MemoryHelperModule, only: create_mem_path + use MemoryManagerModule, only: mem_allocate, mem_deallocate use BndModule, only: BndType use BndExtModule, only: BndExtType use SimModule, only: store_error, store_error_filename, count_errors @@ -37,6 +38,8 @@ module SwfEvpModule type, extends(BndExtType) :: SwfEvpType real(DP), dimension(:), pointer, contiguous :: evaporation => null() !< boundary evaporation array + integer(I4B), pointer :: iflowred => null() !< flag that indicates evaporation will be shut off when depth is less than reduction depth + real(DP), pointer :: reduction_depth => null() !< depth below which evaporation is reduced logical, pointer, private :: read_as_arrays ! pointers to other objects @@ -129,9 +132,13 @@ subroutine evp_allocate_scalars(this) call this%BndExtType%allocate_scalars() ! allocate internal members + call mem_allocate(this%iflowred, 'IFLOWRED', this%memoryPath) + call mem_allocate(this%reduction_depth, 'REDUCTION_DEPTH', this%memoryPath) allocate (this%read_as_arrays) ! Set values + this%iflowred = 1 + this%reduction_depth = DEM6 this%read_as_arrays = .false. end subroutine evp_allocate_scalars @@ -459,15 +466,11 @@ function get_evap_reduce_mult(this, stage, bottom) result(qmult) ! return real(DP) :: qmult ! local - integer(I4B) :: iflowred - real(DP) :: evap_depth real(DP) :: tp - iflowred = 1 qmult = DONE - if (iflowred == 1) then - evap_depth = DEM6 - tp = bottom + evap_depth + if (this%iflowred == 1) then + tp = bottom + this%reduction_depth qmult = sQSaturation(tp, bottom, stage) end if @@ -499,7 +502,6 @@ end subroutine evp_fc !< subroutine evp_da(this) ! modules - use MemoryManagerModule, only: mem_deallocate ! dummy class(SwfEvpType) :: this @@ -507,6 +509,8 @@ subroutine evp_da(this) call this%BndExtType%bnd_da() ! scalars + call mem_deallocate(this%iflowred) + call mem_deallocate(this%reduction_depth) deallocate (this%read_as_arrays) ! arrays