Merge remote-tracking branch 'upstream/develop' (#354)

* refactor(maw): refactor MAW conductance calculation for issue 305 (#310)

Add traps to catch 1) skin factors that are <= 0 when using the SKIN
conductance equation and 2) and negative saturated conductances values.

Closes #305

* feat(tableobj): Add a generic table object for lst file output (#303)

Full implementation for SFR package. Partial implementation for LAK, MAW, and UZF packages.

* feat(tableobj): update MAW package to use tableobj for data output (#315)

Also fix some budget reporting issues in the MAW package

* fix(auxmult): auxmult fix when auxmult and bound are in time series (#316)

* change order of time series interpolation so aux is done first in case it is an auxmult column
* partially addresses #314
* updated release notes

* refactor(BoundaryPackage): Add use of TableObject for print_flows option (#317)

* feat(lnf): update n-point geometry package data

* Revert "feat(lnf): update n-point geometry package data"

This reverts commit 3460cfb.

* fix(lak): revise the way outlet to-mvr flows are stored in budobj (#321)

* fix(lak): revise the way outlet to-mvr flows are stored in budobj

* updated release notes

* refactor(lak/maw/sfr): for CONSTANT features write fixed value instead of DHNOFLO to binary file (#329)

* refactor(lak/maw/sfr): for CONSTANT features write fixed value instead of DHNOFLO to binary file

* update release notes

* fix(maw): new headtable corrected to have mawwells rows instead of maxbound (#332)

* fix(sfr): Remove upstream_fraction check from PACKAGEDATA block (#334)

closes #324

* feat(transport): update develop with some general routines needed for transport (#339)

* feat(transport): update develop with some general routines needed for transport
* fix(mover): correct typo in error message
* closes #336
* refactor(disu): correct type in disu error message
* closes #635

* fix(csub): Allow interbed observations when no interbeds specified (#338)

DNODATA values reported if 1) delay interbed observations are specified
and no delay interbeds are specified and 2) interbed observations are
specified and the number of interbeds is zero. Added test of all
observation types to test_gwf_csub_sk01.py which has no interbeds.

closes #298

* ci(travis): update travis.yml to install miniconda instead of apt python (#346)

* refactor(mf6): standardize output to STDOUT to use istdout variable

* refactor(mf6): standardize output to STDOUT to use istdout variable

Modify mf5to6 to use standard SimVariables.f90 in mf6 src and add
ListType variables in original mf5to6 SimVariables.f90 file to
SimListVariables.f90 and update use statements.

* feat(npf): add save saturation option to NPF (#351)

* changed some things to public in lak, uzf, sfr
* fixed minor bug in uzf
* updated release notes and dfn files
* modified lak to store hlak and clak in the bound array

* refactor(sim_message): Refactor to have independent sim_message (#353)

* Merge remote-tracking branch 'upstream/develop' into refactor-src-mergedev

A few missing components from develop

Co-authored-by: langevin-usgs <langevin@usgs.gov>

autotest/test_gwf_returncodes.py

@@ -147,7 +147,7 @@ def test_mf6_idomain_error():
 
 @raises(RuntimeError)
 def test_unknown_keyword_error():
-    returncode, buff = run_mf6([mf6_exe, 'unknown_keyword'], ws)
+    returncode, buff = run_mf6([mf6_exe, '--unknown_keyword'], ws)
     msg = 'could not run {}'.format('unknown_keyword')
     if returncode != 0:
         err_str = 'mf6: illegal option'

dllsrc/mf6bmi.f90

@@ -10,7 +10,7 @@ module mf6bmi
   use iso_c_binding, only: c_int, c_char, c_double, C_NULL_CHAR, c_loc, c_ptr
   use KindModule, only: DP, I4B
   use ConstantsModule, only: LENORIGIN, LENVARNAME, LENMODELNAME, MAXCHARLEN
-  use SimVariablesModule, only: istdout
+  use SimVariablesModule, only: simstdout, istdout
   use InputOutputModule, only: getunit
   implicit none
 
@@ -30,11 +30,11 @@ function bmi_initialize() result(bmi_status) bind(C, name="initialize")
   !DEC$ ATTRIBUTES DLLEXPORT :: bmi_initialize
     integer(kind=c_int) :: bmi_status
     !
-    ! -- set ISTDUNIT to a physical file unit
+    ! -- set STDOUT to a physical file unit
     istdout = getunit()  
     !
-    ! -- open istdout file mfsim.stdout
-    open(unit=istdout, file='mfsim.stdout')
+    ! -- open stdout file mfsim.stdout
+    open(unit=istdout, file=simstdout)
     !
     ! -- initialize MODFLOW 6
     call Mf6Initialize()

doc/ReleaseNotes/ReleaseNotes.tex

@@ -120,6 +120,7 @@ \section{History}
 \underline{NEW FUNCTIONALITY}
 \begin{itemize}
 \item Added silent command line switch (-s or --silent) that sends all output to the screen (\texttt{STDOUT})  to a file (``mfsim.stdout'').
+\item Added a new option, called SAVE\_SATURATION to the Node Property Flow Package.  When invoked, cell saturation is written to the binary budget file.  The cell saturation can be used by post-processors to determine how much of the cell is saturated without having to know the value for ICELLTYPE or the value for head. If a cell is marked as confined (ICELLTYPE=0) then saturation is always one. If ICELLTYPE is one, then saturation ranges between zero and one.
 \end{itemize}
 
 \underline{BASIC FUNCTIONALITY}

doc/mf6io/gwf/binaryoutput.tex

@@ -397,7 +397,8 @@ \subsubsection{Budget File Contents}
 \texttt{UZF-GWD} & 6 & groundwater discharge to land surface from UZF Package \\
 \texttt{UZF-GWD-TO-MVR} & 6 & groundwater discharge to land surface from UZF Package that is routed to Mover Package\\
 \texttt{FLOW-JA-FACE} & 6 & flow to or from a cell in another GWF Model; TXT1ID1 is the name of the GWF Model described by this budget file, TXT2ID1 is the name of the GWF-GWF Exchange, TXT1ID2 is the name of the connected GWF Model, TXT2ID2 is the name of the GWF-GWF Exchange, and ID2 is the cell or node number of the cell in the connected model \\
-\texttt{DATA-SPDIS} & 6 & specific discharge at the cell center.  The x, y, and z components are stored in auxiliary variables called ``qx'', ``qy'', and ``qz'', respectively.   The flow value written for each cell is zero.  The ``DATA'' prefix on the text identifier can be used by post-processors to recognize that the record does not contain a cell flow budget term.
+\texttt{DATA-SPDIS} & 6 & specific discharge at the cell center.  The x, y, and z components are stored in auxiliary variables called ``qx'', ``qy'', and ``qz'', respectively.   The flow value written for each cell is zero.  The ``DATA'' prefix on the text identifier can be used by post-processors to recognize that the record does not contain a cell flow budget term. \\
+\texttt{DATA-SAT} & 6 & cell saturation.  The cell saturation is stored in an auxiliary variable called ``sat''.   The flow value written for each cell is zero.  The ``DATA'' prefix on the text identifier can be used by post-processors to recognize that the record does not contain a cell flow budget term.  The cell saturation can be used by post-processors to determine how much of the cell is saturated without having to know the value for ICELLTYPE or the value for head.  If a cell is marked as confined (ICELLTYPE=0) then saturation is always one.  If ICELLTYPE is one, then saturation ranges between zero and one.  For Newton GWF simulations, saturation is zero if the head is below the cell bottom.
 \label{table:gwfbud}
 \end{longtable}
 

doc/mf6io/mf6ivar/dfn/gwf-npf.dfn

@@ -6,7 +6,7 @@ type keyword
 reader urword
 optional true
 longname keyword to save NPF flows
-description keyword to indicate that cell-by-cell flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control.
+description keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control.
 
 block options
 name alternative_cell_averaging
@@ -133,7 +133,15 @@ type keyword
 reader urword
 optional true
 longname keyword to save specific discharge
-description keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the cell-by-cell flow file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable.
+description keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable.
+
+block options
+name save_saturation
+type keyword
+reader urword
+optional true
+longname keyword to save saturation
+description keyword to indicate that cell saturation will be written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  Saturation will be saved to the budget file as an auxiliary variable saved with the DATA-SAT text label.  Saturation is a cell variable that ranges from zero to one and can be used by post processing programs to determine how much of a cell volume is saturated.  If ICELLTYPE is 0, then saturation is always one.
 
 block options
 name k22overk

doc/mf6io/mf6ivar/md/mf6ivar.md

@@ -146,7 +146,7 @@
 | GWF | DISU | CELL2D | NCVERT | INTEGER | is the number of vertices required to define the cell.  There may be a different number of vertices for each cell. |
 | GWF | DISU | CELL2D | ICVERT | INTEGER (NCVERT) | is an array of integer values containing vertex numbers (in the VERTICES block) used to define the cell.  Vertices must be listed in clockwise order. |
 | GWF | IC | GRIDDATA | STRT | DOUBLE PRECISION (NODES) | is the initial (starting) head---that is, head at the beginning of the GWF Model simulation.  STRT must be specified for all simulations, including steady-state simulations. One value is read for every model cell. For simulations in which the first stress period is steady state, the values used for STRT generally do not affect the simulation (exceptions may occur if cells go dry and (or) rewet). The execution time, however, will be less if STRT includes hydraulic heads that are close to the steady-state solution.  A head value lower than the cell bottom can be provided if a cell should start as dry. |
-| GWF | NPF | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that cell-by-cell flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control. |
+| GWF | NPF | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control. |
 | GWF | NPF | OPTIONS | ALTERNATIVE_CELL_AVERAGING | STRING | is a text keyword to indicate that an alternative method will be used for calculating the conductance for horizontal cell connections.  The text value for ALTERNATIVE\_CELL\_AVERAGING can be ``LOGARITHMIC'', ``AMT-LMK'', or ``AMT-HMK''.  ``AMT-LMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and logarithmic-mean hydraulic conductivity.  ``AMT-HMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and harmonic-mean hydraulic conductivity. If the user does not specify a value for ALTERNATIVE\_CELL\_AVERAGING, then the harmonic-mean method will be used.  This option cannot be used if the XT3D option is invoked. |
 | GWF | NPF | OPTIONS | THICKSTRT | KEYWORD | indicates that cells having a negative ICELLTYPE are confined, and their cell thickness for conductance calculations will be computed as STRT-BOT rather than TOP-BOT. |
 | GWF | NPF | OPTIONS | VARIABLECV | KEYWORD | keyword to indicate that the vertical conductance will be calculated using the saturated thickness and properties of the overlying cell and the thickness and properties of the underlying cell.  If the DEWATERED keyword is also specified, then the vertical conductance is calculated using only the saturated thickness and properties of the overlying cell if the head in the underlying cell is below its top.  If these keywords are not specified, then the default condition is to calculate the vertical conductance at the start of the simulation using the initial head and the cell properties.  The vertical conductance remains constant for the entire simulation. |
@@ -158,7 +158,8 @@
 | GWF | NPF | OPTIONS | IHDWET | INTEGER | is a keyword and integer flag that determines which equation is used to define the initial head at cells that become wet.  If IHDWET is 0, h = BOT + WETFCT (hm - BOT). If IHDWET is not 0, h = BOT + WETFCT (THRESH). |
 | GWF | NPF | OPTIONS | XT3D | KEYWORD | keyword indicating that the XT3D formulation will be used.  If the RHS keyword is also included, then the XT3D additional terms will be added to the right-hand side.  If the RHS keyword is excluded, then the XT3D terms will be put into the coefficient matrix.  Use of XT3D will substantially increase the computational effort, but will result in improved accuracy for anisotropic conductivity fields and for unstructured grids in which the CVFD requirement is violated.  XT3D requires additional information about the shapes of grid cells.  If XT3D is active and the DISU Package is used, then the user will need to provide in the DISU Package the angldegx array in the CONNECTIONDATA block and the VERTICES and CELL2D blocks. |
 | GWF | NPF | OPTIONS | RHS | KEYWORD | If the RHS keyword is also included, then the XT3D additional terms will be added to the right-hand side.  If the RHS keyword is excluded, then the XT3D terms will be put into the coefficient matrix. |
-| GWF | NPF | OPTIONS | SAVE_SPECIFIC_DISCHARGE | KEYWORD | keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the cell-by-cell flow file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable. |
+| GWF | NPF | OPTIONS | SAVE_SPECIFIC_DISCHARGE | KEYWORD | keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable. |
+| GWF | NPF | OPTIONS | SAVE_SATURATION | KEYWORD | keyword to indicate that cell saturation will be written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  Saturation will be saved to the budget file as an auxiliary variable saved with the DATA-SAT text label.  Saturation is a cell variable that ranges from zero to one and can be used by post processing programs to determine how much of a cell volume is saturated.  If ICELLTYPE is 0, then saturation is always one. |
 | GWF | NPF | OPTIONS | K22OVERK | KEYWORD | keyword to indicate that specified K22 is a ratio of K22 divided by K.  If this option is specified, then the K22 array entered in the NPF Package will be multiplied by K after being read. |
 | GWF | NPF | OPTIONS | K33OVERK | KEYWORD | keyword to indicate that specified K33 is a ratio of K33 divided by K.  If this option is specified, then the K33 array entered in the NPF Package will be multiplied by K after being read. |
 | GWF | NPF | GRIDDATA | ICELLTYPE | INTEGER (NODES) | flag for each cell that specifies how saturated thickness is treated.  0 means saturated thickness is held constant;  $>$0 means saturated thickness varies with computed head when head is below the cell top; $<$0 means saturated thickness varies with computed head unless the THICKSTRT option is in effect.  When THICKSTRT is in effect, a negative value of icelltype indicates that saturated thickness will be computed as STRT-BOT and held constant. |

doc/mf6io/mf6ivar/tex/gwf-npf-desc.tex

@@ -3,7 +3,7 @@
 \item \textbf{Block: OPTIONS}
 
 \begin{description}
-\item \texttt{SAVE\_FLOWS}---keyword to indicate that cell-by-cell flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control.
+\item \texttt{SAVE\_FLOWS}---keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control.
 
 \item \texttt{alternative\_cell\_averaging}---is a text keyword to indicate that an alternative method will be used for calculating the conductance for horizontal cell connections.  The text value for ALTERNATIVE\_CELL\_AVERAGING can be ``LOGARITHMIC'', ``AMT-LMK'', or ``AMT-HMK''.  ``AMT-LMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and logarithmic-mean hydraulic conductivity.  ``AMT-HMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and harmonic-mean hydraulic conductivity. If the user does not specify a value for ALTERNATIVE\_CELL\_AVERAGING, then the harmonic-mean method will be used.  This option cannot be used if the XT3D option is invoked.
 
@@ -27,7 +27,9 @@
 
 \item \texttt{RHS}---If the RHS keyword is also included, then the XT3D additional terms will be added to the right-hand side.  If the RHS keyword is excluded, then the XT3D terms will be put into the coefficient matrix.
 
-\item \texttt{SAVE\_SPECIFIC\_DISCHARGE}---keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the cell-by-cell flow file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable.
+\item \texttt{SAVE\_SPECIFIC\_DISCHARGE}---keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable.
+
+\item \texttt{SAVE\_SATURATION}---keyword to indicate that cell saturation will be written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  Saturation will be saved to the budget file as an auxiliary variable saved with the DATA-SAT text label.  Saturation is a cell variable that ranges from zero to one and can be used by post processing programs to determine how much of a cell volume is saturated.  If ICELLTYPE is 0, then saturation is always one.
 
 \item \texttt{K22OVERK}---keyword to indicate that specified K22 is a ratio of K22 divided by K.  If this option is specified, then the K22 array entered in the NPF Package will be multiplied by K after being read.
 

doc/mf6io/mf6ivar/tex/gwf-npf-options.dat

@@ -7,6 +7,7 @@ BEGIN OPTIONS
   [REWET WETFCT <wetfct> IWETIT <iwetit> IHDWET <ihdwet>]
   [XT3D [RHS]]
   [SAVE_SPECIFIC_DISCHARGE]
+  [SAVE_SATURATION]
   [K22OVERK]
   [K33OVERK]
 END OPTIONS

src/Model/GroundWaterFlow/gwf3.f90

@@ -107,7 +107,7 @@ subroutine gwf_cr(filename, id, modelname, smr)
     use ListsModule,                only: basemodellist
     use BaseModelModule,            only: AddBaseModelToList
     use SimModule,                  only: ustop, store_error, count_errors
-    use InputOutputModule,          only: write_centered
+    use GenericUtilitiesModule,     only: write_centered
     use ConstantsModule,            only: LINELENGTH, LENPACKAGENAME
     use VersionModule,              only: VERSION, MFVNAM, MFTITLE,             &
                                           FMTDISCLAIMER, IDEVELOPMODE
@@ -164,19 +164,20 @@ subroutine gwf_cr(filename, id, modelname, smr)
     call namefile_obj%openlistfile(this%iout)
     !
     ! -- Write title to list file
-    call write_centered('MODFLOW'//MFVNAM, this%iout, 80)
-    call write_centered(MFTITLE, this%iout, 80)
-    call write_centered('GROUNDWATER FLOW MODEL (GWF)', this%iout, 80)
-    call write_centered('VERSION '//VERSION, this%iout, 80)
+    call write_centered('MODFLOW'//MFVNAM, 80, iunit=this%iout)
+    call write_centered(MFTITLE, 80, iunit=this%iout)
+    call write_centered('GROUNDWATER FLOW MODEL (GWF)', 80, iunit=this%iout)
+    call write_centered('VERSION '//VERSION, 80, iunit=this%iout)
     !
     ! -- Write if develop mode
-    if (IDEVELOPMODE == 1) call write_centered('***DEVELOP MODE***',           &
-      this%iout, 80)
+    if (IDEVELOPMODE == 1) then
+      call write_centered('***DEVELOP MODE***', 80, iunit=this%iout)
+    end if
     !
     ! -- Write compiler version
     call get_compiler(compiler)
-    call write_centered(' ', this%iout, 80)
-    call write_centered(trim(adjustl(compiler)), this%iout, 80)
+    call write_centered(' ', 80, iunit=this%iout)
+    call write_centered(trim(adjustl(compiler)), 80, iunit=this%iout)
     !
     ! -- Write disclaimer
     write(this%iout, FMTDISCLAIMER)