\newcommand{\pfkey}[4]{\vspace{0.15in} {\noindent{\normalsize {\sl #1} 
\hspace{0.2in}{\bfseries #2}\hspace{0.2in}[#3]}} \newline \indent #4 \newline 
Example Useage: \newline \vspace{-0.25in}}

%=============================================================================
%=============================================================================

\chapter{ParFlow Files}
\label{ParFlow Files}

In this chapter, we discuss the various file formats used in \parflow{}.
To help simplify the description of these formats, we use a pseudocode
notation composed of {\em fields} and {\em control constructs}.

A field is a piece of data having one of the {\em field types} listed
in Table~\ref{table-field-types}
(note that field types may have one meaning in ASCII files and
another meaning in binary files).
%
\begin{table} \center
\caption{Field types.}
\smallskip
\begin{tabular}{||c||c|c||}
\hline
field type & ASCII & binary \\
\hline\hline
\code{integer} & integer & XDR integer        \\ \hline
\code{real}    & real    & -                  \\ \hline
\code{string}  & string  & -                  \\ \hline
\code{double}  & -       & IEEE 8 byte double \\ \hline
\code{float}   & -       & IEEE 4 byte float   \\ \hline
\end{tabular}
\label{table-field-types}
\end{table}
%
Fields are denoted by enclosing the field name with a \code{<} on the
left and a \code{>} on the right.
The field name is composed of alphanumeric characters and
underscores (\code{_}).
In the defining entry of a field, the field name is also prepended by
its field type and a \code{:}.
The control constructs used in our pseudocode have the keyword names
\code{FOR}, \code{IF}, and \code{LINE},
and the beginning and end of each of these constructs is delimited by the
keywords \code{BEGIN} and \code{END}.

The \code{FOR} construct is used to describe repeated input format patterns.
For example, consider the following file format:
\begin{display}\begin{verbatim}
<integer : num_coordinates>
FOR coordinate = 0 TO <num_coordinates> - 1
BEGIN
   <real : x>  <real : y>  <real : z>
END
\end{verbatim}\end{display}
The field \code{<num_coordinates>} is an integer specifying the number
of coordinates to follow.
The \code{FOR} construct indicates that \code{<num_coordinates>} entries
follow, and each entry is composed of the three real fields,
\code{<x>}, \code{<y>}, and \code{<z>}.
Here is an example of a file with this format:
\begin{display}\begin{verbatim}
3
2.0 1.0 -3.5
1.0 1.1 -3.1
2.5 3.0 -3.7
\end{verbatim}\end{display}

The \code{IF} construct is actually an \code{IF}/\code{ELSE} construct,
and is used to describe input format patterns that appear only under
certain circumstances.
For example, consider the following file format:
\begin{display}\begin{verbatim}
<integer : type>
IF (<type> = 0)
BEGIN
   <real : x>  <real : y>  <real : z>
END
ELSE IF (<type> = 1)
BEGIN
   <integer : i>  <integer : j>  <integer : k>
END
\end{verbatim}\end{display}
The field \code{<type>} is an integer specifying the ``type'' of
input to follow.
The \code{IF} construct indicates that if \code{<type>} has value 0,
then the three real fields, \code{<x>}, \code{<y>}, and \code{<z>}, follow.
If \code{<type>} has value 1, then the three integer fields,
\code{<i>}, \code{<j>}, and \code{<k>}, follow.
Here is an example of a file with this format:
\begin{display}\begin{verbatim}
0
2.0 1.0 -3.5
\end{verbatim}\end{display}

The \code{LINE} construct indicates fields that are on the same line
of a file.
Since input files in \parflow{} are all in ``free format'', it is
used only to describe some output file formats.
For example, consider the following file format:
\begin{display}\begin{verbatim}
LINE
BEGIN
   <real : x>
   <real : y>
   <real : z>
END
\end{verbatim}\end{display}
The \code{LINE} construct indicates that the three real fields,
\code{<x>}, \code{<y>}, and \code{<z>}, are all on the same line.
Here is an example of a file with this format:
\begin{display}\begin{verbatim}
2.0 1.0 -3.5
\end{verbatim}\end{display}

Comment lines may also appear in our file format pseudocode.
All text following a \code{#} character is a comment,
and is not part of the file format.

%=============================================================================
%=============================================================================

\section{Main Input File (.tcl)}
\label{Main Input File (.tcl)}

The main \parflow{} input file is a TCL script.  This might seem
overly combersome at first but the basic input file structure is not
very complicated (although it is somewhat verbose).  For more advanced
users, the TCL scripting means you can very easily create programs to
run \parflow{}.  A simple example is creating a loop to run several
hundred different simulations using different seeds to the random field
generators.  This can be automated from within the \parflow{} input file.

The basic idea behind \parflow{} input is a simple database.  The
database contains entries which have a key and a value associated with
that key.  This is very similiar in nature to the Windows XP/Vista
registry and several other systems.  When \parflow{} runs, it
queries the database you have created by key names to get the values
you have specified.

The command \code{pfset} is used to create the database entries.  A
simple \parflow{} input script contains a long list of
\code{pfset} commands.

It should be noted that the keys are ``dynamic'' in that many are
built up from values of other keys.  For example if you have two wells
named {\em northwell} and {\em southwell} then you will have to set some keys
which specify the parameters for each well.  The keys are built up in
a simple sort of heirarchy.

The following sections contain a description of all of the keys used by
\parflow{}.  For an example of input files you can look at the
\file{test} subdirectory of the \parflow{} distribution.  Looking
over some examples should give you a good feel for how the
file scripts are put together.

Each key's entry has the form:

\pfkey{type}{KeyName}{default value}
{Description}{}

\vspace{0.25in}
The ``type'' is one of integer, double, string, list.  Integer and double are
IEEE numbers.  String is a text string (for example, a filename).
Strings can contain spaces if you use the proper TCL syntax
(i.e. using double quotes).  These types are standard TCL
types.  Lists are strings but they indicate the names of a series of items.
For example you might need to specify the names of the geometries. You
would do this using space seperated names (what we are calling a list)
``layer1 layer2 layer3''.

The descriptions that follow are organized into functional areas.  An
example for each database entry is given.

Note that units used for each physical quantity specified in the input file
must be consistent with units used for all other quantities.  The exact units
used can be any consistent set as \parflow{} does not assume any specific set
of units.  However, it is up to the user to make sure all specifications are
indeed consistent.

%=============================================================================
%=============================================================================

\subsection{Input File Format Number}
\label{Input File Format Number}

\pfkey{integer}{FileVersion}{no default}
{This gives the value of
the input file version number that this file fits.}
\begin{display}\begin{verbatim} 
pfset FileVersion 4
\end{verbatim}\end{display} 

As development of the \parflow{} code continues, the input file format
will vary.  We have thus included an input file format number as a way
of verifying that the correct format type is being used.  The user can
check in the \code{parflow/config/file_versions.h} file to verify that
the format number specified in the input file matches the defined
value of \kbd{PFIN_VERSION}. 

%=============================================================================
%=============================================================================

\subsection{Computing Topology}
\label{Computing Topology}

This section describes how processors are assigned in order to solve the domain
in parallel. “P” allocates the number of processes to the grid-cells in x. “Q” allocates
the number of processes to the grid-cells in y. “R” allocates the number of processes to
the grid-cells in z. Please note “R” should always be 1 if you are running with Solver
Richards \cite{Jones-Woodward01} unless you’re running a totally saturated domain (solver 
IMPES).

\pfkey{integer}{Process.Topology.P}{no default}
{This assigns the process splits in the \emph{x} direction.}
\begin{display}\begin{verbatim} 
pfset Process.Topology.P        2
\end{verbatim}\end{display} 

\pfkey{integer}{Process.Topology.Q}{no default}
{This assigns the process splits in the \emph{y} direction.}
\begin{display}\begin{verbatim} 
pfset Process.Topology.Q        1
\end{verbatim}\end{display} 

\pfkey{integer}{Process.Topology.P}{no default}
{This assigns the process splits in the \emph{z} direction.}
\begin{display}\begin{verbatim} 
pfset Process.Topology.R        1
\end{verbatim}\end{display} 

In addition, you can assign the computing topology when you initiate your parflow script using tcl. 
You must include the topology allocation when using tclsh and the parflow script.

Example Usage: \begin{verbatim}
[from Terminal] tclsh default_single.tcl 2 1 1

[At the top of default_single.tcl you must include the following]
set NP  [lindex $argv 0]
set NQ  [lindex $argv 1]

pfset Process.Topology.P        $NP
pfset Process.Topology.Q        $NQ
pfset Process.Topology.R        1 \end{verbatim}

%=============================================================================
%=============================================================================

\subsection{Computational Grid}
\label{Computational Grid}

The computational grid is briefly described in \S~\ref{Defining the Problem}. The computational grid keys set
the bottom left corner of the domain to a specific point in space. If using 
a .pfsol file, the bottom left corner location of the .pfsol file must be the points
designated in the computational grid. The user can also assign the \emph{x}, \emph{y} 
and \emph{z} location to correspond to a specific coordinate system (i.e. UTM).

\pfkey{double}{ComputationalGrid.Lower.X}{no default}
{This assigns the lower \emph{x} coordinate location for the computational grid.}
\begin{display}\begin{verbatim} 
pfset   ComputationalGrid.Lower.X  0.0
\end{verbatim}\end{display} 

\pfkey{double}{ComputationalGrid.Lower.Y}{no default}
{This assigns the lower \emph{y} coordinate location for the computational grid.}
\begin{display}\begin{verbatim} 
pfset   ComputationalGrid.Lower.Y  0.0
\end{verbatim}\end{display} 

\pfkey{double}{ComputationalGrid.Lower.Z}{no default}
{This assigns the lower \emph{z} coordinate location for the computational grid.}
\begin{display}\begin{verbatim} 
pfset   ComputationalGrid.Lower.Z  0.0
\end{verbatim}\end{display} 

\pfkey{integer}{ComputationalGrid.NX}{no default}
{This assigns the number of grid cells in the \emph{x} direction for the computational grid.}
\begin{display}\begin{verbatim} 
pfset  ComputationalGrid.NX  10 
\end{verbatim}\end{display} 

\pfkey{integer}{ComputationalGrid.NY}{no default}
{This assigns the number of grid cells in the \emph{y} direction for the computational grid.}
\begin{display}\begin{verbatim} 
pfset  ComputationalGrid.NY  10 
\end{verbatim}\end{display} 

\pfkey{integer}{ComputationalGrid.NZ}{no default}
{This assigns the number of grid cells in the \emph{z} direction for the computational grid.}
\begin{display}\begin{verbatim} 
pfset  ComputationalGrid.NZ  10 
\end{verbatim}\end{display} 

\pfkey{real}{ComputationalGrid.DX}{no default}
{This defines the size of grid cells in the \emph{x} direction. Units are \emph{L} and are 
defined by the units of the hydraulic conductivity used in the problem.}
\begin{display}\begin{verbatim} 
pfset  ComputationalGrid.DX  10.0 
\end{verbatim}\end{display} 

\pfkey{real}{ComputationalGrid.DY}{no default}
{This defines the size of grid cells in the \emph{y} direction. Units are \emph{L} and are 
defined by the units of the hydraulic conductivity used in the problem.}
\begin{display}\begin{verbatim} 
pfset  ComputationalGrid.DY  10.0 
\end{verbatim}\end{display} 

\pfkey{real}{ComputationalGrid.DZ}{no default}
{This defines the size of grid cells in the \emph{z} direction. Units are \emph{L} and are 
defined by the units of the hydraulic conductivity used in the problem.}
\begin{display}\begin{verbatim} 
pfset  ComputationalGrid.DZ  1.0 
\end{verbatim}\end{display} 


Example Usage: \begin{verbatim}
#---------------------------------------------------------
# Computational Grid
#---------------------------------------------------------
pfset ComputationalGrid.Lower.X	-10.0
pfset ComputationalGrid.Lower.Y     10.0
pfset ComputationalGrid.Lower.Z	1.0

pfset ComputationalGrid.NX		18
pfset ComputationalGrid.NY		18
pfset ComputationalGrid.NZ		8

pfset ComputationalGrid.DX		8.0
pfset ComputationalGrid.DY		10.0
pfset ComputationalGrid.DZ		1.0
\end{verbatim}


%=============================================================================
%=
\subsection{Geometries}
\label{Geometries}
\pfaddanchor{pfinGeometries}

Here we define all ``geometrical'' information needed by \parflow{}.
For example, the domain (and patches on the domain where boundary
conditions are to be imposed), lithology or hydrostratigraphic units,
faults, initial plume shapes, and so on, are considered geometries.

This input section is a little confusing.  Two items are being
specified, geometry inputs and geometries.  A geometry input is a type
of geometry input (for example a box or an input file).  A geometry input
can contain more than one geometry.  A geometry input of type Box
has a single geometry (the square box defined by the extants of the
two points).  A SolidFile input type can contain several geometries.

\pfkey{list}{GeomInput.Names}{no default}
{
This is a list of the geometry input names which define the containers
for all of the geometries defined for this problem.
}
\begin{display}\begin{verbatim}
pfset GeomInput.Names    "solidinput indinput boxinput"
\end{verbatim}\end{display}

\pfkey{string}{GeomInput.{\em geom\_input\_name}.InputType}{no default}
{
This defines the input type for the geometry input with
{\em geom\_input\_name}.  
This key must be one of: {\bf SolidFile, IndicatorField}, {\bf Box}.
}
\begin{display}\begin{verbatim}
pfset GeomInput.solidinput.InputType  SolidFile
\end{verbatim}\end{display}

\pfkey{list}{GeomInput.{\em geom\_input\_name}.GeomNames}{no default}
{
This is a list of the names of the geometries defined by the geometry
input.  For a geometry input type of Box, the list should contain a
single geometry name.  For the SolidFile geometry type this should
contain a list with the same number of gemetries as were defined
using GMS.  The order of geometries in the SolidFile should match the
names.  For IndicatorField types you need to specify the value in the
input field which matches the name using
GeomInput.{\em geom\_input\_name}.Value.
}
\begin{display}\begin{verbatim}
pfset GeomInput.solidinput.GeomNames "domain bottomlayer \
                                      middlelayer toplayer"
\end{verbatim}\end{display}

\pfkey{string}{GeomInput.{\em geom\_input\_name}.Filename}{no default}
{
For IndicatorField and SolidFile geometry inputs this key specifies
the input filename which contains the field or solid information.
}
\begin{display}\begin{verbatim}
pfset GeomInput.solidinput.FileName   ocwd.pfsol
\end{verbatim}\end{display}


\pfkey{integer}{GeomInput.{\em geometry\_input\_name}.Value}{no default}
{
For IndicatorField geometry inputs you need to specify the mapping between
values in the input file and the geometry names.  The named geometry will
be defined whereever the input file is equal to the specifed value.
}
\begin{display}\begin{verbatim}
pfset GeomInput.sourceregion.Value   11
\end{verbatim}\end{display}

For box geometries you need to specify the location of the box.  This
is done by defining two corners of the the box.

\pfkey{double}{Geom.{\em box\_geom\_name}.Lower.X}{no  default}
{This gives the lower X real space coordinate value
of the previously specified box geometry of name {\em box\_geom\_name}.}
\begin{display}\begin{verbatim}
pfset Geom.background.Lower.X   -1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em box\_geom\_name}.Lower.Y}{no default}
{This gives the lower Y real space coordinate value
of the previously specified box geometry of name {\em box\_geom\_name}.}
\begin{display}\begin{verbatim}
pfset Geom.background.Lower.Y   -1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em box\_geom\_name}.Lower.Z}{no default}
{This gives the lower Z real space coordinate value
of the previously specified box geometry of name {\em box\_geom\_name}.}
\begin{display}\begin{verbatim}
pfset Geom.background.Lower.Z   -1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em box\_geom\_name}.Upper.X}{no default}
{This gives the upper X real space coordinate value
of the previously specified box geometry of name {\em box\_geom\_name}.}
\begin{display}\begin{verbatim}
pfset Geom.background.Upper.X   151.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em box\_geom\_name}.Upper.Y}{no default}
{This gives the upper Y real space coordinate value
of the previously specified box geometry of name {\em box\_geom\_name}.}
\begin{display}\begin{verbatim}
pfset Geom.background.Upper.Y   171.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em box\_geom\_name}.Upper.Z}{no default}
{This gives the upper Z real space coordinate value
of the previously specified box geometry of name {\em box\_geom\_name}.}
\begin{display}\begin{verbatim}
pfset Geom.background.Upper.Z   11.0
\end{verbatim}\end{display}

\pfkey{list}{Geom.{\em geom\_name}.Patches}{no default}
{
Patches are defined on the surfaces of geometries.  Currently you can
only define patches on Box geometries and on the the first geometry in a
SolidFile.  For a Box the order is fixed (left right front back bottom
top) but you can name the sides anything you want.

For SolidFiles the order is printed by the conversion routine that
converts GMS to SolidFile format.  
}
\begin{display}\begin{verbatim}
pfset Geom.background.Patches   "left right front back bottom top"
\end{verbatim}\end{display}

Here is an example geometry input section which has three geometry inputs.

\begin{display}\begin{verbatim}
#---------------------------------------------------------
# The Names of the GeomInputs
#---------------------------------------------------------
pfset GeomInput.Names 			"solidinput indinput boxinput"
#
# For a solid file geometry input type you need to specify the names
# of the gemetries and the filename
#

pfset GeomInput.solidinput.InputType	SolidFile

# The names of the geometries contained in the solid file. Order is
# important and defines the mapping. First geometry gets the first name. 
pfset GeomInput.solidinput.GeomNames	"domain"
#
# Filename that contains the geometry
#

pfset GeomInput.solidinput.FileName 	ocwd.pfsol

#
# An indicator field is a 3D field of values. 
# The values within the field can be mapped 
# to ParFlow geometries. Indicator fields must match the
# computation grid exactly!
#

pfset GeomInput.indinput.InputType 		IndicatorField
pfset GeomInput.indinput.GeomNames    	“sourceregion concenregion”
pfset GeomInput.indinput.FileName		ocwd.pfb

#
# Within the indicator.pfb file, assign the values to each GeomNames
# 
pfset GeomInput.sourceregion.Value 	11
pfset GeomInput.concenregion.Value 	12

#
# A box is just a box defined by two points.
#

pfset GeomInput.boxinput.InputType	Box
pfset GeomInput.boxinput.GeomName	background
pfset Geom.background.Lower.X 		-1.0
pfset Geom.background.Lower.Y 		-1.0
pfset Geom.background.Lower.Z 		-1.0
pfset Geom.background.Upper.X 		151.0
pfset Geom.background.Upper.Y 		171.0
pfset Geom.background.Upper.Z 		11.0

#
# The patch order is fixed in the .pfsol file, but you 
# can call the patch name anything you 
# want (i.e. left right front back bottom top)
#

pfset Geom.domain.Patches             		" z-upper x-lower y-lower \
                                      			x-upper y-upper z-lower"

\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Timing Information}
\label{Timing Information}
\pfaddanchor{pfinTiming}

The data given in the timing section describe all the ``temporal''
information needed by \parflow{}.  The data items are used to
describe time units for later sections, sequence iterations in time,
indicate actual starting and stopping values and give instructions
on when data is printed out.

\pfkey{double}{TimingInfo.BaseUnit}{no default} { This key is used to
  indicate the base unit of time for entering time values.  All time
  should be expressed as a multiple of this value.  This should be set
  to the smallest interval of time to be used in the problem.  For
  example, a base unit of ``1'' means that all times will be integer
  valued.  A base unit of ``0.5'' would allow integers and fractions
  of 0.5 to be used for time input values.

  The rationale behind this restriction is to allow time to be
  discretized on some interval to enable integer arithmetic to be used
  when computing/comparing times.  This avoids the problems associated
  with real value comparisons which can lead to events occurring at
  different timesteps on different architectures or compilers.

  This value is also used when describing ``time cycling data'' in,
  currently, the well and boundary condition sections.  The lengths of
  the cycles in those sections will be integer multiples of this
  value, therefore it needs to be the smallest divisor which produces
  an integral result for every ``real time'' cycle interval length
  needed.  
}
\begin{display}\begin{verbatim}
pfset TimingInfo.BaseUnit      1.0
\end{verbatim}\end{display}

\pfkey{integer}{TimingInfo.StartCount}{no default}
{
This key is used to indicate the time step number that will be associated
with the first advection cycle in a transient problem.  The value
{\bf -1} indicates that advection is not to be done.  The value {\bf 0}
indicates that advection should begin with the given initial conditions.
Values greater than {\bf 0} are intended to mean ``restart'' from some
previous ``checkpoint'' time-step, but this has not yet been implemented.
}
\begin{display}\begin{verbatim}
pfset TimingInfo.StartCount    0
\end{verbatim}\end{display}

\pfkey{double}{TimingInfo.StartTime}{no default}
{
This key is used to indicate the starting time for the simulation.
}
\begin{display}\begin{verbatim}
pfset TimingInfo.StartTime     0.0
\end{verbatim}\end{display}

\pfkey{double}{TimingInfo.StopTime}{no default}
{
This key is used to indicate the stopping time for the simulation.
}
\begin{display}\begin{verbatim}
pfset TimingInfo.StopTime      100.0
\end{verbatim}\end{display}

\pfkey{double}{TimingInfo.DumpInterval}{no default}
{
This key is the real time interval at which time-dependent output should
be written.  A value of {\bf 0} will produce undefined behavior.  If the
value is negative, output will be dumped out every $n$ time steps, where $n$
is the absolute value of the integer part of the value.
}
\begin{display}\begin{verbatim}
pfset TimingInfo.DumpInterval  10.0
\end{verbatim}\end{display}

\pfkey{integer}{TimingInfo.DumpIntervalExecutionTimeLimit}{0}
{
This key is used to indicate a wall clock time to halt the execution
of a run.  At the end of each dump interval the time remaining in the
batch job is compared with the user supplied value, if remaining time
is less than or equal to the supplied value the execution is halted.
Typically used when running on batch systems with time limits to force
a clean shutdown near the end of the batch job.  Time units is
seconds, a value of {\bf 0} (the default) disables the check.

Currently only supported on SLURM based systems, ``--with-slurm'' must be specified
at configure time to enable.
}
\begin{display}\begin{verbatim}
pfset TimingInfo.DumpIntervalExecutionTimeLimit 360
\end{verbatim}\end{display}

\vspace{0.5in}

For {\em Richards' equation cases only} input is collected for time step
selection.  Input for this section is given as follows:

\pfkey{list}{TimeStep.Type}{no default}
{
This key must be one of: {\bf Constant} or {\bf Growth}.  The value
{\bf Constant} defines a constant time step.  The value {\bf Growth}
defines a time step that starts as $dt_0$ and is defined for
other steps as $dt^{new} = \gamma dt^{old}$ such that $dt^{new} \leq 
dt_{max}$ and $dt^{new} \geq dt_{min}$.
}
\begin{display}\begin{verbatim}
pfset TimeStep.Type      Constant
\end{verbatim}\end{display}

\pfkey{double}{TimeStep.Value}{no default}
{
This key is used only if a constant time step is selected and indicates
the value of the time step for all steps taken.
}
\begin{display}\begin{verbatim}
pfset TimeStep.Value      0.001
\end{verbatim}\end{display}

\pfkey{double}{TimeStep.InitialStep}{no default}
{
This key specifies the initial time step $dt_0$ if the {\bf Growth} type time
step is selected. 
}
\begin{display}\begin{verbatim}
pfset TimeStep.InitialStep    0.001
\end{verbatim}\end{display}

\pfkey{double}{TimeStep.GrowthFactor}{no default}
{
This key specifies the growth factor $\gamma$ by which a time step will be
multiplied to get the new time step when the {\bf Growth} type time step is
selected. 
}
\begin{display}\begin{verbatim}
pfset TimeStep.GrowthFactor      1.5
\end{verbatim}\end{display}

\pfkey{double}{TimeStep.MaxStep}{no default}
{
This key specifies the maximum time step allowed, $dt_{max}$, when the {\bf
Growth} type time step is selected.
}
\begin{display}\begin{verbatim}
pfset TimeStep.MaxStep      86400
\end{verbatim}\end{display}

\pfkey{double}{TimeStep.MinStep}{no default}
{
This key specifies the minimum time step allowed, $dt_{min}$, when the {\bf
Growth} type time step is selected.
}
\begin{display}\begin{verbatim}
pfset TimeStep.MinStep      1.0e-3
\end{verbatim}\end{display}

Here is a detailed example of how timing keys might be used in a simualtion.
\begin{display}\begin{verbatim}
#-----------------------------------------------------------------------------
# Setup timing info [hr]
# 8760 hours in a year. Dumping files every 24 hours. Hourly timestep
#-----------------------------------------------------------------------------
pfset TimingInfo.BaseUnit		1.0
pfset TimingInfo.StartCount		0
pfset TimingInfo.StartTime		0.0
pfset TimingInfo.StopTime		8760.0
pfset TimingInfo.DumpInterval	-24

## Timing constant example
pfset TimeStep.Type			Constant
pfset TimeStep.Value			1.0

## Timing growth example
pfset TimeStep.Type			Growth
pfset TimeStep.InitialStep		0.0001
TimeStep.GrowthFactor		1.4
TimeStep.MaxStep			1.0
TimeStep.MinStep			0.0001
\end{verbatim}\end{display}

%=============================================================================

%=============================================================================
%=============================================================================

\subsection{Time Cycles}
\label{Time Cycles}
\pfaddanchor{pfinTimeCycles}

The data given in the time cycle section describe how time intervals are created and named to be used for time-dependent boundary and well information needed by \parflow{}.  All the time cycles are synched to the {\bf TimingInfo.BaseUnit} key described above and are {\em integer multipliers} of that value.

\pfkey{list}{CycleNames}{no default}
{
This key is used to specify the named time cycles to be used in a simulation.  It is a list of names and each name defines a time cycle and the number of items determines the total number of time cycles specified.  Each named cycle is described using a number of keys defined below.
}
\begin{display}\begin{verbatim}
pfset Cycle.Names constant onoff
\end{verbatim}\end{display}

\pfkey{list}{Cycle.{\em cycle\_name}.Names}{no default}
{
This key is used to specify the named time intervals for each cycle.  It is a list of names and each name defines a time interval when a specific boundary condition is applied and the number of items determines the total number of intervals in that time cycle. 
}
\begin{display}\begin{verbatim}
pfset Cycle.onoff.Names "on off"
\end{verbatim}\end{display}

\pfkey{integer}{Cycle.{\em cycle\_name.interval\_name}.Length}{no default}
{
This key is used to specify the length of a named time intervals.  It is an {\em integer multiplier} of the value set for the {\bf TimingInfo.BaseUnit} key described above. The total length of a given time cycle is the sum of all the intervals multiplied by the base unit.
}
\begin{display}\begin{verbatim}
pfset Cycle.onoff.on.Length             10
\end{verbatim}\end{display}

\pfkey{integer}{Cycle.{\em cycle\_name}.Repeat}{no default}
{
This key is used to specify the how many times a named time interval repeats.  A positive value specifies a number of repeat cycles a value of -1 specifies that the cycle repeat for the entire simulation.
}
\begin{display}\begin{verbatim}
pfset Cycle.onoff.Repeat               -1
\end{verbatim}\end{display}

Here is a detailed example of how time cycles might be used in a simualtion.
\begin{display}\begin{verbatim}

#-----------------------------------------------------------------------------
# Time Cycles
#-----------------------------------------------------------------------------
pfset Cycle.Names 			"constant rainrec"
pfset Cycle.constant.Names		"alltime"
pfset Cycle.constant.alltime.Length	8760
pfset Cycle.constant.Repeat		-1

# Creating a rain and recession period for the rest of year
pfset Cycle.rainrec.Names		"rain rec"
pfset Cycle.rainrec.rain.Length	10
pfset Cycle.rainrec.rec.Length	8750
pfset Cycle.rainrec.Repeat              	-1
\end{verbatim}\end{display}


%=============================================================================
%=============================================================================

\subsection{Domain}
\label{Domain}
\pfaddanchor{pfinDomain}

The domain may be represented by any of the solid types in
\S~\ref{Geometries} above that allow the definition of surface patches.
These surface patches are used to define boundary conditions in
\S~\ref{Boundary Conditions: Pressure} and
\S~\ref{Boundary Conditions: Saturation} below.
Subsequently, it is required that the union (or combination) of the defined surface
patches equal the entire domain surface.  NOTE: This requirement is NOT checked
in the code.

\pfkey{string}{Domain.GeomName}{no default}
{
This key specifies which of the named geometries is the problem domain.
}
\begin{display}\begin{verbatim}
pfset Domain.GeomName    domain
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Phases and Contaminants}
\label{Phases and Contaminants}

\pfkey{list}{Phase.Names}{no default}
{This specifies the names of phases to be modeled. 
Currently only 1 or 2 phases may be modeled.}
\begin{display}\begin{verbatim}
pfset Phase.Names    "water"
\end{verbatim}\end{display}

\pfkey{list}{Contaminant.Names}{no default}
{This specifies the names of contaminants to be advected.}
\begin{display}\begin{verbatim}
pfset Contaminants.Names   "tce"
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Gravity, Phase Density and Phase Viscosity}
\label{Gravity, Phase Density and Phase Viscosity}
\pfaddanchor{pfinGravityPhaseDensityPhaseViscosity}

\pfkey{double}{Gravity}{no default}
{Specifies the gravity constant to be used.}
\begin{display}\begin{verbatim}
pfset Gravity	1.0
\end{verbatim}\end{display}

\pfkey{string}{Phase.{\em phase\_name}.Density.Type}{no default}
{This key specifies whether density will be a constant value or if 
it will be given by an equation of state of the form $(rd)exp(cP)$, 
where $P$ is pressure, $rd$ is
the density at atmospheric pressure, and $c$ is the phase
compressibility constant.
This key must be either {\bf Constant} or {\bf EquationOfState}.
}
\begin{display}\begin{verbatim}
pfset Phase.water.Density.Type	 Constant
\end{verbatim}\end{display}

\pfkey{double}{Phase.{\em phase\_name}.Density.Value}{no default}
{This specifies the value of density if this phase was specified to have a
constant density value for the phase {\em phase\_name}.}
\begin{display}\begin{verbatim}
pfset Phase.water.Density.Value   1.0
\end{verbatim}\end{display}

\pfkey{double}{Phase.{\em phase\_name}.Density.ReferenceDensity}{no default}
{This key specifies the reference density if an equation of state density
function is specified for the phase {\em phase\_name}.}
\begin{display}\begin{verbatim}
pfset Phase.water.Density.ReferenceDensity   1.0
\end{verbatim}\end{display}

\pfkey{double}{Phase.{\em phase\_name}.Density.CompressibilityConstant}
{no default}
{This key specifies the phase compressibility constant 
if an equation of state density function is specified for the phase
{\em phase|-name}.}
\begin{display}\begin{verbatim}
pfset Phase.water.Density.CompressibilityConstant   1.0
\end{verbatim}\end{display}

\pfkey{string}{Phase.{\em phase\_name}.Viscosity.Type}{Constant}
{This key specifies whether viscosity will be a constant value.
Currently, the only choice for this key is {\bf Constant}.}
\begin{display}\begin{verbatim}
pfset Phase.water.Viscosity.Type   Constant
\end{verbatim}\end{display}

\pfkey{double}{Phase.{\em phase\_name}.Viscosity.Value}{no default}
{This specifies the value of viscosity if this phase was specified to have a
constant viscosity value.}
\begin{display}\begin{verbatim}
pfset Phase.water.Viscosity.Value   1.0
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Chemical Reactions}
\label{Chemical Reactions}
\pfaddanchor{pfinChemicalReactions}

\pfkey{double}{Contaminants.{\em contaminant\_name}.Degradation.Value}{no default}
{
This key specifies the half-life decay rate of the named contaminant,
{\em contaminant\_name}.  At present only first order decay reactions
are implemented and it is assumed that one contaminant cannot decay
into another.
}
\begin{display}\begin{verbatim}
pfset Contaminants.tce.Degradation.Value         0.0
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Permeability}
\label{Permeability}
\pfaddanchor{pfinPermeability}

In this section, permeability property values are assigned to grid
points within geometries (specified in \S~\ref{Geometries} above)
using one of the methods described below.  Permeabilities are assumed to be a
diagonal tensor with entries given as,
\[ 
\left( 
\begin{array}{ccc}
k_x({\bf x}) & 0 & 0 \\
0 & k_y({\bf x}) & 0 \\
0 & 0 & k_z({\bf x}) 
\end{array} \right) 
K({\bf x}),
\]
where $K({\bf x})$ is the permeability field given below.  Specification of the
tensor entries ($k_x, k_y$ and $k_z$) will be given at the end of this section.

The random field routines ({\em turning bands} and {\em pgs}) can use
conditioning data if the user so desires.  It is not necessary to use
conditioning as
\parflow{} automatically defaults to not use conditioning data, but
if conditioning is desired, the following key should be set:

\pfkey{string}{Perm.Conditioning.FileName}{``NA''}
{
This key specifies the name of the file that contains the conditioning
data.  The default string {\bf NA} indicates that conditioning data is
not applicable.
}
\begin{display}\begin{verbatim}
pfset Perm.Conditioning.FileName   "well_cond.txt"
\end{verbatim}\end{display}

The file that contains the conditioning data is a simple ascii file
containing points and values. The format is:

\begin{display}\begin{verbatim}
nlines
x1 y1 z1 value1
x2 y2 z2 value2
.  .  .    .
.  .  .    .
.  .  .    .
xn yn zn valuen
\end{verbatim}\end{display}

The value of {\em nlines} is just the number
of lines to follow in the file, which 
is equal to the number of data points.

The variables {\em xi,yi,zi} are the real space coordinates (in the
units used for the given parflow run) of
a point at which a fixed permeability value
is to be assigned. The variable {\em valuei} is the actual permeability
value that is known. 

Note that the coordinates
are not related to the grid in any way. Conditioning
does not require that fixed values be on a 
grid. The PGS algorithm
will map the given value to the closest grid point
and that will be fixed. This is done for speed reasons.
The conditioned turning bands algorithm does not do
this; conditioning is done for every grid point using
the given conditioning data at the location given.
Mapping to grid points for that algorithm does not give
any speedup, so there is no need to do it.

NOTE: The given values should be the actual measured
values - adjustment in the conditioning for the lognormal
distribution that is assumed is taken care of in the algorithms.

The general format for the permeability input is as follows:

\pfkey{list}{Geom.Perm.Names}{no default}
{
This key specifies all of the geometries to which a permeability field 
will be assigned.  These geometries must cover the entire computational
domain.
}
\begin{display}\begin{verbatim}
pfset GeomInput.Names   "background domain concen_region"
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\rm geometry\_name}.Perm.Type}{no default}
{
This key specifies which method is to be used to assign permeability
data to the named geometry, {\em geometry\_name}.  It must be either
{\bf Constant}, {\bf TurnBands}, {\bf ParGuass}, or {\bf PFBFile}.  The
{\bf Constant} value indicates that a constant is to be assigned to all
grid cells within a geometry.  The {\bf TurnBand} value indicates that
Tompson's Turning Bands method is to be used to assign permeability data
to all grid cells within a geometry \cite{TAG89}.  
The {\bf ParGauss} value indicates
that a Parallel Gaussian Simulator method is to be used to assign
permeability data to all grid cells within a geometry.  The
{\bf PFBFile} value indicates that premeabilities are to be read from
the ``ParFlow Binary'' file.  Both the Turning Bands and Parallel
Gaussian Simulators generate a random field with correlation lengths
in the $3$ spatial directions given by $\lambda_x$, $\lambda_y$, and
$\lambda_z$ with the geometric mean of the log normal field given by
$\mu$ and the standard deviation of the normal field given by $\sigma$.
In generating the field both of these methods can be made to stratify
the data, that is follow the top or bottom surface.  The generated field
can also be made so that the data is normal or log normal, with or
without bounds truncation.  Turning Bands uses a line process, the
number of lines used and the resolution of the process can be changed as
well as the maximum normalized frequency $K_{\rm max}$ and the normalized
frequency increment $\delta K$.  The Parallel Gaussian Simulator uses
a search neighborhood, the number of simulated points and the number of
conditioning points can be changed.
}
\begin{display}\begin{verbatim}
pfset Geom.background.Perm.Type   Constant
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.Value}{no default}
{
This key specifies the value assigned to all points in the named
geometry, {\em geometry\_name}, if the type was set to constant.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.Value   1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.LambdaX}{no default}
{
This key specifies the x correlation length, $\lambda_x$, of the field
generated for the named geometry, {\em geometry\_name},
if either the Turning Bands or Parallel Gaussian Simulator are
chosen.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.LambdaX   200.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.LambdaY}{no default}
{
This key specifies the y correlation length, $\lambda_y$, of the field
generated for the named geometry, {\em geometry\_name}, if either the
Turning Bands or Parallel Gaussian Simulator are chosen.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.LambdaY   200.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.LambdaZ}{no default}
{
This key specifies the z correlation length, $\lambda_z$, of the field
generated for the named geometry, {\em geometry\_name}, if either the
Turning Bands or Parallel Gaussian Simulator are chosen.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.LambdaZ   10.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.GeomMean}{no default}
{
This key specifies the geometric mean, $\mu$, of the log normal field
generated for the named geometry, {\em geometry\_name}, if either the
Turning Bands or Parallel Gaussian Simulator are chosen.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.GeomMean   4.56
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.Sigma}{no default}
{
This key specifies the standard deviation, $\sigma$, of the normal field
generated for the named geometry, {\em geometry\_name}, if either the
Turning Bands or Parallel Gaussian Simulator are chosen.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.Sigma   2.08
\end{verbatim}\end{display}

\pfkey{integer}{Geom.{\em geometry\_name}.Perm.Seed}{1}
{
This key specifies the initial seed for the random number generator used
to generate the field for the named geometry, {\em geometry\_name}, if
either the Turning Bands or Parallel Gaussian Simulator are chosen.
This number must be positive.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.Seed   1
\end{verbatim}\end{display}

\pfkey{integer}{Geom.{\em geometry\_name}.Perm.NumLines}{100}
{
This key specifies the number of lines to be used in the Turning Bands
algorithm for the named geometry, {\em geometry\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.NumLines   100
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.RZeta}{5.0}
{
This key specifies the resolution of the line processes, in terms of the
minimum grid spacing, to be used in the Turning Bands algorithm for the
named geometry, {\em geometry\_name}.  Large values imply high resolution.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.RZeta   5.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.KMax}{100.0}
{
This key specifies the the maximum normalized frequency, $K_{\rm max}$,
to be used in the Turning Bands algorithm for the named geometry,
{\em geometry\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.KMax   100.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.DelK}{0.2}
{
This key specifies the normalized frequency increment, $\delta K$, to be
used in the Turning Bands algorithm for the named geometry,
{\em geometry\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.DelK   0.2
\end{verbatim}\end{display}

\pfkey{integer}{Geom.{\em geometry\_name}.Perm.MaxNPts}{no default}
{
This key sets limits on the number of simulated points in the search
neighborhood to be used in the Parallel Gaussian Simulator for the named
geometry, {\em geometry\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.MaxNPts   5
\end{verbatim}\end{display}

\pfkey{integer}{Geom.{\em geometry\_name}.Perm.MaxCpts}{no default}
{
This key sets limits on the number of external conditioning points in the
search neighborhood to be used in the Parallel Gaussian Simulator for
the named geometry, {\em geometry\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.MaxCpts   200
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.Perm.LogNormal}{"LogTruncated"}
{
The key specifies when a normal, log normal, truncated normal or
truncated log normal field is to be generated by the method for the
named geometry, {\em geometry\_name}.  This value must be one of
{\bf Normal}, {\bf Log}, {\bf NormalTruncated} or {\bf LogTruncate}
and can be used with either Turning Bands or the Parallel Gaussian
Simulator.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.LogNormal   "LogTruncated"
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.Perm.StratType}{"Bottom"}
{
This key specifies the stratification of the permeability field
generated by the method for the named geometry, {\em geometry\_name}.
The value must be one of {\bf Horizontal}, {\bf Bottom} or {\bf Top} and
can be used with either the Turning Bands or the Parallel Gaussian
Simulator.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.StratType  "Bottom"
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.LowCutoff}{no default}
{
This key specifies the low cutoff value for truncating the generated
field for the named geometry, {\em geometry\_name}, when either the
NormalTruncated or LogTruncated values are chosen.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.LowCutoff   0.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.HighCutoff}{no default}
{
This key specifies the high cutoff value for truncating the generated
field for the named geometry, {\em geometry\_name}, when either the
NormalTruncated or LogTruncated values are chosen.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.HighCutoff   100.0
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.Perm.FileName}{no default}
{
This key specifies that permeability values for the specified geometry,
{\em geometry\_name}, are given according to a user-supplied description
in the ``ParFlow Binary'' file whose filename is given as the value.
For a description of the ParFlow Binary file format, see
\S~\ref{ParFlow Binary Files (.pfb)}.  The ParFlow Binary file
associated with the named geometry must contain a collection of
permeability values corresponding in a one-to-one manner to the entire
computational grid.  That is to say, when the contents of the file are
read into the simulator, a complete permeability description for the
entire domain is supplied.  Only those values associated with
computational cells residing within the geometry (as it is represented
on the computational grid) will be copied into data structures used
during the course of a simulation.  Thus, the values associated with
cells outside of the geounit are irrelevant.  For clarity, consider a
couple of different scenarios.  For example, the user may create a file
for each geometry such that appropriate permeability values are given
for the geometry and ``garbage" values (e.g., some flag value) are given
for the rest of the computational domain. In this case, a separate
binary file is specified for each geometry.  Alternatively, one may
place all values representing the permeability field on the union of the
geometries into a single binary file.  Note that the permeability values
must be represented in precisely the same configuration as the
computational grid. Then, the same file could be specified for each
geounit in the input file.  Or, the computational domain could be
described as a single geouint (in the ParFlow input file) in which case
the permeability values would be read in only once.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.FileName "domain_perm.pfb"
\end{verbatim}\end{display}

\pfkey{string}{Perm.TensorType}{no default}
{
This key specifies whether the permeability tensor entries $k_x, k_y$ and $k_z$
will be specified as three constants within a set of regions covering the
domain or whether the entries will be specified cell-wise by files.  
The choices
for this key are {\bf TensorByGeom} and {\bf TensorByFile}.  
}
\begin{display}\begin{verbatim}
pfset Perm.TensorType     TensorByGeom
\end{verbatim}\end{display}

\pfkey{string}{Geom.Perm.TensorByGeom.Names}{no default}
{
This key specifies all of the geometries to which permeability tensor entries
will be assigned.  These geometries must cover the entire computational
domain.
}
\begin{display}\begin{verbatim}
pfset Geom.Perm.TensorByGeom.Names   "background domain" 
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.TensorValX}{no default}
{
This key specifies the value of $k_x$ for the geometry given by {\em
geometry\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.TensorValX   1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.TensorValY}{no default}
{
This key specifies the value of $k_y$ for the geometry given by {\em
geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.TensorValY   1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Perm.TensorValZ}{no default}
{
This key specifies the value of $k_z$ for the geometry given by {\em
geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.TensorValZ   1.0
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.Perm.TensorFileX}{no default}
{
This key specifies that $k_x$ values for the specified geometry,
{\em geometry\_name}, are given according to a user-supplied description
in the ``ParFlow Binary'' file whose filename is given as the value.
The only choice for the value of {\em geometry\_name} is ``domain''.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.TensorByFileX   "perm_x.pfb"
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.Perm.TensorFileY}{no default}
{
This key specifies that $k_y$ values for the specified geometry,
{\em geometry\_name}, are given according to a user-supplied description
in the ``ParFlow Binary'' file whose filename is given as the value.
The only choice for the value of {\em geometry\_name} is ``domain''.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.TensorByFileY   "perm_y.pfb"
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.Perm.TensorFileZ}{no default}
{
This key specifies that $k_z$ values for the specified geometry,
{\em geometry\_name}, are given according to a user-supplied description
in the ``ParFlow Binary'' file whose filename is given as the value.
The only choice for the value of {\em geometry\_name} is ``domain''.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Perm.TensorByFileZ   "perm_z.pfb"
\end{verbatim}\end{display}

%Here are examples of each of the above input types:
%\begin{display}\begin{verbatim}
%constant field:
%  17 0
%  0.83
%turning bands:
%  43 1
%  1.0 1.0 1.0  2.8 1.0
%  50 5.0 100.0 0.2 1 1
%  1
%parallel Gaussian simulator:
%  43 2
%  1.0 1.0 1.0  2.8 1.0
%  1 200 5 1
%  0
%input field:
%  15 3
%  17 
%  perm_geounit1.pfb
%\end{verbatim}\end{display}


%=============================================================================
%=============================================================================

\subsection{Porosity}
\label{Porosity}
\pfaddanchor{pfinPorosity}

Here, porosity values are assigned within geounits (specified in
\S~\ref{Geometries} above) using one of the methods described below.

The format for this section of input is:

\pfkey{list}{Geom.Porosity.GeomNames}{no default}
{
This key specifies all of the geometries on which a porosity will be 
assigned.  These geometries must cover the entire computational domain.
}
\begin{display}\begin{verbatim}
pfset Geom.Porosity.GeomNames   "background"
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.Porosity.Type}{no default}
{
This key specifies which method is to be used to assign porosity data
to the named geometry, {\em geometry\_name}.  The only choice currently
available is {\bf Constant} which indicates that a constant is to be
assigned to all grid cells within a geometry.
}
\begin{display}\begin{verbatim}
pfset Geom.background.Porosity.Type   Constant
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.Porosity.Value}{no default}
{
This key specifies the value assigned to all points in the named
geometry, {\em geometry\_name}, if the type was set to constant.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Porosity.Value   1.0
\end{verbatim}\end{display}

\subsection{Specific Storage}
\label{Specific Storage}
\pfaddanchor{pfinSpecStorage}

Here, specific storage ($S_s$ in Equation \ref{eq:richard}) values are assigned within geounits (specified in
\S~\ref{Geometries} above) using one of the methods described below.

The format for this section of input is:

\pfkey{list}{Specific Storage.GeomNames}{no default}
{
This key specifies all of the geometries on which a different specific storage value will be 
assigned.  These geometries must cover the entire computational domain.
}
\begin{display}\begin{verbatim}
pfset SpecificStorage.GeomNames       "domain"
\end{verbatim}\end{display}

\pfkey{string}{SpecificStorage.Type}{no default}
{
This key specifies which method is to be used to assign specific storage data.  The only choice currently
available is {\bf Constant} which indicates that a constant is to be
assigned to all grid cells within a geometry.
}
\begin{display}\begin{verbatim}
pfset SpecificStorage.Type            Constant
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.SpecificStorage.Value}{no default}
{
This key specifies the value assigned to all points in the named
geometry, {\em geometry\_name}, if the type was set to constant.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.SpecificStorage.Value 1.0e-4
\end{verbatim}\end{display}

%%
%% == @RMM dZ Multipliers
%%
\subsection{dZMultipliers}
\label{dZ Multipliers}
\pfaddanchor{pfindZMult}

Here, dZ multipliers ($\delta Z * m$) values are assigned within geounits (specified in
\S~\ref{Geometries} above) using one of the methods described below.

The format for this section of input is:

\pfkey{string}{ Solver.Nonlinear.VariableDz}{False}
{
This key specifies whether dZ multipliers are to be used, the default is False.
The default indicates a false or non-active variable dz and each layer thickness is 1.0 [L].
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.VariableDz     True
\end{verbatim}\end{display}

\pfkey{list}{dzScale.GeomNames}{no default}
{
This key specifies which problem domain is being applied a variable dz subsurface.
These geometries must cover the entire computational domain.
}
\begin{display}\begin{verbatim}
pfset dzScale.GeomNames domain
\end{verbatim}\end{display}


\pfkey{string}{dzScale.Type}{no default}
{
This key specifies which method is to be used to assign variable vertical grid spacing.  The choices currently
available are {\bf Constant} which indicates that a constant is to be
assigned to all grid cells within a geometry, {\bf nzList} which assigns all layers of a given model to a list value,
and {\bf PFBFile} which reads in values from a distributed pfb file.
}
\begin{display}\begin{verbatim}
pfset dzScale.Type            Constant
\end{verbatim}\end{display}

\pfkey{list}{Specific dzScale.GeomNames}{no default}
{
This key specifies all of the geometries on which a different dz scaling value will be 
assigned.  These geometries must cover the entire computational domain.
}
\begin{display}\begin{verbatim}
pfset dzScale.GeomNames       "domain"
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.dzScale.Value}{no default}
{
This key specifies the value assigned to all points in the named
geometry, {\em geometry\_name}, if the type was set to constant.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.dzScale.Value 1.0
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.dzScale.FileName}{no default}
{
This key specifies file to be read in for variable dz values for the given 
geometry, {\em geometry\_name}, if the type was set to {\bf PFBFile}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.dzScale.FileName  vardz.pfb
\end{verbatim}\end{display}

\pfkey{integer}{dzScale.nzListNumber}{no default}
{
This key indicates the number of layers with variable dz in the subsurface.
This value is the same as the \emph{ComputationalGrid.NZ} key.
}
\begin{display}\begin{verbatim}
pfset dzScale.nzListNumber  10
\end{verbatim}\end{display}


\pfkey{double}{Cell.{\em nzListNumber}.dzScale.Value}{no default}
{
This key assigns the thickness of each layer defined by nzListNumber.
ParFlow assigns the layers from the bottom-up (i.e. the bottom 
of the domain is layer 0, the top is layer NZ-1). The total domain depth
(\emph{Geom.domain.Upper.Z}) does not change with variable dz. The layer thickness
is calculated by \emph{ComputationalGrid.DZ *dZScale}.
}
\begin{display}\begin{verbatim}
pfset Cell.0.dzScale.Value 1.0
\end{verbatim}\end{display}

Example Usage: 
\begin{display}\begin{verbatim}

#--------------------------------------------
# Variable dz Assignments
#------------------------------------------
# Set VariableDz to be true
# Indicate number of layers (nzlistnumber), which is the same as nz
# (1) There is nz*dz = total depth to allocate,  
# (2) Each layer’s thickness is dz*dzScale, and
# (3) Assign the layer thickness from the bottom up.
# In this example nz = 5; dz = 10; total depth 40;
# Layers 	Thickness [m]
# 0 		15 			Bottom layer
# 1		15
# 2		5
# 3		4.5			
# 4 		0.5			Top layer
pfset Solver.Nonlinear.VariableDz     True
pfset dzScale.GeomNames            domain
pfset dzScale.Type            nzList
pfset dzScale.nzListNumber       5
pfset Cell.0.dzScale.Value 1.5
pfset Cell.1.dzScale.Value 1.5
pfset Cell.2.dzScale.Value 0.5
pfset Cell.3.dzScale.Value 0.45
pfset Cell.4.dzScale.Value 0.05
\end{verbatim}\end{display}



%%
%% == mannings roughness
%%
\subsection{Manning's Roughness Values}
\label{Manning's Roughness Values}
\pfaddanchor{pfinMannings}

Here, Manning's roughness values ($n$ in Equations \ref{eq:manningsx} and \ref{eq:manningsy}) are assigned to the upper boundary of the domain using one of the methods described below.

The format for this section of input is:

\pfkey{list}{Mannings.GeomNames}{no default}
{
This key specifies all of the geometries on which a different Mannings roughness value will be 
assigned.  Mannings values may be assigned by {\bf PFBFile} or as {\bf Constant} by geometry.  These geometries must cover the entire upper surface of the computational domain.
}
\begin{display}\begin{verbatim}
pfset Mannings.GeomNames       "domain"
\end{verbatim}\end{display}

\pfkey{string}{Mannings.Type}{no default}
{
This key specifies which method is to be used to assign Mannings roughness data.  The choices currently
available are {\bf Constant} which indicates that a constant is to be
assigned to all grid cells within a geometry and {\bf PFBFile} which indicates that all values are read in from a distributed, grid-based \parflow{} binary file.
}
\begin{display}\begin{verbatim}
pfset Mannings.Type "Constant"
\end{verbatim}\end{display}

\pfkey{double}{Mannings.Geom.{\em geometry\_name}.Value}{no default}
{
This key specifies the value assigned to all points in the named
geometry, {\em geometry\_name}, if the type was set to constant.
}
\begin{display}\begin{verbatim}
pfset Mannings.Geom.domain.Value 5.52e-6
\end{verbatim}\end{display}

\pfkey{double}{Mannings.FileName}{no default}
{
This key specifies the value assigned to all points be read in from a \parflow{} binary file.
}
\begin{display}\begin{verbatim}
pfset Mannings.FileName roughness.pfb
\end{verbatim}\end{display}

Complete example of setting Mannings roughness $n$ values by geometry:
\begin{display}\begin{verbatim}
pfset Mannings.Type "Constant"
pfset Mannings.GeomNames "domain"
pfset Mannings.Geom.domain.Value 5.52e-6
\end{verbatim}\end{display}

%%==
%%==
%%==
\subsection{Topographical Slopes}
\label{Topographical Slopes}
\pfaddanchor{pfintoposlopes}

Here, topographical slope values ($S_{f,x}$ and $S_{f,y}$ in Equations \ref{eq:manningsx} and \ref{eq:manningsy}) are assigned to the upper boundary of the domain using one of the methods described below. Note that due to the negative sign in these equations $S_{f,x}$ and $S_{f,y}$ take a sign in the direction \emph{opposite} of the direction of the slope.  That is, negative slopes point "downhill" and positive slopes "uphill".

The format for this section of input is:

\pfkey{list}{ToposlopesX.GeomNames}{no default}
{
This key specifies all of the geometries on which a different $x$ topographic slope values will be 
assigned.  Topographic slopes may be assigned by {\bf PFBFile} or as {\bf Constant} by geometry.  These geometries must cover the entire upper surface of the computational domain.
}
\begin{display}\begin{verbatim}
pfset ToposlopesX.GeomNames       "domain"
\end{verbatim}\end{display}

\pfkey{list}{ToposlopesY.GeomNames}{no default}
{
This key specifies all of the geometries on which a different $y$ topographic slope values will be 
assigned.  Topographic slopes may be assigned by {\bf PFBFile} or as {\bf Constant} by geometry.  These geometries must cover the entire upper surface of the computational domain.
}
\begin{display}\begin{verbatim}
pfset ToposlopesY.GeomNames       "domain"
\end{verbatim}\end{display}

\pfkey{string}{ToposlopesX.Type}{no default}
{
This key specifies which method is to be used to assign topographic slopes.  The choices currently
available are {\bf Constant} which indicates that a constant is to be
assigned to all grid cells within a geometry and {\bf PFBFile} which indicates that all values are read in from a distributed, grid-based \parflow{} binary file.
}
\begin{display}\begin{verbatim}
pfset ToposlopesX.Type "Constant"
\end{verbatim}\end{display}

\pfkey{double}{ToposlopeX.Geom.{\em geometry\_name}.Value}{no default}
{
This key specifies the value assigned to all points in the named
geometry, {\em geometry\_name}, if the type was set to constant.
}
\begin{display}\begin{verbatim}
pfset ToposlopeX.Geom.domain.Value 0.001
\end{verbatim}\end{display}

\pfkey{double}{ToposlopesX.FileName}{no default}
{
This key specifies the value assigned to all points be read in from a \parflow{} binary file.
}
\begin{display}\begin{verbatim}
pfset TopoSlopesX.FileName lw.1km.slope_x.pfb
\end{verbatim}\end{display}

\pfkey{double}{ToposlopesY.FileName}{no default}
{
This key specifies the value assigned to all points be read in from a \parflow{} binary file.
}
\begin{display}\begin{verbatim}
pfset TopoSlopesY.FileName lw.1km.slope_y.pfb
\end{verbatim}\end{display}

Example of setting $x$ and $y$ slopes by geometry:
\begin{display}\begin{verbatim}
pfset TopoSlopesX.Type "Constant"
pfset TopoSlopesX.GeomNames "domain"
pfset TopoSlopesX.Geom.domain.Value 0.001

pfset TopoSlopesY.Type "Constant"
pfset TopoSlopesY.GeomNames "domain"
pfset TopoSlopesY.Geom.domain.Value -0.001
\end{verbatim}\end{display}

Example of setting $x$ and $y$ slopes by file:
\begin{display}\begin{verbatim}
pfset TopoSlopesX.Type "PFBFile"
pfset TopoSlopesX.GeomNames "domain"
pfset TopoSlopesX.FileName lw.1km.slope_x.pfb

pfset TopoSlopesY.Type "PFBFile"
pfset TopoSlopesY.GeomNames "domain"
pfset TopoSlopesY.FileName lw.1km.slope_y.pfb
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Retardation}
\label{Retardation}
\pfaddanchor{pfinRetardation}

Here, retardation values are assigned for contaminants within
geounits (specified in \S~\ref{Geometries} above) using one of
the functions described below.  The format for this section of
input is:

\pfkey{list}{Geom.Retardation.GeomNames}{no default}
{
This key specifies all of the geometries to which the contaminants will
have a retardation function applied.
}
\begin{display}\begin{verbatim}
pfset GeomInput.Names   "background"
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geometry\_name}.{\em contaminant\_name}.Retardation.Type}{no default}
{
This key specifies which function is to be used to compute the
retardation for the named contaminant, {\em contaminant\_name}, in the
named geometry, {\em geometry\_name}.  The only choice currently
available is {\bf Linear} which indicates that a simple linear
retardation function is to be used to compute the retardation.
}
\begin{display}\begin{verbatim}
pfset Geom.background.tce.Retardation.Type   Linear
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geometry\_name}.{\em contaminant\_name}.Retardation.Value}{no default}
{
This key specifies the distribution coefficient for the linear function
used to compute the retardation of the named contaminant,
{\em contaminant\_name}, in the named geometry, {\em geometry\_name}.
The value should be scaled by the density of the material in the
geometry.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Retardation.Value   0.2
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Full Multiphase Mobilities}

Here we define phase mobilities by specifying the
relative permeability function.  Input is specified differently depending on
what problem is being specified.  For full multi-phase problems, the following
input keys are used.  See the next section for the correct Richards' equation 
input format.

\pfkey{string}{Phase.{\em phase\_name}.Mobility.Type}{no default}
{
This key specifies whether the mobility for {\em phase\_name} will be a given
constant or a polynomial of the form,  
$(S - S_0)^{a}$, where $S$ is saturation, $S_0$ is
irreducible saturation, and $a$ is some exponent.
The possibilities for this key are {\bf Constant} and {\bf Polynomial}.
}
\begin{display}\begin{verbatim}
pfset Phase.water.Mobility.Type   Constant
\end{verbatim}\end{display}

\pfkey{double}{Phase.{\em phase\_name}.Mobility.Value}{no default}
{
This key specifies the constant mobility value for phase {\em phase\_name}.
}
\begin{display}\begin{verbatim}
pfset Phase.water.Mobility.Value   1.0
\end{verbatim}\end{display}

\pfkey{double}{Phase.{\em phase\_name}.Mobility.Exponent}{2.0}
{
This key specifies the exponent used in a polynomial representation of the
relative permeability.  Currently, only a value of $2.0$ is allowed for this
key. 
}
\begin{display}\begin{verbatim}
pfset Phase.water.Mobility.Exponent   2.0
\end{verbatim}\end{display}

\pfkey{double}{Phase.{\em phase\_name}.Mobility.IrreducibleSaturation}
{0.0}
{
This key specifies the irreducible saturation 
used in a polynomial representation of the relative permeability.
Currently, only a value of 0.0 is allowed for this key.
}
\begin{display}\begin{verbatim}
pfset Phase.water.Mobility.IrreducibleSaturation   0.0
\end{verbatim}\end{display}


\subsection{Richards' Equation Relative Permeabilities}
\label{Richards RelPerm}

The following keys are used to describe relative permeability input for the
Richards' equation implementation.  They will be ignored if a full two-phase
formulation is used.

\pfkey{string}{Phase.RelPerm.Type}{no default}
{
This key specifies the type of relative permeability function that will be used
on all specified geometries.  Note that only one type of relative permeability
may be used for the entire problem.  However, parameters may be different for
that type in different geometries.  For instance, if the problem consists of
three geometries, then {\bf VanGenuchten} may be specified with three different
sets of parameters for the three different goemetries.  However, once {\bf
VanGenuchten} is specified, one geometry cannot later be specified to have {\bf
Data} as its relative permeability.  The possible values for this key 
are {\bf Constant, VanGenuchten, Haverkamp, Data,} and {\bf  Polynomial}.
}
\begin{display}\begin{verbatim}
pfset Phase.RelPerm.Type   Constant
\end{verbatim}\end{display}

The various possible functions are defined as follows.
The {\bf Constant} specification means that the relative permeability will be
constant on the specified geounit.  The {\bf VanGenuchten} specification means
that the relative permeability will be given as a Van Genuchten function
\cite{VanGenuchten80} with the form,
\begin{eqnarray}
k_r(p) = \frac{(1 - \frac{(\alpha p)^{n-1}}{(1 + (\alpha p)^n)^m})^2}
{(1 + (\alpha p)^n)^{m/2}},
\end{eqnarray}
where $\alpha$ and $n$ are soil parameters and $m = 1 - 1/n$, on each region.
The {\bf Haverkamp} specification means that the relative permeability will be
given in the following form \cite{Haverkamp-Vauclin81},
\begin{eqnarray}
k_r(p) = \frac{A}{A + p^{\gamma}},
\end{eqnarray}
where $A$ and $\gamma$ are soil parameters, on each region.
The {\bf Data} specification is currently unsupported but will later mean that
data points for the relative permeability curve will be given and \parflow{}
will set up the proper interpolation coefficients to get values between the
given data points.
The {\bf Polynomial} specification
defines a polynomial relative permeability function for each 
region of the form,
\begin{eqnarray}
k_r(p) = \sum_{i=0}^{degree} c_ip^i.
\end{eqnarray}

\pfkey{list}{Phase.RelPerm.GeomNames}{no default}
{This key specifies the geometries on which relative permeability will be
given.  The union of these geometries must cover the entire computational
domain.}
\begin{display}\begin{verbatim}
pfset Phase.RelPerm.Geonames   domain
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.RelPerm.Value}{no default}
{This key specifies the constant relative permeability value on the specified
geometry. }
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.Value    0.5
\end{verbatim}\end{display}

\pfkey{integer}{Phase.RelPerm.VanGenuchten.File}{0}
{This key specifies whether soil parameters for the VanGenuchten function are
specified in a pfb file or by region.  The options are either 0 for
specification by region, or 1 for specification in a file.  Note that either
all parameters are specified in files (each has their own input file) or none
are specified by files.  Parameters specified by files are: $\alpha$ and N.}
\begin{display}\begin{verbatim}
pfset Phase.RelPerm.VanGenuchten.File   1
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.RelPerm.Alpha.Filename}{no default}
{This key specifies a pfb filename containing the alpha parameters for the
VanGenuchten function cell-by-cell.  The ONLY option for {\em geom\_name} is
``domain''.}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.Alpha.Filename   alphas.pfb
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.RelPerm.N.Filename}{no default}
{This key specifies a pfb filename containing the N parameters for the
VanGenuchten function cell-by-cell.  The ONLY option for {\em geom\_name} is
``domain''.}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.N.Filename   Ns.pfb
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.RelPerm.Alpha}{no default}
{This key specifies the $\alpha$ parameter for the Van Genuchten function
specified on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.Alpha  0.005
\end{verbatim}\end{display}


\pfkey{double}{Geom.{\em geom\_name}.RelPerm.N}{no default}
{This key specifies the $N$ parameter for the Van Genuchten function specified
on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.N   2.0
\end{verbatim}\end{display}

\pfkey{int}{Geom.{\em geom\_name}.RelPerm.NumSamplePoints}{0}
{This key specifies the number of sample points for a spline base interpolation table
for the Van Genuchten function specified on {\em geom\_name}.   If this number is 0 (the default) 
then the function is evaluated directly.   Using the interpolation table is faster but is less
accurate.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.NumSamplePoints  20000
\end{verbatim}\end{display}

\pfkey{int}{Geom.{\em geom\_name}.RelPerm.MinPressureHead}{no default}
 {This key specifies the lower value for a spline base interpolation
 table for the Van Genuchten function specified on {\em geom\_name}.
  The upper value of the range is 0.  This value is used only when the
 table lookup method is used ({\em NumSamplePoints} is greater than 0).
  }
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.MinPressureHead -300
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.RelPerm.A}{no default}
{This key specifies the $A$ parameter for the Haverkamp relative permeability
on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.A  1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.RelPerm.Gamma}{no default}
{This key specifies the the $\gamma$ parameter for the Haverkamp relative
permeability on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.Gamma  1.0
\end{verbatim}\end{display}

\pfkey{integer}{Geom.{\em geom\_name}.RelPerm.Degree}{no default}
{This key specifies the degree of the polynomial for the Polynomial relative
permeability given on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.Degree  1
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.RelPerm.Coeff.{\em coeff\_number}}
{no default}
{This key specifies the {\em coeff\_number}th coefficient of the Polynomial
relative permeability given on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.RelPerm.Coeff.0  0.5
pfset Geom.domain.RelPerm.Coeff.1  1.0
\end{verbatim}\end{display}

NOTE: For all these cases, if only one region is to be used (the domain),
the background region should NOT be set as that single region.  
Using the background will prevent the
upstream weighting from being correct near Dirichlet boundaries.

%=============================================================================
%=============================================================================

\subsection{Phase Sources}
\label{Phase Sources}

The following keys are used to specify phase source terms.
The units of the source term are $1/T$.
So, for example, to specify a region with constant flux rate of $L^3/T$,
one must be careful to convert this rate to the proper units by
dividing by the volume of the enclosing region.
For {\em Richards' equation} input, the source term must be given as a flux
multiplied by density.

\pfkey{string}{PhaseSources.{\em phase\_name}.Type}
{no default}
{This key specifies the type of source to use for phase {\em phase\_name}.  
Possible values for this key are {\bf Constant} and {\bf PredefinedFunction}.
{\bf Constant} type phase sources specify a constant phase source 
value for a given set of regions.  
{\bf PredefinedFunction} type phase sources use a preset function
(choices are listed below) to specify the source.  Note that the
{\bf PredefinedFunction} type can only be used to set a single source 
over the entire domain and not separate sources over different regions.
}
\begin{display}\begin{verbatim}
pfset PhaseSources.water.Type   Constant
\end{verbatim}\end{display}


\pfkey{list}{PhaseSources.{\em phase\_name}.GeomNames}{no default}
{This key specifies the names of the geometries on which source terms will be
specified.  This is used only for {\bf Constant} type phase sources.
Regions listed later ``overlay'' regions listed earlier.
}
\begin{display}\begin{verbatim}
pfset PhaseSources.water.GeomNames   "bottomlayer middlelayer toplayer"
\end{verbatim}\end{display}

\pfkey{double}{PhaseSources.{\em phase\_name}.Geom.{\em geom\_name}.Value}
{no default}
{This key specifies the value of a constant source term applied to phase
{\em phase \_name} on geometry {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset PhaseSources.water.Geom.toplayer.Value   1.0
\end{verbatim}\end{display}


\pfkey{string}{PhaseSources.{\em phase\_name}.PredefinedFunction}{no default}
{This key specifies which of the predefined functions will be used for the
source.
Possible values for this key are {\bf X, XPlusYPlusZ, X3Y2PlusSinXYPlus1,}
\newline {\bf X3Y4PlusX2PlusSinXYCosYPlus1, XYZTPlus1} and {\bf
XYZTPlus1PermTensor}.
}
\begin{display}\begin{verbatim}
pfset PhaseSources.water.PredefinedFunction   XPlusYPlusZ
\end{verbatim}\end{display}

The choices for this key correspond to sources as follows:
\begin{description}
\item[{\bf X}: ] ${\rm source}\; = 0.0$

\item[{\bf XPlusYPlusX}: ] ${\rm source}\; = 0.0$

\item[{\bf X3Y2PlusSinXYPlus1}:]
${\rm source}\; = -(3x^2 y^2 + y\cos(xy))^2 - (2x^3 y + x\cos(xy))^2 
- (x^3 y^2 + \sin(xy) + 1) (6x y^2 + 2x^3 -(x^2 +y^2) \sin(xy))$ \\
This function type specifies that the source applied over the entire domain is
as noted above.  This corresponds to $p=x^{3}y^{2}+\sin(xy)+1$ in the problem 
$-\nabla\cdot (p\nabla p)=f$.

\item[{\bf X3Y4PlusX2PlusSinXYCosYPlus1}:] 
${\rm source}\; = -(3x^22 y^4 + 2x + y\cos(xy)\cos(y))^2 
- (4x^3 y^3 + x\cos(xy)\cos(y) - \sin(xy)\sin(y))^2 
- (x^3 y^4 + x^2 + \sin(xy)\cos(y) + 1)
(6xy^4 + 2 - (x^2 + y^2 + 1)\sin(xy)\cos(y) 
+ 12x^3 y^2 - 2x\cos(xy)\sin(y))$ \\
This function type specifies that the source applied over the entire domain is
as noted above.  This corresponds to $p=x^{3}y^{4}+x^{2}+\sin (xy)\cos(y) +1$ 
in the problem $-\nabla\cdot (p\nabla p)=f$.

\item[{\bf XYZTPlus1}: ]
${\rm source}\; = xyz - t^2 (x^2 y^2 +x^2 z^2 +y^2 z^2)$ \\
This function type specifies that the source applied over the entire domain is
as noted above.  This corresponds to $p = xyzt + 1$ in the problem 
$\frac{\partial p}{\partial t}-\nabla\cdot (p\nabla p)=f$.

\item[{\bf XYZTPlus1PermTensor}: ]
${\rm source}\; = xyz - t^2 (x^2 y^2 3 + x^2 z^2 2 + y^2 z^2)$ \\
This function type specifies that the source applied over the entire domain is
as noted above.  This corresponds to $p = xyzt + 1$ in the problem 
$\frac{\partial p}{\partial t}-\nabla\cdot (Kp\nabla p)=f$, where 
$K = diag(1 \;\; 2 \;\; 3)$.

\end{description}


%=============================================================================
%=============================================================================
\subsection{Capillary Pressures}
\label{Capillary Pressures}
\pfaddanchor{pfinCapillaryPressures}

Here we define capillary pressure.
Note: this section needs to be defined {\em only} for multi-phase flow
and should not be defined for single phase and Richards' equation cases.
The format for this section of input is:

\pfkey{string}{CapPressure.{\em phase\_name}.Type}{"Constant"}
{
This key specifies the capillary pressure between phase $0$ and the named
phase, {\em phase\_name}. The only choice available is {\bf Constant}
which indicates that a constant capillary pressure exists between the
phases.
}
\begin{display}\begin{verbatim}
pfset CapPressure.water.Type   Constant
\end{verbatim}\end{display}

\pfkey{list}{CapPressure.{\em phase\_name}.GeomNames}{no default}
{
This key specifies the geometries that capillary pressures will be
computed for in the named phase, {\em phase\_name}.  Regions listed
later ``overlay'' regions listed earlier.  Any geometries not listed
will be assigned $0.0$ capillary pressure by \parflow{}.
}
\begin{display}\begin{verbatim}
pfset CapPressure.water.GeomNames   "domain"
\end{verbatim}\end{display}


\pfkey{double}{Geom.{\em geometry\_name}.CapPressure.{\em phase\_name}.Value}
{0.0}
{
This key specifies the value of the capillary pressure in the named
geometry, {\em geometry\_name}, for the named phase, {\em phase\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.CapPressure.water.Value   0.0
\end{verbatim}\end{display}


{\sl Important note}: the code currently works only for capillary
pressure equal zero.

%=============================================================================
%=============================================================================

%=============================================================================
%=============================================================================

\subsection{Saturation}
\label{Saturation}

This section is {\em only} relevant to the Richards' equation cases.  All keys
relating to this section will be ignored for other cases.
The following keys are used to define the saturation-pressure curve.

\pfkey{string}{Phase.Saturation.Type}{no default}
{
This key specifies the type of saturation function that will be used
on all specified geometries.  Note that only one type of saturation
may be used for the entire problem.  However, parameters may be different for
that type in different geometries.  For instance, if the problem consists of
three geometries, then {\bf VanGenuchten} may be specified with three different
sets of parameters for the three different goemetries.  However, once {\bf
VanGenuchten} is specified, one geometry cannot later be specified to have {\bf
Data} as its saturation.  The possible values for this key 
are {\bf Constant, VanGenuchten, Haverkamp, Data, Polynomial} and 
{\bf PFBFile}.
}
\begin{display}\begin{verbatim}
pfset Phase.Saturation.Type   Constant
\end{verbatim}\end{display}

The various possible functions are defined as follows.
The {\bf Constant} specification means that the saturation will be
constant on the specified geounit.  The {\bf VanGenuchten} specification means
that the saturation will be given as a Van Genuchten function 
\cite{VanGenuchten80} with the form,
\begin{eqnarray}
s(p) = \frac{s_{sat} - s_{res}}{(1 + (\alpha p)^n)^m} + s_{res},
\end{eqnarray}
where $s_{sat}$ is the saturation at saturated conditions, 
$s_{res}$ is the residual saturation, and
$\alpha$ and $n$ are soil parameters with $m = 1 - 1/n$, on each region.
The {\bf Haverkamp} specification means that the saturation will be
given in the following form \cite{Haverkamp-Vauclin81},
\begin{eqnarray}
s(p) = \frac{\alpha(s_{sat} - s_{res})}{A + p^{\gamma}} + s_{res},
\end{eqnarray}
where $A$ and $\gamma$ are soil parameters, on each region.
The {\bf Data} specification is currently unsupported but will later mean that
data points for the saturation curve will be given and \parflow{}
will set up the proper interpolation coefficients to get values between the
given data points.
The {\bf Polynomial} specification
defines a polynomial saturation function for each region of the form,
\begin{eqnarray}
s(p) = \sum_{i=0}^{degree} c_ip^i.
\end{eqnarray}
The {\bf PFBFile} specification means that the saturation will be taken as a
spatially varying but constant in pressure function given by data in a
\parflow{} binary (.pfb) file.  

\pfkey{list}{Phase.Saturation.GeomNames}{no default}
{This key specifies the geometries on which saturation will be given.
The union of these geometries must cover the entire computational domain.}
\begin{display}\begin{verbatim}
pfset Phase.Saturation.Geonames   domain
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.Saturation.Value}{no default}
{This key specifies the constant saturation value on the {\em geom\_name} 
region. }
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.Value    0.5
\end{verbatim}\end{display}

\pfkey{integer}{Phase.Saturation.VanGenuchten.File}{0}
{This key specifies whether soil parameters for the VanGenuchten function are
specified in a pfb file or by region.  The options are either 0 for
specification by region, or 1 for specification in a file.  Note that either
all parameters are specified in files (each has their own input file) or none
are specified by files.  Parameters specified by files are $\alpha$, N, SRes,
and SSat.}
\begin{display}\begin{verbatim}
pfset Phase.Saturation.VanGenuchten.File   1
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.Saturation.Alpha.Filename}{no default}
{This key specifies a pfb filename containing the alpha parameters for the
VanGenuchten function cell-by-cell.  The ONLY option for {\em geom\_name} is
``domain''.}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.Filename   alphas.pfb
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.Saturation.N.Filename}{no default}
{This key specifies a pfb filename containing the N parameters for the
VanGenuchten function cell-by-cell.  The ONLY option for {\em geom\_name} is
``domain''.}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.N.Filename   Ns.pfb
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.Saturation.SRes.Filename}{no default}
{This key specifies a pfb filename containing the SRes parameters for the
VanGenuchten function cell-by-cell.  The ONLY option for {\em geom\_name} is
``domain''.}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.SRes.Filename   SRess.pfb
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.Saturation.SSat.Filename}{no default}
{This key specifies a pfb filename containing the SSat parameters for the
VanGenuchten function cell-by-cell.  The ONLY option for {\em geom\_name} is
``domain''.}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.SSat.Filename   SSats.pfb
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.Saturation.Alpha}{no default}
{This key specifies the $\alpha$ parameter for the Van Genuchten function
specified on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.Alpha  0.005
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.Saturation.N}{no default}
{This key specifies the $N$ parameter for the Van Genuchten function specified
on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.N   2.0
\end{verbatim}\end{display}

Note that if both a Van Genuchten saturation and relative permeability are
specified, then the soil parameters should be the same for each in order
to have a consistent problem.

\pfkey{double}{Geom.{\em geom\_name}.Saturation.SRes}{no default}
{This key specifies the residual saturation on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.SRes   0.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.Saturation.SSat}{no default}
{This key specifies the saturation at saturated conditions on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.SSat   1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.Saturation.A}{no default}
{This key specifies the $A$ parameter for the Haverkamp saturation
on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.A   1.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.Saturation.Gamma}{no default}
{This key specifies the the $\gamma$ parameter for the Haverkamp saturation 
on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.Gamma   1.0
\end{verbatim}\end{display}

\pfkey{integer}{Geom.{\em geom\_name}.Saturation.Degree}{no default}
{This key specifies the degree of the polynomial for the Polynomial 
saturation given on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.Degree   1
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.Saturation.Coeff.{\em coeff\_number}}
{no default}
{This key specifies the {\em coeff\_number}th coefficient of the Polynomial
saturation given on {\em geom\_name}.
}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.Coeff.0   0.5
pfset Geom.domain.Saturation.Coeff.1   1.0
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.Saturation.FileName}{no default}
{This key specifies the name of the file containing saturation values for the
domain.  It is assumed that {\em geom\_name} is ``domain'' for this key.}
\begin{display}\begin{verbatim}
pfset Geom.domain.Saturation.FileName  "domain_sats.pfb"
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Internal Boundary Conditions}
\label{Internal Boundary Conditions}
\pfaddanchor{pfinInternalBoundaryConditions}

In this section, we define internal Dirichlet boundary conditions
by setting the pressure at points in the domain.
The format for this section of input is:

\pfkey{string}{InternalBC.Names}{no default}
{
This key specifies the names for the internal boundary conditions.
At each named point, ${\rm x}$, ${\rm y}$ and ${\rm z}$ will specify
the coordinate locations and ${\rm h}$ will specify the hydraulic head
value of the condition.  This real location is ``snapped'' to the
nearest gridpoint in \parflow{}.

NOTE: Currently, \parflow{} assumes that internal boundary conditions
and pressure wells are separated by at least one cell from any external
boundary.  The user should be careful of this when defining the input
file and grid.
}
\begin{display}\begin{verbatim}
pfset InternalBC.Names   "fixedvalue"
\end{verbatim}\end{display}

\pfkey{double}{InternalBC.{\em internal\_bc\_name}.X}{no default}
{
This key specifies the x-coordinate, ${\rm x}$, of the named,
{\em internal\_bc\_name}, condition.
}
\begin{display}\begin{verbatim}
pfset InternalBC.fixedheadvalue.X   40.0
\end{verbatim}\end{display}

\pfkey{double}{InternalBC.{\em internal\_bc\_name}.Y}{no default}
{
This key specifies the y-coordinate, ${\rm y}$, of the named,
{\em internal\_bc\_name}, condition.
}
\begin{display}\begin{verbatim}
pfset InternalBC.fixedheadvalue.Y   65.2
\end{verbatim}\end{display}

\pfkey{double}{InternalBC.{\em internal\_bc\_name}.Z}{no default}
{
This key specifies the z-coordinate, ${\rm z}$, of the named,
{\em internal\_bc\_name}, condition.
}
\begin{display}\begin{verbatim}
pfset InternalBC.fixedheadvalue.Z   12.1
\end{verbatim}\end{display}

\pfkey{double}{InternalBC.{\em internal\_bc\_name}.Value}{no default}
{
This key specifies the value of the named,
{\em internal\_bc\_name}, condition.
}
\begin{display}\begin{verbatim}
pfset InternalBC.fixedheadvalue.Value   100.0
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Boundary Conditions: Pressure}
\label{Boundary Conditions: Pressure}

Here we define the pressure boundary conditions.
The Dirichlet conditions below are hydrostatic conditions, and it
is assumed that at each phase interface the pressure is constant.
{\em It is also assumed here that all phases
are distributed within the domain at all times such that the lighter phases 
are vertically higher than the heavier phases.}

Boundary condition input is associated with domain patches
(see \S~\ref{Domain}).  Note that different patches may have different types of
boundary conditions on them.


\pfkey{list}{BCPressure.PatchNames}{no default}
{This key specifies the names of patches on which pressure boundary conditions
will be specified.  Note that these must all be patches on the external
boundary of the domain and these patches must ``cover'' that external boundary.
}
\begin{display}\begin{verbatim}
pfset BCPressure.PatchNames    "left right front back top bottom"
\end{verbatim}\end{display}

\pfkey{string}{Patch.{\em patch\_name}.BCPressure.Type}{no default}
{This key specifies the type of boundary condition data given for patch 
{\em patch\_name}.  Possible values for this key are {\bf DirEquilRefPatch, 
DirEquilPLinear, FluxConst, FluxVolumetric, PressureFile, FluxFile, OverlandFow, OverlandFlowPFB} and 
{\bf ExactSolution}.  The choice {\bf DirEquilRefPatch} specifies that the
pressure on the specified patch will be in hydrostatic equilibrium with a
constant reference pressure given on a reference patch.
The choice {\bf DirEquilPLinear} specifies that the pressure on the specified
patch will be in hydrostatic equilibrium with pressure given along a piecewise
line at elevation $z=0$.  The choice {\bf FluxConst} defines a constant normal 
flux boundary condition through the domain patch.  
This flux must be specified in units of $[L]/[T]$.
For {\em Richards' equation}, fluxes must be specified as a mass flux and given
as the above flux multiplied by the density.  Thus, this choice of input type
for a Richards' equation problem has units of $([L]/[T])([M]/[L]^3)$.
The choice {\bf FluxVolumetric} defines a volumetric flux boundary condition
through the domain patch.  The units should be consistent with all other
user input for the problem.
For {\em Richards' equation} fluxes must be specified as a mass flux and given
as the above flux multiplied by the density.
The choice {\bf PressureFile} defines a hydraulic head boundary condition that
is read from a properly distributed .pfb file.  Only the values needed for the
patch are used.
The choice {\bf FluxFile} defines a flux boundary condition that is read form a
properly distributed .pfb file defined on a grid consistent with the pressure
field grid.  Only the values needed for the patch are used.
The choices {\bf OverlandFlow} and {\bf OverlandFlowPFB} both turn on 
fully-coupled overland flow routing as described in \cite{KM06} and in \S~\ref{Overland Flow}.  
The key {\bf OverlandFlow} corresponds to a {\bf Value} key with a positive or negative value, 
to indicate uniform fluxes (such as rainfall or evapotranspiration) over the entire domain 
while the key {\bf OverlandFlowPFB} allows a \code{.pfb} file to contain grid-based, 
spatially-variable fluxes.
The choice {\bf ExactSolution} specifies that an exact known solution 
is to be applied as a
Dirichlet boundary condition on the respective patch.  Note that this does not
change according to any cycle.  Instead, time dependence is handled by
evaluating at the time the boundary condition value is desired.  
The solution is specified by using a predefined function (choices are described
below).  NOTE: These last three types of boundary condition input is for 
{\em Richards' equation cases only!} 
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.Type  DirEquilRefPatch
\end{verbatim}\end{display}

\pfkey{string}{Patch.{\em patch\_name}.BCPressure.Cycle}{no default}
{This key specifies the time cycle to which boundary condition data for patch
{\em patch\_name} corresponds. 
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.Cycle   Constant
\end{verbatim}\end{display}

\pfkey{string}{Patch.{\em patch\_name}.BCPressure.RefGeom}{no default}
{This key specifies the name of the solid on which the reference patch for the
{\bf DirEquilRefPatch} boundary condition data is given.
Care should be taken to
make sure the correct solid is specified in cases of layered domains.  
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.RefGeom   domain
\end{verbatim}\end{display}

\pfkey{string}{Patch.{\em patch\_name}.BCPressure.RefPatch}{no default}
{This key specifies the reference patch on which the 
{\bf DirEquilRefPatch} boundary condition data is given.  This patch must be on
the reference solid specified by the Patch.{\em patch\_name}.BCPressure.RefGeom
key.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.RefPatch    bottom
\end{verbatim}\end{display}

\pfkey{double}{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.Value}
{no default}
{This key specifies the reference pressure value for the 
{\bf DirEquilRefPatch} boundary condition or the constant flux value for the
{\bf FluxConst} boundary condition, or the constant volumetric flux for the
{\bf FluxVolumetric} boundary condition.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.Value  -14.0
\end{verbatim}\end{display}

\pfkey{double}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.{\em phase\_name}.IntValue}
{no default}
{Note that the reference conditions for types {\bf DirEquilPLinear} and 
{\bf DirEquilRefPatch}
boundary conditions are for phase 0 {\em only}.  
This key specifies the constant pressure value
along the interface with phase {\em phase\_name} for cases with two phases
present.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.water.IntValue   -13.0
\end{verbatim}\end{display}

\pfkey{double}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.XLower}
{no default}
{This key specifies the lower $x$ coordinate of a line in the xy-plane.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.XLower  0.0
\end{verbatim}\end{display}

\pfkey{double}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.YLower}
{no default}
{This key specifies the lower $y$ coordinate of a line in the xy-plane.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.YLower  0.0
\end{verbatim}\end{display}

\pfkey{double}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.XUpper}
{no default}
{This key specifies the upper $x$ coordinate of a line in the xy-plane.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.XUpper  1.0
\end{verbatim}\end{display}

\pfkey{double}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.YUpper}
{no default}
{This key specifies the upper $y$ coordinate of a line in the xy-plane.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.YUpper  1.0
\end{verbatim}\end{display}

\pfkey{integer}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.NumPoints}
{no default}
{This key specifies the number of points on which pressure data is given along
the line used in the type {\bf DirEquilPLinear} boundary conditions.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.NumPoints   2
\end{verbatim}\end{display}

\pfkey{double}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.{\em point\_number}.Location}
{no default}
{This key specifies a number between 0 and 1 which represents
the location of a point on the line on which data is given for 
type {\bf DirEquilPLinear} boundary conditions.
Here 0 corresponds to the lower end of the line, and 1 corresponds to
the upper end.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.0.Location   0.0
\end{verbatim}\end{display}

\pfkey{double}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.{\em point\_number}.Value}
{no default}
{This key specifies the pressure value for phase 0 at point number 
{\em point\_number} and $z=0$ for
type {\bf DirEquilPLinear} boundary conditions.
All pressure values on the patch are determined by first projecting the
boundary condition coordinate onto the line, then linearly interpolating
between the neighboring point pressure values on the line.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.0.Value   14.0
\end{verbatim}\end{display}

\pfkey{string}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.FileName}
{no default}
{This key specifies the name of a properly distributed \code{.pfb} file 
that contains boundary data to be read for types {\bf PressureFile} 
and {\bf FluxFile}.  For flux data, the data must be defined over a grid 
consistent with the pressure field.  In both cases, only the values needed 
for the patch will be used.  The rest of the data is ignored.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.FileName   ocwd_bc.pfb
\end{verbatim}\end{display}



\pfkey{string}
{Patch.{\em patch\_name}.BCPressure.{\em interval\_name}.PredefinedFunction}
{no default}
{This key specifies the predefined function that will be used to specify
Dirichlet boundary conditions on patch {\em patch\_name}.
Note that this does not change according to any cycle.
Instead, time dependence is handled by
evaluating at the time the boundary condition value is desired.
Choices for this key include {\bf X, XPlusYPlusZ, X3Y2PlusSinXYPlus1, 
X3Y4PlusX2PlusSinXYCosYPlus1, XYZTPlus1} and {\bf XYZTPlus1PermTensor}.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCPressure.alltime.PredefinedFunction  XPlusYPlusZ
\end{verbatim}\end{display}
The choices for this key correspond to pressures as follows.
\begin{description}
\item[{\bf X}: ] $p = x$
\item[{\bf XPlusYPlusZ}: ] $p = x + y + z$
\item[{\bf X3Y2PlusSinXYPlus1}: ] $p = x^3 y^2 + \sin(xy) + 1$
\item[{\bf X3Y4PlusX2PlusSinXYCosYPlus1}: ] 
$p = x^3 y^4 + x^2 + \sin(xy)\cos y + 1$
\item[{\bf XYZTPlus1}: ] $p = xyzt + 1$
\item[{\bf XYZTPlus1PermTensor}: ] $p = xyzt + 1$
\end{description}

Example Script:
\begin{display}\begin{verbatim}

#---------------------------------------------------------
# Initial conditions: water pressure [m]
#---------------------------------------------------------
# Using a patch is great when you are not using a box domain
# If using a box domain HydroStaticDepth is fine
# If your RefPatch is z-lower (bottom of domain), the pressure is positive.
# If your RefPatch is z-upper (top of domain), the pressure is negative.
### Set water table to be at the bottom of the domain, the top layer is initially dry
pfset ICPressure.Type				HydroStaticPatch
pfset ICPressure.GeomNames		domain
pfset Geom.domain.ICPressure.Value	2.2

pfset Geom.domain.ICPressure.RefGeom	domain
pfset Geom.domain.ICPressure.RefPatch	z-lower

### Using a .pfb to initialize
pfset ICPressure.Type                                   PFBFile
pfset ICPressure.GeomNames		 "domain"
pfset Geom.domain.ICPressure.FileName	press.00090.pfb

pfset Geom.domain.ICPressure.RefGeom	domain
pfset Geom.domain.ICPressure.RefPatch	z-upper
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Boundary Conditions: Saturation}
\label{Boundary Conditions: Saturation}
\pfaddanchor{pfinBoundaryConditionsSaturation}

Note: this section needs to be defined {\em only} for multi-phase flow
and should {\em not} be defined for the single phase and Richards'
equation cases. 

Here we define the boundary conditions for the saturations.  Boundary
condition input is associated with domain patches (see \S~\ref{Domain}).
Note that different patches may have different types of boundary
conditions on them.

\pfkey{list}{BCSaturation.PatchNames}{no default}
{
This key specifies the names of patches on which saturation boundary
conditions will be specified.  Note that these must all be patches on
the external boundary of the domain and these patches must ``cover''
that external boundary.
}
\begin{display}\begin{verbatim}
pfset BCSaturation.PatchNames    "left right front back top bottom"
\end{verbatim}\end{display}

\pfkey{string}{Patch.{\em patch\_name}.BCSaturation.{\em phase\_name}.Type}{no default}
{
This key specifies the type of boundary condition data given for the
given phase, {\em phase\_name}, on the given patch {\em patch\_name}.
Possible values for this key are {\bf DirConstant},
{\bf ConstantWTHeight} and {\bf PLinearWTHeight}.  The choice
{\bf DirConstant} specifies that the saturation is constant on the whole
patch.  The choice {\bf ConstantWTHeight} specifies a constant height of
the water-table on the whole patch.  The choice {\bf PLinearWTHeight}
specifies that the height of the water-table on the patch will be given
by a piecewise linear function.

Note: the types {\bf ConstantWTHeight} and {\bf PLinearWTHeight} assume
we are running a 2-phase problem where phase 0 is the water phase.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCSaturation.water.Type  ConstantWTHeight
\end{verbatim}\end{display}

\pfkey{double}{Patch.{\em patch\_name}.BCSaturation.{\em phase\_name}.Value}{no default}
{
This key specifies either the constant saturation value if
{\bf DirConstant} is selected or the constant water-table height if
{\bf ConstantWTHeight} is selected.
}
\begin{display}\begin{verbatim}
pfset Patch.top.BCSaturation.air.Value 1.0
\end{verbatim}\end{display}


\pfkey{double}{Patch.{\em patch\_name}.BCSaturation.{\em phase\_name}.XLower}{no default}
{
This key specifies the lower $x$ coordinate of a line in the xy-plane if
type {\bf PLinearWTHeight} boundary conditions are specified.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCSaturation.water.XLower -10.0
\end{verbatim}\end{display}

\pfkey{double}{Patch.{\em patch\_name}.BCSaturation.{\em phase\_name}.YLower}{no default}
{
This key specifies the lower $y$ coordinate of a line in the xy-plane if
type {\bf PLinearWTHeight} boundary conditions are specified.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCSaturation.water.YLower 5.0
\end{verbatim}\end{display}

\pfkey{double}{Patch.{\em patch\_name}.BCSaturation.{\em phase\_name}.XUpper}{no default}
{
This key specifies the upper $x$ coordinate of a line in the xy-plane if
type {\bf PLinearWTHeight} boundary conditions are specified.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCSaturation.water.XUpper  125.0
\end{verbatim}\end{display}

\pfkey{double}{Patch.{\em patch\_name}.BCSaturation.{\em phase\_name}.YUpper}{no default}
{
This key specifies the upper $y$ coordinate of a line in the xy-plane if
type {\bf PLinearWTHeight} boundary conditions are specified.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCSaturation.water.YUpper  82.0
\end{verbatim}\end{display}

\pfkey{integer}{Patch.{\em patch\_name}.BCPressure.{\em phase\_name}.NumPoints}{no default}
{
This key specifies the number of points on which saturation data is
given along the line used for type {\bf DirEquilPLinear} boundary
conditions.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCPressure.water.NumPoints 2
\end{verbatim}\end{display}

\pfkey{double}{Patch.{\em patch\_name}.BCPressure.{\em phase\_name}.{\em point\_number}.Location}{no default}
{
This key specifies a number between 0 and 1 which represents the
location of a point on the line for which data is given in type
{\bf DirEquilPLinear} boundary conditions.  The line is parameterized
so that  0 corresponds to the lower end of the line, and 1 corresponds
to the upper end.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCPressure.water.0.Location 0.333
\end{verbatim}\end{display}

\pfkey{double}{Patch.{\em patch\_name}.BCPressure.{\em phase\_name}.{\em point\_number}.Value}{no default}
{
This key specifies the water-table height for the given point if type
{\bf DirEquilPLinear} boundary conditions are selected.  All saturation
values on the patch are determined by first projecting the water-table 
height value onto the line, then linearly interpolating between the
neighboring water-table height values onto the line.
}
\begin{display}\begin{verbatim}
pfset Patch.left.BCPressure.water.0.Value  4.5
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Initial Conditions: Phase Saturations}
\label{Initial Conditions: Phase Saturations}
\pfaddanchor{pfinInitialConditionsPhaseSaturations}

Note: this section needs to be defined {\em only} for multi-phase flow
and should {\em not} be defined for single phase and Richards' equation cases.

Here we define initial phase saturation conditions.
The format for this section of input is:

\pfkey{string}{ICSaturation.{\em phase\_name}.Type}{no default}
{
This key specifies the type of initial condition that will be applied
to different geometries for given phase, {\em phase\_name}.  The only
key currently available is {\bf Constant}.  The choice {\bf Constant}
will apply constants values within geometries for the phase. 
}
\begin{display}\begin{verbatim}
ICSaturation.water.Type Constant
\end{verbatim}\end{display}

\pfkey{string}{ICSaturation.{\em phase\_name}.GeomNames}{no default}
{
This key specifies the geometries on which an initial condition will be
given if the type is set to {\bf Constant}.

Note that geometries listed later ``overlay'' geometries listed earlier.
}
\begin{display}\begin{verbatim}
ICSaturation.water.GeomNames "domain"
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_input\_name}.ICSaturation.{\em phase\_name}.Value}{no default}
{
This key specifies the initial condition value assigned to all points in
the named geometry, {\em geom\_input\_name}, if the type was set to
{\bf Constant}.
}
\begin{display}\begin{verbatim}
Geom.domain.ICSaturation.water.Value 1.0
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Initial Conditions: Pressure}
\label{Initial Conditions: Pressure}

The keys in this section are used to specify pressure initial conditions for 
Richards' equation cases {\em only}.  These keys will be ignored if any other
case is run.

\pfkey{string}{ICPressure.Type}{no default}
{This key specifies the type of initial condition given.  The choices for this
key are {\bf Constant, HydroStaticDepth, HydroStaticPatch} and {\bf PFBFile}.
The choice {\bf Constant} specifies that the initial pressure will be constant
over the regions given.  The choice {\bf HydroStaticDepth} specifies that the
initial pressure within a region 
will be in hydrostatic equilibrium with a given pressure
specified at a given depth.  The choice {\bf HydroStaticPatch} specifies that
the initial pressure within a region 
will be in hydrostatic equilibrium with a given pressure on a specified patch.
Note that all regions must have the same type of initial data - different
regions cannot have different types of initial data.  However, the parameters
for the type may be different.
The {\bf PFBFile} specification means that the initial pressure 
will be taken as a spatially varying function given by data in a
\parflow{} binary (.pfb) file.  
}
\begin{display}\begin{verbatim}
pfset ICPressure.Type   Constant
\end{verbatim}\end{display}

\pfkey{list}{ICPressure.GeomNames}{no default}
{This key specifies the geometry names on which the initial pressure data will
be given.  These geometries must comprise the entire domain.
Note that conditions for regions that overlap other regions will have
unpredictable results.  The regions given must be disjoint.
}
\begin{display}\begin{verbatim}
pfset ICPressure.GeomNames   "toplayer middlelayer bottomlayer"
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.ICPressure.Value}{no default}
{This key specifies the initial pressure value for type {\bf Constant} initial
pressures and the reference pressure value for types {\bf HydroStaticDepth} and
{\bf HydroStaticPatch}.
}
\begin{display}\begin{verbatim}
pfset Geom.toplayer.ICPressure.Value  -734.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.ICPressure.RefElevation}{no default}
{This key specifies the reference elevation on which the reference pressure is
given for type {\bf HydroStaticDepth} initial pressures.
}
\begin{display}\begin{verbatim}
pfset Geom.toplayer.ICPressure.RefElevation  0.0
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.ICPressure.RefGeom}{no default}
{This key specifies the geometry on which the reference patch resides for type
{\bf HydroStaticPatch} initial pressures.
}
\begin{display}\begin{verbatim}
pfset Geom.toplayer.ICPressure.RefGeom   bottomlayer
\end{verbatim}\end{display}

\pfkey{double}{Geom.{\em geom\_name}.ICPressure.RefPatch}{no default}
{This key specifies the patch on which the reference pressure is given for type
{\bf HydorStaticPatch} initial pressures.
}
\begin{display}\begin{verbatim}
pfset Geom.toplayer.ICPressure.RefPatch   bottom
\end{verbatim}\end{display}

\pfkey{string}{Geom.{\em geom\_name}.ICPressure.FileName}{no default}
{This key specifies the name of the file containing pressure values for the
domain.  It is assumed that {\em geom\_name} is ``domain'' for this key.}
\begin{display}\begin{verbatim}
pfset Geom.domain.ICPressure.FileName  "ic_pressure.pfb"
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Initial Conditions: Phase Concentrations}
\label{Initial Conditions: Phase Concentrations}
\pfaddanchor{pfinInitialConditionsPhaseConcentrations}

Here we define initial concentration conditions for contaminants.
The format for this section of input is:

\pfkey{string}{PhaseConcen.{\em phase\_name}.{\em contaminant\_name}.Type}{no default}
{
This key specifies the type of initial condition that will be applied
to different geometries for given phase, {\em phase\_name}, and the
given contaminant, {\em contaminant\_name}.  The choices for this key
are {\bf Constant} or {\bf PFBFile}.  The choice {\bf Constant} will
apply constants values to different geometries.  The choice
{\bf PFBFile}  will read values from a ``ParFlow Binary'' file
(see \S~\ref{ParFlow Binary Files (.pfb)}).
}
\begin{display}\begin{verbatim}
PhaseConcen.water.tce.Type Constant
\end{verbatim}\end{display}

\pfkey{string}{PhaseConcen.{\em phase\_name}.GeomNames}{no default}
{
This key specifies the geometries on which an initial condition will be
given, if the type was set to {\bf Constant}.

Note that geometries listed later ``overlay'' geometries listed earlier.
}
\begin{display}\begin{verbatim}
PhaseConcen.water.GeomNames "ic_concen_region"
\end{verbatim}\end{display}

\pfkey{double}{PhaseConcen.{\em phase\_name}.{\em contaminant\_name}.{\em geom\_input\_name}.Value}{no default}
{
This key specifies the initial condition value assigned to all points in
the named geometry, {\em geom\_input\_name}, if the type was set to
{\bf Constant}.
}
\begin{display}\begin{verbatim}
PhaseConcen.water.tce.ic_concen_region.Value 0.001
\end{verbatim}\end{display}

\pfkey{string}{PhaseConcen.{\em phase\_name}.{\em contaminant\_name}.FileName}{no default}
{
This key specifies the name of the ``ParFlow Binary'' file which
contains the initial condition values if the type was set to {\bf PFBFile}.
}
\begin{display}\begin{verbatim}
PhaseConcen.water.tce.FileName "initial_concen_tce.pfb"
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Known Exact Solution}
\label{ExactSolution}

For {\em Richards equation cases only} we allow specification of an exact
solution to be used for testing the code.  
Only types that have been coded and predefined are allowed.
Note that if this is speccified as something other than no known solution,
corresponding boundary conditions and phase sources should also be specified.

\pfkey{string}{KnownSolution}{no default}
{This specifies the predefined function that will be used as the known
solution.  Possible choices for this key are {\bf NoKnownSolution, Constant, 
X, XPlusYPlusZ, X3Y2PlusSinXYPlus1, X3Y4PlusX2PlusSinXYCosYPlus1, XYZTPlus1}
and {\bf XYZTPlus1PermTensor}.
}
\begin{display}\begin{verbatim}
pfset KnownSolution  XPlusYPlusZ
\end{verbatim}\end{display}
Choices for this key correspond to solutions as follows.
\begin{description}
\item[{\bf NoKnownSolution}: ] No solution is known for this problem.
\item[{\bf Constant}: ] $p = {\rm constant}$
\item[{\bf X}: ] $p = x$
\item[{\bf XPlusYPlusZ}: ] $p = x + y + z$
\item[{\bf X3Y2PlusSinXYPlus1}: ] $p = x^3 y^2 + sin(xy) + 1$
\item[{\bf X3Y4PlusX2PlusSinXYCosYPlus1}: ] 
$p = x^3 y^4 + x^2 + \sin(xy)\cos y + 1$
\item[{\bf XYZTPlus1}: ] $p = xyzt + 1$
\item[{\bf XYZTPlus1}: ] $p = xyzt + 1$
\end{description}

\pfkey{double}{KnownSolution.Value}{no default}
{This key specifies the constant value of the known solution for type 
{\bf Constant} known solutions.
}
\begin{display}\begin{verbatim}
pfset KnownSolution.Value  1.0
\end{verbatim}\end{display}

Only for known solution test cases will information on the $L^2$-norm
of the pressure error be printed.

%=============================================================================
%=============================================================================

\subsection{Wells}
\label{Wells}
\pfaddanchor{pfinWells}

Here we define wells for the model.  The format for this section of
input is:

\pfkey{string}{Wells.Names}{no default}
{
This key specifies the names of the wells for which input data will be
given.
}
\begin{display}\begin{verbatim}
Wells.Names "test_well inj_well ext_well"
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.InputType}{no default}
{
This key specifies the type of well to be defined for the given well, 
{\em well\_name}.  This key can be either {\bf Vertical} or
{\bf Recirc}.  The value {\bf Vertical} indicates that this is a
single segmented well whose action will be specified by the user.
The value {\bf Recirc} indicates that this is a dual segmented,
recirculating, well with one segment being an extraction well and
another being an injection well.  The extraction well filters out a
specified fraction of each contaminant and recirculates the remainder
to the injection well where the diluted fluid is injected back in.  The
phase saturations at the extraction well are passed without modification
to the injection well.

Note with the recirculating well, several input options are not needed
as the extraction well will provide these values to the injection well.
}
\begin{display}\begin{verbatim}
Wells.test_well.InputType Vertical
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.Action}{no default}
{
This key specifies the pumping action of the well.  This key can be
either {\bf Injection} or {\bf Extraction}.  A value of {\bf Injection}
indicates that this is an injection well.  A value of {\bf Extraction}
indicates that this is an extraction well.
}
\begin{display}\begin{verbatim}
Wells.test_well.Action Injection
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.Type}{no default}
{
This key specfies the mechanism by which the well works (how \parflow{}
works with the well data) if the input type key is set to {\bf Vectical}.
This key can be either {\bf Pressure} or {\bf Flux}.  A value of
{\bf Pressure} indicates that the data provided for the well is in terms
of hydrostatic pressure and \parflow{} will ensure that the computed
pressure field satisfies this condition in the computational cells which
define the well.  A value of {\bf Flux} indicates that the data provided
is in terms of volumetric flux rates and \parflow{} will ensure that the
flux field satisfies this condition in the computational cells which
define the well.
}
\begin{display}\begin{verbatim}
Wells.test_well.Type Flux
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.ExtractionType}{no default}
{
This key specfies the mechanism by which the extraction well works (how
\parflow{} works with the well data) if the input type key is set to
{\bf Recirc}.  This key can be either {\bf Pressure} or {\bf Flux}.  A
value of {\bf Pressure} indicates that the data provided for the well is
in terms of hydrostatic pressure and \parflow{} will ensure that the
computed pressure field satisfies this condition in the computational
cells which define the well.  A value of {\bf Flux} indicates that the
data provided is in terms of volumetric flux rates and \parflow{} will
ensure that the flux field satisfies this condition in the computational
cells which define the well.
}
\begin{display}\begin{verbatim}
Wells.ext_well.ExtractionType Pressure
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.InjectionType}{no default}
{
This key specfies the mechanism by which the injection well works (how
\parflow{} works with the well data) if the input type key is set to
{\bf Recirc}.  This key can be either {\bf Pressure} or {\bf Flux}.  A
value of {\bf Pressure} indicates that the data provided for the well is
in terms of hydrostatic pressure and \parflow{} will ensure that the
computed pressure field satisfies this condition in the computational
cells which define the well.  A value of {\bf Flux} indicates that the
data provided is in terms of volumetric flux rates and \parflow{} will
ensure that the flux field satisfies this condition in the computational
cells which define the well.
}
\begin{display}\begin{verbatim}
Wells.inj_well.InjectionType Flux
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.X}{no default}
{
This key specifies the x location of the vectical well if the input
type is set to {\bf Vectical} or of both the extraction and injection
wells if the input type is set to {\bf Recirc}.
}
\begin{display}\begin{verbatim}
Wells.test_well.X 20.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.Y}{no default}
{
This key specifies the y location of the vectical well if the input
type is set to {\bf Vectical} or of both the extraction and injection
wells if the input type is set  to {\bf Recirc}.
}
\begin{display}\begin{verbatim}
Wells.test_well.Y 36.5
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.ZUpper}{no default}
{
This key specifies the z location of the upper extent of a vectical well
if the input type is set to {\bf Vectical}.
}
\begin{display}\begin{verbatim}
Wells.test_well.ZUpper 8.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.ExtractionZUpper}{no default}
{
This key specifies the z location of the upper extent of a extraction
well if the input type is set to {\bf Recirc}.
}
\begin{display}\begin{verbatim}
Wells.ext_well.ExtractionZUpper 3.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.InjectionZUpper}{no default}
{
This key specifies the z location of the upper extent of a injection
well if the input type is set to {\bf Recirc}.
}
\begin{display}\begin{verbatim}
Wells.inj_well.InjectionZUpper 6.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.ZLower}{no default}
{
This key specifies the z location of the lower extent of a vectical well
if the input type is set to {\bf Vectical}.
}
\begin{display}\begin{verbatim}
Wells.test_well.ZLower 2.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.ExtractionZLower}{no default}
{
This key specifies the z location of the lower extent of a extraction
well if the input type is set to {\bf Recirc}.
}
\begin{display}\begin{verbatim}
Wells.ext_well.ExtractionZLower 1.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.InjectionZLower}{no default}
{
This key specifies the z location of the lower extent of a injection
well if the input type is set to {\bf Recirc}.
}
\begin{display}\begin{verbatim}
Wells.inj_well.InjectionZLower 4.0
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.Method}{no default}
{
This key specifies a method by which pressure or flux for a vertical
well will be weighted before assignment to computational cells.  This
key can only be {\bf Standard} if the type key is set to {\bf Pressure};
or this key can be either {\bf Standard}, {\bf Weighted} or
{\bf Patterned} if the type key is set to {\bf Flux}.  A value of
{\bf Standard} indicates that the pressure or flux data will be used as
is.  A value of {\bf Weighted} indicates that the flux data is to be
weighted by the cells permeability divided by the sum of all cell
permeabilities which define the well.  The value of {\bf Patterned} is
not implemented.
}
\begin{display}\begin{verbatim}
Wells.test_well.Method Weighted
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.ExtractionMethod}{no default}
{
This key specifies a method by which pressure or flux for an extraction
well will be weighted before assignment to computational cells.  This
key can only be {\bf Standard} if the type key is set to {\bf Pressure};
or this key can be either {\bf Standard}, {\bf Weighted} or
{\bf Patterned} if the type key is set to {\bf Flux}.  A value of
{\bf Standard} indicates that the pressure or flux data will be used as
is.  A value of {\bf Weighted} indicates that the flux data is to be
weighted by the cells permeability divided by the sum of all cell
permeabilities which define the well.  The value of {\bf Patterned} is
not implemented.
}
\begin{display}\begin{verbatim}
Wells.ext_well.ExtractionMethod Standard
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.InjectionMethod}{no default}
{
This key specifies a method by which pressure or flux for an injection
well will be weighted before assignment to computational cells.  This
key can only be {\bf Standard} if the type key is set to {\bf Pressure};
or this key can be either {\bf Standard}, {\bf Weighted} or
{\bf Patterned} if the type key is set to {\bf Flux}.  A value of
{\bf Standard} indicates that the pressure or flux data will be used as
is.  A value of {\bf Weighted} indicates that the flux data is to be
weighted by the cells permeability divided by the sum of all cell
permeabilities which define the well.  The value of {\bf Patterned} is
not implemented.
}
\begin{display}\begin{verbatim}
Wells.inj_well.InjectionMethod Standard
\end{verbatim}\end{display}

\pfkey{string}{Wells.{\em well\_name}.Cycle}{no default}
{
This key specifies the time cycles to which data for the well
{\em well\_name} corresponds. 
}
\begin{display}\begin{verbatim}
Wells.test_well.Cycle "all_time"
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Pressure.Value}{no default}
{
This key specifies the hydrostatic pressure value for a vectical well
if the type key is set to {\bf Pressure}.

Note This value gives the pressure of the primary phase (water) at
$z=0$.  The other phase pressures (if any) are computed from the physical
relationships that exist between the phases.
}
\begin{display}\begin{verbatim}
Wells.test_well.all_time.Pressure.Value 6.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Extraction.Pressure.Value}{no default}
{
This key specifies the hydrostatic pressure value for an extraction well
if the extraction type key is set to {\bf Pressure}.

Note This value gives the pressure of the primary phase (water) at
$z=0$.  The other phase pressures (if any) are computed from the physical
relationships that exist between the phases.
}
\begin{display}\begin{verbatim}
Wells.ext_well.all_time.Extraction.Pressure.Value 4.5
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Injection.Pressure.Value}{no default}
{
This key specifies the hydrostatic pressure value for an injection well
if the injection type key is set to {\bf Pressure}.

Note This value gives the pressure of the primary phase (water) at
$z=0$.  The other phase pressures (if any) are computed from the physical
relationships that exist between the phases.
}
\begin{display}\begin{verbatim}
Wells.inj_well.all_time.Injection.Pressure.Value 10.2
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Flux.{\em phase\_name}.Value}{no default}
{
This key specifies the volumetric flux for a vectical well if the type
key is set to {\bf Flux}.

Note only a positive number should be entered, \parflow{} assignes the
correct sign based on the chosen action for the well.
}
\begin{display}\begin{verbatim}
Wells.test_well.all_time.Flux.water.Value 250.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Extraction.Flux.{\em phase\_name}.Value}{no default}
{
This key specifies the volumetric flux for an extraction well if the
extraction type key is set to {\bf Flux}.

Note only a positive number should be entered, \parflow{} assignes the
correct sign based on the chosen action for the well.
}
\begin{display}\begin{verbatim}
Wells.ext_well.all_time.Extraction.Flux.water.Value 125.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Injection.Flux.{\em phase\_name}.Value}{no default}
{
This key specifies the volumetric flux for an injection well if the
injection type key is set to {\bf Flux}.

Note only a positive number should be entered, \parflow{} assignes the
correct sign based on the chosen action for the well.
}
\begin{display}\begin{verbatim}
Wells.inj_well.all_time.Injection.Flux.water.Value 80.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Saturation.{\em phase\_name}.Value}{no default}
{
This key specifies the saturation value of a vertical well.
}
\begin{display}\begin{verbatim}
Wells.test_well.all_time.Saturation.water.Value 1.0
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Concentration.{\em phase\_name}.{\em contaminant\_name}.Value}{no default}
{
This key specifies the contaminant value of a vertical well.
}
\begin{display}\begin{verbatim}
Wells.test_well.all_time.Concentration.water.tce.Value 0.0005
\end{verbatim}\end{display}

\pfkey{double}{Wells.{\em well\_name}.{\em interval\_name}.Injection.Concentration.{\em phase\_name}.{\em contaminant\_name}.Fraction}{no default}
{
This key specifies the fraction of the extracted contaminant which gets
resupplied to the injection well.
}
\begin{display}\begin{verbatim}
Wells.inj_well.all_time.Injection.Concentration.water.tce.Fraction 0.01
\end{verbatim}\end{display}

\noindent
Multiple wells assigned to one grid location can occur in several
instances.  The current actions taken by the code are as follows:

\begin{itemize}
\item If multiple pressure wells are assigned to one grid cell, the code
      retains only the last set of overlapping well values entered.
\item If multiple flux wells are assigned to one grid cell, the code sums the
      contributions of all overlapping wells to get one effective well flux.
\item If multiple pressure and flux wells are assigned to one grid cell, the
      code retains the last set of overlapping hydrostatic pressure values
      entered and sums all the overlapping flux well values to get an
      effective pressure/flux well value.
\end{itemize}

%=============================================================================
%=============================================================================

\subsection{Code Parameters}
\label{Code Parameters}
\pfaddanchor{pfinCodeParameters}

In addition to input keys related to the physics capabilities
and modeling specifics there are some key values used by various
algorithms and general control flags for \parflow{}.  These
are described next :

\pfkey{string}{Solver.Linear}{PCG}
{This key specifies the linear solver used for solver {\bf IMPES}.  Choices for this key are {\bf MGSemi, PPCG, PCG} and {\bf   CGHS}.  The choice
%% == we need to check the solvers used here to match the descriptions
{\bf MGSemi} is an algebraic mulitgrid linear solver (not a preconditioned conjugate gradient) which may be less robust than {\bf PCG} as described in \cite{Ashby-Falgout90}.  The choice {\bf PPCG} is a preconditioned conjugate gradient solver.
The choice {\bf PCG} is a conjugate gradient solver with a multigrid preconditioner.  The choice {\bf CGHS} is a conjugate gradient solver.
}
\begin{display}\begin{verbatim}
pfset Solver.Linear   MGSemi
\end{verbatim}\end{display}

\pfkey{integer}{Solver.SadvectOrder}{2}
{
This key controls the order of the explicit method used in
advancing the saturations.  This value can be either 1 for
a standard upwind first order or 2 for a second order
Godunov method.
}
\begin{display}\begin{verbatim}
pfset Solver.SadvectOrder 1
\end{verbatim}\end{display}

\pfkey{integer}{Solver.AdvectOrder}{2}
{
This key controls the order of the explicit method used in
advancing the concentrations.  This value can be either 1 for
a standard upwind first order or 2 for a second order
Godunov method.
}
\begin{display}\begin{verbatim}
pfset Solver.AdvectOrder 2
\end{verbatim}\end{display}

\pfkey{double}{Solver.CFL}{0.7}
{
This key gives the value of the weight put on the computed
CFL limit before computing a global timestep value.  Values
greater than 1 are not suggested and in fact because this is
an approximation, values slightly less than 1 can also produce
instabilities.
}
\begin{display}\begin{verbatim}
pfset Solver.CFL 0.7
\end{verbatim}\end{display}

\pfkey{integer}{Solver.MaxIter}{1000000}
{
This key gives the maximum number of iterations that will
be allowed for time-stepping.  This is to prevent a run-away
simulation.
}
\begin{display}\begin{verbatim}
pfset Solver.MaxIter 100
\end{verbatim}\end{display}



\pfkey{double}{Solver.RelTol}{1.0}
{
This value gives the relative tolerance for the linear
solve algorithm.
}
\begin{display}\begin{verbatim}
pfset Solver.RelTol 1.0
\end{verbatim}\end{display}

\pfkey{double}{Solver.AbsTol}{1E-9}
{
This value gives the absolute tolerance for the linear
solve algorithm.
}
\begin{display}\begin{verbatim}
pfset Solver.AbsTol 1E-8
\end{verbatim}\end{display}

\pfkey{double}{Solver.Drop}{1E-8}
{
This key gives a clipping value for data written to PFSB
files.  Data values greater than the negative of this
value and less than the value itself are treated as zero
and not written to PFSB files.
}
\begin{display}\begin{verbatim}
pfset Solver.Drop 1E-6
\end{verbatim}\end{display}

\pfkey{string}{Solver.PrintSubsurf}{True}
{
This key is used to turn on printing of the subsurface data,
Permeability and Porosity.  The data is printed after it is
generated and before the main time stepping loop - only once
during the run.  The data is written as a PFB file.
}
\begin{display}\begin{verbatim}
pfset Solver.PrintSubsurf False
\end{verbatim}\end{display}

\pfkey{string}{Solver.PrintPressure}{True}
{
This key is used to turn on printing of the pressure data.
The printing of the data is controlled by values in the
timing information section.  The data is written as a PFB
file.
}
\begin{display}\begin{verbatim}
pfset Solver.PrintPressure False
\end{verbatim}\end{display}

\pfkey{string}{Solver.PrintVelocities}{False}
{
This key is used to turn on printing of the x, y and z
velocity data.  The printing of the data is controlled by
values in the timing information section.  The data is
written as a PFB file.
}
\begin{display}\begin{verbatim}
pfset Solver.PrintVelocities True
\end{verbatim}\end{display}

\pfkey{string}{Solver.PrintSaturation}{True}
{
This key is used to turn on printing of the saturation data.
The printing of the data is controlled by values in the
timing information section.  The data is written as a PFB file.
}
\begin{display}\begin{verbatim}
pfset Solver.PrintSaturation False
\end{verbatim}\end{display}

\pfkey{string}{Solver.PrintConcentration}{True}
{
This key is used to turn on printing of the concentration data.
The printing of the data is controlled by values in the
timing information section.  The data is written as a PFSB file.
}
\begin{display}\begin{verbatim}
pfset Solver.PrintConcentration False
\end{verbatim}\end{display}

\pfkey{string}{Solver.PrintWells}{True}
{
This key is used to turn on collection and printing of the
well data.  The data is collected at intervals given by values
in the timing information section.  Printing occurs at the
end of the run when all collected data is written. 
}
\begin{display}\begin{verbatim}
pfset Solver.PrintWells False
\end{verbatim}\end{display}

\pfkey{string}{Solver.PrintLSMSink}{False}
{
This key is used to turn on printing of the
flux array passed from \code{CLM} to \parflow{}.   Printing occurs at each {\bf DumpInterval} time. 
}
\begin{display}\begin{verbatim}
pfset Solver.PrintLSMSink True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloSubsurfData}{False}
{
This key is used to specify printing of the subsurface data,
Permeability and Porosity in silo binary file format.  The data is printed after it is
generated and before the main time stepping loop - only once
during the run.  This data may be read in by VisIT and other visualization packages.
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloSubsurfData True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloPressure}{False}
{
This key is used to specify printing of the saturation data in silo binary format.
The printing of the data is controlled by values in the
timing information section. This data may be read in by VisIT and other visualization packages.
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloPressure True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloSaturation}{False}
{
This key is used to specify printing of the saturation data using silo binary format.
The printing of the data is controlled by values in the
timing information section.  
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloSaturation True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloConcentration}{False}
{
This key is used to specify printing of the concentration data in silo binary format.
The printing of the data is controlled by values in the
timing information section.
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloConcentration True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloVelocities}{False}
{
This key is used to specify printing of the x, y and z
velocity data in silo binary format.  The printing of the data is controlled by
values in the timing information section.
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloVelocities True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloSlopes}{False}
{
This key is used to specify printing of the x and y slope data using silo binary format.
The printing of the data is controlled by values in the
timing information section.  
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloSlopes  True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloMannings}{False}
{
This key is used to specify printing of the Manning's roughness data in silo binary format.
The printing of the data is controlled by values in the
timing information section.
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloMannings True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloSpecificStorage}{False}
{
This key is used to specify printing of the specific storage data in silo binary format.
The printing of the data is controlled by values in the
timing information section.
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloSpecificStorage True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloMask}{False}
{
This key is used to specify printing of the mask data using silo binary format.  The mask contains values equal to one for active cells and zero for inactive cells. 
The printing of the data is controlled by values in the
timing information section.  
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloMask  True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloEvapTrans}{False}
{
This key is used to specify printing of the evaporation and rainfall flux data using silo binary format.  This data comes from either \code{clm} or from external calls to \parflow{} such as \emph{WRF}.  This data is in units of $[L^3 T^{-1}]$.
The printing of the data is controlled by values in the
timing information section.  
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloEvapTrans  True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloEvapTransSum}{False}
{
This key is used to specify printing of the evaporation and rainfall flux data using silo binary format as a running, cumulative amount.  This data comes from either \code{clm} or from external calls to \parflow{} such as \emph{WRF}.  This data is in units of $[L^3]$.
The printing of the data is controlled by values in the
timing information section.  
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloEvapTransSum  True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloOverlandSum}{False}
{
This key is used to specify calculation and printing of the total overland outflow from the domain using silo binary format as a running cumulative amount.   This is integrated along all domain boundaries and is calculated any location that slopes at the edge of the domain point outward. This data is in units of $[L^3]$.
The printing of the data is controlled by values in the
timing information section.  
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloOverlandSum  True
\end{verbatim}\end{display}

\pfkey{string}{Solver.TerrainFollowingGrid}{False}
{
This key specifies that a terrain-following coordinate transform is used for solver Richards. This key sets x and y subsurface slopes to be the same as the Topographic slopes (a value of False sets these subsurface slopes to zero).
These slopes are used in the Darcy fluxes to add a density, gravity -dependent term.  This key will not change the output files (that is the output is still orthogonal) or the geometries (they will still follow the computational grid)-- these two things are both to do items.  This key only changes solver Richards, not solver Impes. 
}

\begin{display}\begin{verbatim}
pfset Solver.TerrainFollowingGrid                        True
\end{verbatim}\end{display}


%=============================================================================
%=============================================================================

\subsection{SILO Options}
\label{SILO Options}

The following keys are used to control how SILO writes data.  SILO
 allows writing to PDB and HDF5 file formats.  SILO also allows data
 compression to be used, which can save signicant amounts of disk
 space for some problems.  

\pfkey{string}{SILO.Filetype}{PDB}
{
This key is used to specify the SILO filetype.   Allowed values are PDB and HDF5.  Note that you must have configured SILO with HDF5 in order to use that option.}
\begin{display}\begin{verbatim}
pfset SILO.Filetype  PDB
\end{verbatim}\end{display}


\pfkey{string}{SILO.CompressionOptions}{}
{
This key is used to specify the SILO compression options. See the SILO manual for the
 DB\_SetCompression command for information on available options.   NOTE: the options avaialable are highly dependent on the configure options when building SILO.}
\begin{display}\begin{verbatim}
pfset SILO.CompressionOptions  ``METHOD=GZIP''
\end{verbatim}\end{display}


%=============================================================================
%=============================================================================

\subsection{Richards' Equation Solver Parameters}
\label{RE Solver Parameters}

The following keys are used to specify various parameters used by the linear
and nonlinear solvers in the Richards' equation implementation.
For information about these solvers, see 
\cite{Woodward98} and \cite{Ashby-Falgout90}.

\pfkey{double}{Solver.Nonlinear.ResidualTol}{1e-7}
{This key specifies the tolerance that measures how much the relative 
reduction in the nonlinear residual should be before nonlinear 
iterations stop.  The magnitude of the residual is measured with 
the $l^1$ (max) norm.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.ResidualTol   1e-4
\end{verbatim}\end{display}

\pfkey{double}{Solver.Nonlinear.StepTol}{1e-7}
{This key specifies the tolerance that measures how small the 
difference between two consecutive nonlinear steps can be before 
nonlinear iterations stop.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.StepTol   1e-4
\end{verbatim}\end{display}

\pfkey{integer}{Solver.Nonlinear.MaxIter}{15}
{This key specifies the maximum number of nonlinear iterations allowed before
iterations stop with a convergence failure.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.MaxIter   50
\end{verbatim}\end{display}


\pfkey{integer}{Solver.Linear.KrylovDimension}{10}
{This key specifies the maximum number of vectors to be used in setting up the
Krylov subspace in the GMRES iterative solver.  These vectors are of problem
size and it should be noted that large increases in this parameter can limit
problem sizes.  However, increasing this parameter can sometimes help nonlinear
solver convergence.
}
\begin{display}\begin{verbatim}
pfset Solver.Linear.KrylovDimension   15
\end{verbatim}\end{display}

\pfkey{integer}{Solver.Linear.MaxRestarts}{0}
{This key specifies the number of restarts allowed to the GMRES solver.
Restarts start the development of the Krylov subspace over using the current
iterate as the initial iterate for the next pass.
}
\begin{display}\begin{verbatim}
pfset Solver.Linear.MaxRestarts   2
\end{verbatim}\end{display}

\pfkey{integer}{Solver.MaxConvergencFailures}{3}
{This key gives the maximum number of convergence failures
allowed.   Each convergence failure cuts the timestep 
in half and the solver tries to advance the solution with the reduced
timestep.

The default value is 3.   

Note that setting this value to a value greater than 9 may result in
errors in how time cycles are calculated.  Time is discretized in
terms of the base time unit and if the solver begins to take very
small timesteps \(smaller than base time unit \/ 1000\) the values based
on time cycles will be change at slightly incorrect times. If the
problem is failing converge so poorly that a large number of restarts are
required, consider setting the timestep to a smaller value.
}
\begin{display}\begin{verbatim}
pfset Solver.MaxConvergenceFailures 4
\end{verbatim}\end{display}

\pfkey{string}{Solver.Nonlinear.PrintFlag}{HighVerbosity}
{This key specifies the amount of informational data that is printed to the
\code{*.out.kinsol.log} file.  Choices for this key are {\bf NoVerbosity, 
LowVerbosity, NormalVerbosity} and {\bf  HighVerbosity}.  The choice 
{\bf NoVerbosity} prints no statistics about the nonlinear convergence
process.  The choice {\bf LowVerbosity} outputs the nonlinear iteration count, 
the scaled norm of the nonlinear function, and the number of function calls.
The choice {\bf NormalVerbosity} prints the same as for {\bf LowVerbosity} 
and also the global strategy statistics.  The choice {\bf HighVerbosity} prints
the same as for {\bf NormalVerbosity} with the addition of further Krylov 
iteration statistics.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.PrintFlag   NormalVerbosity
\end{verbatim}\end{display}

\pfkey{string}{Solver.Nonlinear.EtaChoice}{Walker2}
{This key specifies how the linear system tolerance will be selected.  
The linear system is solved until a relative residual reduction of $\eta$ 
is achieved.  Linear residuall norms are measured in the $l^2$ norm.
Choices for this key include {\bf EtaConstant, Walker1} and 
{\bf Walker2}.  If the choice {\bf EtaConstant} is specified, then $\eta$ will
be taken as constant.  The choices {\bf Walker1} and {\bf Walker2} specify
choices for $\eta$ developed by Eisenstat and Walker \cite{EW96}.  The choice
{\bf Walker1} specifies that $\eta$ will be given by
$| \|F(u^k)\| - \|F(u^{k-1}) + J(u^{k-1})*p \|  |  / \|F(u^{k-1})\|$. The
choice {\bf Walker2} specifies that $\eta$ will be given by
$\gamma \|F(u^k)\| / \|F(u^{k-1})\|^{\alpha}$.  For both of the last two
choices, $\eta$ is never allowed to be less than 1e-4.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.EtaChoice   EtaConstant
\end{verbatim}\end{display}


\pfkey{double}{Solver.Nonlinear.EtaValue}{1e-4}
{This key specifies the constant value of $\eta$ for the EtaChoice key
{\bf EtaConstant}.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.EtaValue   1e-7
\end{verbatim}\end{display}

\pfkey{double}{Solver.Nonlinear.EtaAlpha}{2.0}
{This key specifies the value of $\alpha$ for the case of EtaChoice being 
{\bf Walker2}. 
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.EtaAlpha   1.0
\end{verbatim}\end{display}

\pfkey{double}{Solver.Nonlinear.EtaGamma}{0.9}
{This key specifies the value of $\gamma$ for the case of EtaChoice being
{\bf Walker2}.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.EtaGamma   0.7
\end{verbatim}\end{display}

\pfkey{string}{Solver.Nonlinear.UseJacobian}{False}
{This key specifies whether the Jacobian will be used in matrix-vector products
or whether a matrix-free version of the code will run.  Choices for this key 
are {\bf False} and {\bf True}.
Using the Jacobian will most likely decrease the number of nonlinear iterations
but require more memory to run.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.UseJacobian   True
\end{verbatim}\end{display}

\pfkey{double}{Solver.Nonlinear.DerivativeEpsilon}{1e-7}
{This key specifies the value of $\epsilon$ used in approximating the action of
the Jacobian on a vector with approximate directional derivatives of the
nonlinear function.  This parameter is only used when the UseJacobian key is
{\bf False}.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.DerivativeEpsilon   1e-8
\end{verbatim}\end{display}

\pfkey{string}{Solver.Nonlinear.Globalization}{LineSearch}
{This key specifies the type of global strategy to use.  Possible choices for
this key are {\bf InexactNewton} and {\bf LineSearch}.  The choice {\bf
InexactNewton} specifies no global strategy, and the choice {\bf LineSearch}
specifies that a line search strategy should be used where the nonlinear step
can be lengthened or decreased to satisfy certain criteria.
}
\begin{display}\begin{verbatim}
pfset Solver.Nonlinear.Globalization   LineSearch
\end{verbatim}\end{display}

\pfkey{string}{Solver.Linear.Preconditioner}{MGSemi}
{This key specifies which preconditioner to use.  Currently, the three choices
are {\bf NoPC, MGSemi, PFMG, PFMGOctree} and {\bf SMG}.  The choice {\bf NoPC} specifies that no
preconditioner should be used.  The choice {\bf MGSemi} specifies 
a semi-coarsening multigrid algorithm which uses a point relaxation method. 
The choice {\bf SMG} specifies a semi-coarsening multigrid algorithm which uses
plane relaxations.  This method is more robust than {\bf MGSemi}, but generally
requires more memory and compute time. The choice {\bf PFMGOctree} can be more efficient for problems with large numbers of inactive cells.
}
\begin{display}\begin{verbatim}
pfset Solver.Linear.Preconditioner   MGSemi
\end{verbatim}\end{display}


\pfkey{string}{Solver.Linear.Preconditioner.SymmetricMat}{Symmetric}
{This key specifies whether the preconditioning matrix is symmetric.
Choices fo rthis key are {\bf Symmetric} and {\bf Nonsymmetric}.
The choice {\bf Symmetric} specifies that the symmetric part of the Jacobian
will be used as the preconditioning matrix.  The choice {\bf Nonsymmetric}
specifies that the full Jacobian will be used as the preconditioning matrix.
NOTE: ONLY {\bf Symmetric} CAN BE USED IF MGSemi IS THE SPECIFIED
PRECONDITIONER! 
}
\begin{display}\begin{verbatim}
pfset Solver.Linear.Preconditioner.SymmetricMat     Symmetric
\end{verbatim}\end{display}

\pfkey{integer}{Solver.Linear.Preconditioner.{\em precond\_method}.MaxIter}{1}
{This key specifies the maximum number of iterations to take in solving the 
preconditioner system with {\em precond\_method} solver.
}
\begin{display}\begin{verbatim}
pfset Solver.Linear.Preconditioner.SMG.MaxIter    2
\end{verbatim}\end{display}

\pfkey{integer}{Solver.Linear.Preconditioner.SMG.NumPreRelax}
{1}
{This key specifies the number of relaxations to take before coarsening in the
specified preconditioner method.  Note that this key is only relevant to
the SMG multigrid preconditioner.
}
\begin{display}\begin{verbatim}
pfset Solver.Linear.Preconditioner.SMG.NumPreRelax    2
\end{verbatim}\end{display}

\pfkey{integer}
{Solver.Linear.Preconditioner.SMG.NumPostRelax}{1}
{This key specifies the number of relaxations to take after coarsening in the
specified preconditioner method.  Note that this key is only relevant to
the SMG multigrid preconditioner.
}
\begin{display}\begin{verbatim}
pfset Solver.Linear.Preconditioner.SMG.NumPostRelax    0
\end{verbatim}\end{display}

\pfkey{string}{Solver.Linear.Preconditioner.PFMG.RAPType}{NonGalerkin} {For the
  PFMG solver, this key specifies the {\em Hypre} RAP type.  Valid values
  are {\bf Galerkin} or {\bf NonGalerkin} }
\begin{display}\begin{verbatim}
pfset Solver.Linear.Preconditioner.PFMG.RAPType    Galerkin
\end{verbatim}\end{display}


\pfkey{logical}{Solver.EvapTransFile}{False}
{This key specifies specifies that the Flux terms for Richards' equation are read in from a \file{.pfb} file.  This file has $[T^-1]$
units. Note this key is for a steady-state flux and should \emph{not} be used in conjunction with the transient key below.}
\begin{display}\begin{verbatim}
pfset Solver.EvapTransFile    True
\end{verbatim}\end{display}

\pfkey{logical}{Solver.EvapTransFileTransient}{False}
{This key specifies specifies that the Flux terms for Richards' equation are read in from a series of flux \file{.pfb} file.  Each file has $[T^-1]$
	units. Note this key should not be used with the key above, only one of these keys should be set to \code{True} at a time, not both.
}
\begin{display}\begin{verbatim}
pfset Solver.EvapTransFileTransient    True
\end{verbatim}\end{display}


\pfkey{string}{Solver.EvapTrans.FileName}{no default}
{This key specifies specifies filename for the distributed \file{.pfb} file that contains the flux values for Richards' equation.  This file has $[T^-1]$
units.  For the steady-state option (\emph{Solver.EvapTransFile=\bf{True}}) this key should be the complete filename.  For the transient option 
(\emph{Solver.EvapTransFileTransient=\bf{True}} then the filename is a header and \parflow{} will load one file per timestep, with the form \file{filename.00000.pfb}.}
\begin{display}\begin{verbatim}
pfset Solver.EvapTrans.FileName   evap.trans.test.pfb
\end{verbatim}\end{display}

\pfkey{string}{Solver.LSM}{none}
{This key specifies whether a land surface model, such as \code{CLM}, will be called each solver timestep.
Choices for this key include {\bf none} and {\bf CLM}. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.LSM CLM
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\subsection{Spinup Options}
\label{Spinup Options}

These keys allow for \emph{reduced or dampened physics} during model spinup or initialization. They are {\bf only} intended
for these initialization periods, {\bf not} for regular runtime. 

\pfkey{integer}{OverlandFlowSpinUp}{0}
{This key specifies that a \emph{simplified} form of the overland flow boundary condition (Equation \ref{eq:overland_bc}) be used in place of the full equation.
This formulation \emph{removes lateral flow} and drives and ponded water pressures to zero.  While this can be helpful in spinning up the subsurface,
this is no longer coupled subsurface-surface flow.  If set to zero (the default) this key behaves normally. 
}
\begin{display}\begin{verbatim}
pfset OverlandFlowSpinUp   1
\end{verbatim}\end{display}

\pfkey{double}{OverlandFlowSpinUpDampP1}{0.0}
{This key sets $P_1$ and provides exponential dampening to the pressure relationship in the overland flow equation by adding the 
following term: $P_2*exp(\psi*P_2)$}
\begin{display}\begin{verbatim}
pfset OverlandSpinupDampP1  10.0
\end{verbatim}\end{display}

\pfkey{double}{OverlandFlowSpinUpDampP2}{0.0}
{This key sets $P_2$ and provides exponential dampening to the pressure relationship in the overland flow equation adding the 
following term: $P_2*exp(\psi*P_2)$ }
\begin{display}\begin{verbatim}
pfset OverlandSpinupDampP2  0.1
\end{verbatim}\end{display}

%======
% new keys for CLM distributed forcing, restarts, etc here  @RMM
%=======
\subsection{CLM Solver Parameters}
\label{CLM Solver Parameters}

\pfkey{string}{Solver.CLM.Print1dOut}{False}
{This key specifies whether the \code{CLM} one dimensional (averaged over each processor) output file is written or not. 
Choices for this key include {\bf True} and {\bf False}. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.Print1dOut   False
\end{verbatim}\end{display}


\pfkey{integer}{Solver.CLM.IstepStart}{1}
{This key specifies the value of the counter, {\em istep} in \code{CLM}.  This key primarily determines the start of the output counter for \code{CLM}.It is used to restart a run by setting the key to the ending step of the previous run plus one. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.IstepStart     8761      
\end{verbatim}\end{display}


\pfkey{String}{Solver.CLM.MetForcing}{no default}
{This key specifies defines whether 1D (uniform over  
the domain), 2D (spatially distributed) or 3D (spatially distributed with multiple timesteps per \code{.pfb} forcing file) forcing data is used.  Choices for this key are {\bf 1D}, {\bf 2D} and {\bf 3D}. This key  
has no default so the user {\em must} set it to 1D, 2D or 3D. Failure to set this key will cause \code{CLM} to still be run but with unpredictable values causing \code{CLM} to eventually crash. 1D meteorological forcing files are text files with single columns for each variable and each timestep per row, while 2D forcing files are distributed \parflow{} binary files, one for each variable and timestep.  File names are specified in the {\bf Solver.CLM.MetFileName} variable below. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.MetForcing   2D       
\end{verbatim}\end{display}

\pfkey{String}{Solver.CLM.MetFileName}{no default}
{This key specifies defines the file name for 1D, 2D or 3D forcing data.  1D meteorological forcing files are text files with single columns for each variable and each timestep per row, while 2D and 3D forcing files are distributed \parflow{} binary files, one for each variable and timestep (2D) or one for each variable and \emph{multiple} timesteps (3D). Behavior of this key is different for 1D and 2D and 3D cases, as sepcified by the {\bf Solver.CLM.MetForcing} key above. For 1D cases, it is the {\em FULL FILE NAME}. Note that in this configuration, this forcing file is {\bf not} distributed,  the user does not provide copies such as \code{narr.1hr.txt.0}, \code{narr.1hr.txt.1} for each processor.  \parflow{} only needs the single original file ({\em e.g.} narr.1hr.txt).  For 2D cases, this key is the {\em BASE FILE NAME} for the 2D forcing files, currently set to NLDAS, with individual files determined as follows \code{NLDAS.<variable>.<time step>.pfb}. Where the \code{<variable>} is the forcing variable and \code{<timestep>} is the integer file counter corresponding to {\em istep} above. Forcing is needed for following variables:
\begin{description}
\item[{\bf DSWR}: ] Downward Visible or Short-Wave radiation $[W/m^2]$.
\item[{\bf DLWR}: ] Downward Infa-Red or Long-Wave radiation $[W/m^2]$
\item[{\bf APCP}: ] Precipitation rate $[mm/s]$
\item[{\bf Temp}: ] Air temperature $[K]$
\item[{\bf UGRD}: ] West-to-East or U-component of wind $[m/s]$
\item[{\bf VGRD}: ] South-to-North or V-component of wind $[m/s]$
\item[{\bf Press}: ] Atmospheric Pressure $[pa]$
\item[{\bf SPFH}: ] Water-vapor specific humidity $[kg/kg]$
\label{clm_forcing}
\end{description}
Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.MetFileName                             narr.1hr.txt
\end{verbatim}\end{display}

\pfkey{String}{Solver.CLM.MetFilePath}{no default}
{This key specifies defines the location of 1D, 2D or 3D forcing data. For 1D cases, this is the path to a single forcing file ({\em e.g.}  
\code{narr.1hr.txt}). For 2D and 3D cases, this is the path to the directory  
containing all forcing files. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.MetFilePath		"path/to/met/forcing/data/"
\end{verbatim}\end{display}

\pfkey{integer}{Solver.CLM.MetFileNT}{no default}
{This key specifies the number of timesteps per file for 3D forcing data. 
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.MetFileNT	24	
\end{verbatim}\end{display}

%====
% @BH Forcing the vegetation in CLM
%=====
\pfkey{string}{Solver.CLM.ForceVegetation}{False}
{This key specifies whether vegetation should be forced in \code{CLM}. Currently this option only works for 1D and 3D forcings, as specified by the key \code{Solver.CLM.MetForcing}. Choices for this key include {\bf True} and {\bf False}. Forced vegetation variables are :
\begin{description}
\item[{\bf LAI}: ] Leaf Area Index $[-]$
\item[{\bf SAI}: ] Stem Area Index $[-]$
\item[{\bf Z0M}: ] Aerodynamic roughness length $[m]$
\item[{\bf DISPLA}: ] Displacement height $[m]$
\label{clm_forcing}
\end{description}
In the case of 1D meteorological forcings, \code{CLM} requires four files for vegetation time series and one vegetation map. The four files should be named respectively \code{lai.dat}, \code{sai.dat}, \code{z0m.dat}, \code{displa.dat}. They are ASCII files and contain 18 time-series columns (one per IGBP vegetation class, and each timestep per row). The vegetation map should be a properly distributed 2D \parflow{} binary file (\code{.pfb}) which contains vegetation indices (from 1 to 18). The vegetation map filename is \code{veg\_map.pfb}. \parflow{} uses the vegetation map to pass to \code{CLM} a 2D map for each vegetation variable at each time step. 
In the case of 3D meteorological forcings, \parflow{} expects four distincts properly distributed \parflow{} binary file (\code{.pfb}), the third dimension being the timesteps. The files should be named \code{LAI.pfb}, \code{SAI.pfb}, \code{Z0M.pfb}, \code{DISPLA.pfb}. No vegetation map is needed in this case.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.ForceVegetation  True
\end{verbatim}\end{display}

\pfkey{string}{Solver.WriteSiloCLM}{False}
{This key specifies whether the \code{CLM} writes two dimensional binary output files to a silo binary format. This data may be read in by VisIT and other visualization packages.  Note that \code{CLM} and silo must be compiled and linked at runtime for this option to be active. These files are all written according to the standard format used for all \parflow{} variables, using the {\em runname}, and {\em istep}.  Variables are either two-dimensional or over the number of \code{CLM} layers (default of ten). 
}
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloCLM True
\end{verbatim}\end{display}
The output variables are:
\begin{description}
\item \file{eflx_lh_tot} for latent heat flux total $[W/m^2]$ using the silo variable {\em LatentHeat};
\item \file{eflx_lwrad_out} for outgoing long-wave radiation $[W/m^2]$ using the silo variable {\em LongWave};
\item \file{eflx_sh_tot} for sensible heat flux total $[W/m^2]$ using the silo variable {\em SensibleHeat};
\item \file{eflx_soil_grnd} for ground heat flux $[W/m^2]$ using the silo variable {\em GroundHeat};
\item \file{qflx_evap_tot} for total evaporation $[mm/s]$ using the silo variable {\em EvaporationTotal}; 
\item \file{qflx_evap_grnd} for ground evaporation without condensation $[mm/s]$ using the silo variable {\em EvaporationGroundNoSublimation};
\item \file{qflx_evap_soi} for soil evaporation $[mm/s]$ using the silo variable {\em EvaporationGround};
\item \file{qflx_evap_veg} for vegetation evaporation $[mm/s]$ using the silo variable {\em EvaporationCanopy};
\item \file{qflx_tran_veg} for vegetation transpiration $[mm/s]$ using the silo variable {\em Transpiration}; 
\item \file{qflx_infl} for soil infiltration $[mm/s]$ using the silo variable {\em Infiltration};
\item \file{swe_out} for snow water equivalent $[mm]$ using the silo variable {\em SWE};
\item \file{t_grnd} for ground surface temperature $[K]$ using the silo variable {\em TemperatureGround}; and 
\item \file{t_soil} for soil temperature over all layers $[K]$ using the silo variable {\em TemperatureSoil}.
\end{description}

\pfkey{string}{Solver.PrintCLM}{False}
{This key specifies whether the \code{CLM} writes two dimensional binary output files to a \file{PFB} binary format. Note that \code{CLM} must be compiled and linked at runtime for this option to be active. These files are all written according to the standard format used for all \parflow{} variables, using the {\em runname}, and {\em istep}.  Variables are either two-dimensional or over the number of \code{CLM} layers (default of ten). 
}
\begin{display}\begin{verbatim}
pfset Solver.PrintCLM True
\end{verbatim}\end{display}
The output variables are:
\begin{description}
\item \file{eflx_lh_tot} for latent heat flux total $[W/m^2]$ using the silo variable {\em LatentHeat};
\item \file{eflx_lwrad_out} for outgoing long-wave radiation $[W/m^2]$ using the silo variable {\em LongWave};
\item \file{eflx_sh_tot} for sensible heat flux total $[W/m^2]$ using the silo variable {\em SensibleHeat};
\item \file{eflx_soil_grnd} for ground heat flux $[W/m^2]$ using the silo variable {\em GroundHeat};
\item \file{qflx_evap_tot} for total evaporation $[mm/s]$ using the silo variable {\em EvaporationTotal}; 
\item \file{qflx_evap_grnd} for ground evaporation without sublimation $[mm/s]$ using the silo variable {\em EvaporationGroundNoSublimation};
\item \file{qflx_evap_soi} for soil evaporation $[mm/s]$ using the silo variable {\em EvaporationGround};
\item \file{qflx_evap_veg} for vegetation evaporation $[mm/s]$ using the silo variable {\em EvaporationCanopy};
\item \file{qflx_tran_veg} for vegetation transpiration $[mm/s]$ using the silo variable {\em Transpiration}; 
\item \file{qflx_infl} for soil infiltration $[mm/s]$ using the silo variable {\em Infiltration};
\item \file{swe_out} for snow water equivalent $[mm]$ using the silo variable {\em SWE};
\item \file{t_grnd} for ground surface temperature $[K]$ using the silo variable {\em TemperatureGround}; and 
\item \file{t_soil} for soil temperature over all layers $[K]$ using the silo variable {\em TemperatureSoil}.
\end{description}

\pfkey{string}{Solver.WriteCLMBinary}{True}
{This key specifies whether the \code{CLM} writes two dimensional binary output files in a generic binary format.  
 Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.WriteCLMBinary False
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.BinaryOutDir}{True}
{This key specifies whether the \code{CLM} writes each set of two dimensional binary output files to a corresponding directory.  
These directories my be created before \parflow{} is run (using the tcl script, for example).  Choices for this key include {\bf True} and {\bf False}. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.BinaryOutDir True
\end{verbatim}\end{display}
These directories are: 
\begin{description}
\item \file{/qflx_top_soil} for soil flux;
\item \file{/qflx_infl} for infiltration;
\item \file{/qflx_evap_grnd} for ground evaporation;
\item \file{/eflx_soil_grnd} for ground heat flux;
\item \file{/qflx_evap_veg} for vegetation evaporation;
\item \file{/eflx_sh_tot} for sensible heat flux;
\item \file{/eflx_lh_tot} for latent heat flux; 
\item \file{/qflx_evap_tot} for total evaporation; 
\item \file{/t_grnd} for ground surface temperature; 
\item \file{/qflx_evap_soi} for soil evaporation; 
\item \file{/qflx_tran_veg} for vegetation transpiration; 
\item \file{/eflx_lwrad_out} for outgoing long-wave radiation;
\item \file{/swe_out} for snow water equivalent; and
\item \file{/diag_out} for diagnostics. 
\end{description}

\pfkey{string}{Solver.CLM.CLMFileDir}{no default}
{This key specifies what directory all output from the \code{CLM} is written to.  This key may be set to \file{"./"} or \file{""} to write output to the \parflow{} run directory.  This directory must be created before \parflow{} is run. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.CLMFileDir "CLM_Output/"
\end{verbatim}\end{display}

\pfkey{integer}{Solver.CLM.CLMDumpInterval}{1}
{This key specifies how often output from the \code{CLM} is written.  This key is in integer multipliers of the \code{CLM} timestep. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.CLMDumpInterval 2
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.EvapBeta}{Linear}
{This key specifies the form of the bare soil evaporation $\beta$ parameter in \code{CLM}.   The valid types for this key are {\bf None, Linear, Cosine}. \begin{description}
\item[{\bf None}: ] No beta formulation, $\beta=1$.
\item[{\bf Linear}: ] $\beta=\frac{\phi S-\phi S_{res}}{\phi-\phi S_{res}}$
\item[{\bf Cosine}: ] $\beta=\frac{1}{2}(1-\cos(\frac{(\phi -\phi S_{res})}{(\phi S-\phi S_{res})}\pi)$
\end{description}
Note that $S_{res}$ is specified by the key \code{Solver.CLM.ResSat} below, that $\beta$ is limited between zero and one and also that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.EvapBeta Linear
\end{verbatim}\end{display}

\pfkey{double}{Solver.CLM.ResSat}{0.1}
{This key specifies the residual saturation for the $\beta$ function in \code{CLM} specified above.  Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.ResSat  0.15
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.VegWaterStress}{Saturation}
{This key specifies the form of the plant water stress function $\beta_t$ parameter in \code{CLM}.   The valid types for this key are {\bf None, Saturation, Pressure}. \begin{description}
\item[{\bf None}: ] No transpiration water stress formulation, $\beta_t=1$.
\item[{\bf Saturation}: ] $\beta_t=\frac{\phi S -\phi S_{wp}}{\phi S_{fc}-\phi S_{wp}}$
\item[{\bf Pressure}: ] $\beta_t=\frac{P - P_{wp}}{P_{fc}-P_{wp}}$
\end{description}
Note that the wilting point, $S_{wp}$ or $p_{wp}$, is specified by the key \code{Solver.CLM.WiltingPoint} below, that the field capacity, $S_{fc}$ or $p_{fc}$, is specified by the key \code{Solver.CLM.FieldCapacity} below, that $\beta_t$ is limited between zero and one and also that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.VegWaterStress  Pressure
\end{verbatim}\end{display}

\pfkey{double}{Solver.CLM.WiltingPoint}{0.1}
{This key specifies the wilting point for the $\beta_t$ function in \code{CLM} specified above. Note that the units for this function are pressure $[m]$ for a {\bf Pressure} formulation and saturation $[-]$ for a {\bf Saturation} formulation. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.WiltingPoint  0.15
\end{verbatim}\end{display}

\pfkey{double}{Solver.CLM.FieldCapacity}{1.0}
{This key specifies the field capacity for the $\beta_t$ function in \code{CLM} specified above. Note that the units for this function are pressure $[m]$ for a {\bf Pressure} formulation and saturation $[-]$ for a {\bf Saturation} formulation. Note that \code{CLM} must be compiled and linked at runtime for this option to be active.
}
\begin{display}\begin{verbatim}
pfset Solver.CLM.FieldCapacity  0.95
\end{verbatim}\end{display}
%====
% @RMM Irrigation Keys here
%=====
\pfkey{string}{Solver.CLM.IrrigationTypes}{none}
{This key specifies the form of the irrigation in \code{CLM}.  The valid types for this key are {\bf none, Spray, Drip, Instant}.
} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.IrrigationTypes Drip
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.IrrigationCycle}{Constant}
{This key specifies the cycle of the irrigation in \code{CLM}.  The valid types for this key are {\bf Constant, Deficit}. Note only {\bf Constant} is currently implemented.  Constant cycle applies irrigation each day from IrrigationStartTime to IrrigationStopTime in GMT.
} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.IrrigationCycle Constant
\end{verbatim}\end{display}

\pfkey{double}{Solver.CLM.IrrigationRate}{no default}
{This key specifies the rate of the irrigation in \code{CLM} in [mm/s].  
} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.IrrigationRate 10.
\end{verbatim}\end{display}

\pfkey{double}{Solver.CLM.IrrigationStartTime}{no default}
{This key specifies the start time of the irrigation in \code{CLM} GMT.  
} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.IrrigationStartTime 0800
\end{verbatim}\end{display}

\pfkey{double}{Solver.CLM.IrrigationStopTime}{no default}
{This key specifies the stop time of the irrigation in \code{CLM} GMT.  
} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.IrrigationStopTime 1200
\end{verbatim}\end{display}

\pfkey{double}{Solver.CLM.IrrigationThreshold}{0.5}
{This key specifies the threshold value for the irrigation in \code{CLM} [-].  
} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.IrrigationThreshold 0.2
\end{verbatim}\end{display}

\pfkey{integer}{Solver.CLM.ReuseCount}{1}
{How many times to reuse a \code{CLM} atmospheric forcing file input. For example timestep=1, 
reuse =1 is normal behavior but reuse=2 and timestep=0.5 subdivides 
the time step using the same \code{CLM} input for both halves instead of needing two files.
This is particually useful for large, distributed runs when the user wants to run \parflow{}
at a smaller timestep than the \code{CLM} forcing.  Forcing files will be re-used and
total fluxes adjusted accordingly without needing duplicate files.} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.ReuseCount      5
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.WriteLogs}{True}
{When {\bf False}, this disables writing of the CLM output log files for each processor. 
For example, in the clm.tcl test case, if this flag is added {\bf False}, 
washita.output.txt.\emph{p} and washita.para.out.dat.\emph{p} (were \emph{p} is the processor \#)
are not created, assuming \emph{washita} is the run name.} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.WriteLogs    False
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.WriteLastRST}{False}
{Controls whether CLM restart files are sequentially written or whether a single
file \emph{restart file name}.00000.\emph{p} is overwritten each time the restart file is output,
where \emph{p} is the processor number. If "True" only one file is written/overwritten and if 
"False" outputs are written more frequently. Compatible with DailyRST and ReuseCount;
for the latter, outputs are written every n steps where n is the value of ReuseCount.} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.WriteLastRST   True
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.DailyRST}{True}
{Controls whether CLM writes daily restart files (default) or at every time step when set 
to False; outputs are numbered according to the istep from ParFlow. If {\bf ReuseCount=n}, 
with n greater than 1, the output will be written every n steps (i.e. it still writes hourly 
restart files if your time step is 0.5 or 0.25, etc...). Fully compatible with {\bf WriteLastRST=False}
so that each daily output is overwritten to time 00000 in \emph{restart file name}.00000.p 
where \emph{p} is the processor number.} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.DailyRST    False
\end{verbatim}\end{display}

\pfkey{string}{Solver.CLM.SingleFile}{False}
{Controls whether \parflow{} writes all \code{CLM} output variables as a single file per time step.
When "True", this combines the output of all the CLM output variables into a special multi-layer 
PFB with the file extension ".C.pfb". The first 13 layers correspond to the 2-D CLM
outputs and the remaining layers are the soil temperatures in each layer. For example, 
a model with 4 soil layers will create a SingleFile CLM output with 17 layers at each time step. The
file pseudo code is given below in \S~\ref{ParFlow Binary Files (.c.pfb)} and the variables and units
are as specified in the multiple \code{PFB} and \code{SILO} formats as above.} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.SingleFile   True
\end{verbatim}\end{display}

\pfkey{integer}{Solver.CLM.RootZoneNZ}{10}
{This key sets the number of soil layers the \parflow{} expects from
  \code{CLM}.  It will allocate and format all the arrays for passing
  variables to and from \code{CLM} accordingly.  This value now sets 
  the \code{CLM} number as well so recompilation is not required anymore.
  Most likely the key \code{Solver.CLM.SoiLayer}, described
  below, will also need to be changed.}
\begin{display}\begin{verbatim}
pfset Solver.CLM.RootZoneNZ      4
\end{verbatim}\end{display}

\pfkey{integer}{Solver.CLM.SoiLayer}{7}
{This key sets the soil layer, and thus the soil depth, that \code{CLM} uses for the seasonal temperature adjustment for all leaf and stem area indices.} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.SoiLayer      4
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================
\section{ParFlow NetCDF4 Parallel I/O}
\label{ParFlow NetCDF4 Parallel I/O}
NetCDF4 parallel I/O is being implemented in ParFlow. As of now only output capability is implemented. Input functionality will be added in later version. Currently user has option of printing 3-D time varying pressure or saturation or both in a single NetCDF file containing multiple time steps. User should configure ParFlow(pfsimulatior part) "- -with-netcdf" option and link the appropriate NetCDF4 library.
Naming convention of output files is analogues to binary file names. Following options are available for NetCDF4 output along with various performance tuning options. User is advised to explore NetCDF4 chunking and ROMIO hints option for better I/O performance.

\textbf{\textit{HDF5 Library version 1.8.16 or higher is required for NetCDF4 parallel I/O}}

\pfkey{integer}{NetCDF.NumStepsPerFile}{ }
{This key sets number of time steps user wishes to output in a NetCDF4 file. Once the time step count increases beyond this number, a new file is automatically created.} 
\begin{display}\begin{verbatim}
pfset NetCDF.NumStepsPerFile    5
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WritePressure}{False}
{This key sets pressure variable to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WritePressure    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteSaturation}{False}
{This key sets saturation variable to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteSaturation    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteMannings}{False}
{This key sets Mannings coefficients to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteMannings	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteSubsurface}{False}
{This key sets subsurface data(permeabilities, porosity, specific storage) to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteSubsurface	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteSlopes}{False}
{This key sets x and y slopes to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteSlopes	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteMask}{False}
{This key sets mask to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteMask	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteDZMultiplier}{False}
{This key sets DZ multipliers to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteDZMultiplier	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteEvapTrans}{False}
{This key sets Evaptrans to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteEvapTrans	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteEvapTransSum}{False}
{This key sets Evaptrans sum to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteEvapTransSum	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteOverlandSum}{False}
{This key sets overland sum to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteOverlandSum	    True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteOverlandBCFlux}{False}
{This key sets overland bc flux to be written in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteOverlandBCFlux	    True
\end{verbatim}\end{display}

\subsection{NetCDF4 Chunking}
Chunking may have significant impact on I/O. If this key is not set, default chunking scheme will be used by NetCDF library. Chunks are hypercube(hyperslab) of any dimension. When chunking is used, chunks are written in single write operation which can reduce access times. For more information on chunking, refer to NetCDF4 user guide.

\pfkey{string}{NetCDF.Chunking}{False}
{This key sets chunking for each time varying 3-D variable in NetCDF4 file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.Chunking    True
\end{verbatim}\end{display}

Following keys are used only when \textbf{NetCDF.Chunking} is set to true. These keys are used to set chunk sizes in x, y and z direction. A typical size of chunk in each direction should be equal to number of grid points in each direction for each processor. e.g. If we are using a grid of 400(x)X400(y)X30(z) with 2-D domain decomposition of 8X8, then each core has 50(x)X50(y)X30(z) grid points. These values can be used to set chunk sizes each direction. For unequal distribution, chunk sizes should as large as largest value of grid points on the processor. e.g. If one processor has grid distribution of 40(x)X40(y)X30(z) and another has 50(x)X50(y)X30(z), the later values should be used to set chunk sizes in each direction.
 
\pfkey{integer}{NetCDF.ChunkX}{None}
{This key sets chunking size in x-direction.} 
\begin{display}\begin{verbatim}
pfset NetCDF.ChunkX    50
\end{verbatim}\end{display}

\pfkey{integer}{NetCDF.ChunkY}{None}
{This key sets chunking size in y-direction.} 
\begin{display}\begin{verbatim}
pfset NetCDF.ChunkY    50
\end{verbatim}\end{display}

\pfkey{integer}{NetCDF.ChunkZ}{None}
{This key sets chunking size in z-direction.} 
\begin{display}\begin{verbatim}
pfset NetCDF.ChunkZ    30
\end{verbatim}\end{display}

\subsection{ROMIO Hints}
ROMIO is a poratable MPI-IO implementation developed at Argonne National Laboratory, USA. Currently it is released as a part of MPICH. ROMIO sets hints to optimize I/O operations for MPI-IO layer through MPI\_Info object. This object is passed on to NetCDF4 while creating a file. ROMIO hints are set in a text file in "key" and "value" pair. \textit{For correct settings contact your HPC site administrator}. As in chunking, ROMIO hints can have significant performance impact on I/O.

\pfkey{string}{NetCDF.ROMIOhints}{None}
{This key sets ROMIO hints file to be passed on to NetCDF4 interface.If this key is set, the file must be present and readable in experiment directory.} 
\begin{display}\begin{verbatim}
pfset NetCDF.ROMIOhints    romio.hints
\end{verbatim}\end{display}
An example ROMIO hints file looks as follows.
\begin{display}\begin{verbatim}
romio_ds_write disable
romio_ds_read disable
romio_cb_write enable
romio_cb_read enable
cb_buffer_size 33554432
\end{verbatim}\end{display}

\subsection{Node Level Collective I/O}
A node level collective strategy has been implemented for I/O. One process on each compute node gathers the data, indices and counts from the participating processes on same compute node. All the root processes from each compute node open a parallel NetCDF4 file and write the data. e.g. If ParFlow is running on 3 compute nodes where each node consists of 24 processors(cores); only 3 I/O streams to filesystem would be opened by each root processor each compute node. This strategy could be particularly useful when ParFlow is running on large number of processors and every processor participating in I/O may create a bottleneck.
\textit{\textbf{Node level collective I/O is currently implemented for 2-D domain decomposition and variables Pressure and Saturation only. All the other ParFlow NetCDF output Tcl flags should be set to false(default value). CLM output is independently handled and not affected by this key.  Moreover on speciality architectures, this may not be a portable feature. Users are advised to test this feature on their machine before putting into production.}}

\pfkey{string}{NetCDF.NodeLevelIO}{False}
{This key sets flag for node level collective I/O.} 
\begin{display}\begin{verbatim}
pfset NetCDF.NodeLevelIO   True
\end{verbatim}\end{display}

\subsection{NetCDF4 Initial Conditions: Pressure}
Analogues to ParFlow binary files, NetCDF4 based option can be used to set the initial conditions for pressure to be read from an ``nc" file containing single time step of pressure. The name of the variable in ``nc" file should be ``pressure". A sample NetCDF header of an initial condition file looks as follows. The names of the dimensions are not important. The order of dimensions is important e.g. \textit{(time, lev, lat, lon) or (time,z, y, x)}
\begin{display}\begin{verbatim}
netcdf initial_condition {
dimensions:
	x = 200 ;
	y = 200 ;
	z = 40 ;
	time = UNLIMITED ; // (1 currently)
variables:
	double time(time) ;
	double pressure(time, z, y, x) ;
}
\end{verbatim}\end{display}
\textit{\textbf{Node level collective I/O is currently not implemented for setting initial conditions.}}

\pfkey{string}{ICPressure.Type}{no default}
{This key sets flag for initial conditions to be read from a NetCDF file.} 
\begin{display}\begin{verbatim}
pfset ICPressure.Type   NCFile
pfset Geom.domain.ICPressure.FileName   "initial_condition.nc"
\end{verbatim}\end{display}

\subsection{NetCDF4 Slopes}
NetCDF4 based option can be used slopes to be read from an ``nc" file containing single time step of slope values. The name of the variable in ``nc" file should be ``slopex" and ``slopey" A sample NetCDF header of slope file looks as follows. The names of the dimensions are not important. The order of dimensions is important e.g. \textit{(time, lat, lon) or (time, y, x)}
\begin{display}\begin{verbatim}
netcdf slopex {
dimensions:
	time = UNLIMITED ; // (1 currently)
	lon = 41 ;
	lat = 41 ;
variables:
  	double time(time) ;
	double slopex(time, lat, lon) ;
}
netcdf slopey {
dimensions:
	time = UNLIMITED ; // (1 currently)
	lon = 41 ;
	lat = 41 ;
variables:
	double time(time) ;
	double slopey(time, lat, lon) ;
}
\end{verbatim}\end{display}
The two NetCDF files can be merged into one single file and can be used with tcl flags. The variable names should be exactly as mentioned above. Please refer to ``slopes.nc" under Little Washita test case.
\textit{\textbf{Node level collective I/O is currently not implemented for setting initial conditions.}}

\pfkey{string}{TopoSlopesX.Type}{no default}
{This key sets flag for slopes in x direction to be read from a NetCDF file.} 
\begin{display}\begin{verbatim}
pfset TopoSlopesX.Type   NCFile
pfset TopoSlopesX.FileName   "slopex.nc"
\end{verbatim}\end{display}
\pfkey{string}{TopoSlopesY.Type}{no default}
{This key sets flag for slopes in y direction to be read from a NetCDF file.} 
\begin{display}\begin{verbatim}
pfset TopoSlopesY.Type   NCFile
pfset TopoSlopesy.FileName   "slopey.nc"
\end{verbatim}\end{display}

\subsection{NetCDF4 Transient EvapTrans Forcing}
Following keys can be used for NetCDF4 based transient evaptrans forcing. The file should contain forcing for all time steps. For a given time step, if the forcing is null, zero values could be filled for the given time step in the ``.nc" file. The format of the sample file looks as follows. The names of the dimensions are not important. The order of dimensions is important e.g. \textit{(time, lev, lat, lon) or (time,z, y, x)}
\begin{display}\begin{verbatim}
netcdf evap_trans {
dimensions:
	time = UNLIMITED ; // (1000 currently)
	x = 72 ;
	y = 72 ;
	z = 3 ;
variables:
	double evaptrans(time, z, y, x) ;
}
\end{verbatim}\end{display}
\textit{\textbf{Node level collective I/O is currently not implemented for transient evaptrans forcing.}}

\pfkey{string}{NetCDF.EvapTransFileTransient}{False}
{This key sets flag for transient evaptrans forcing to be read from a NetCDF file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.EvapTransFileTransient   True
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.EvapTrans.FileName}{no default}
{This key sets the name of the NetCDF transient evaptrans forcing file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.EvapTrans.FileName         "evap_trans.nc"
\end{verbatim}\end{display}

\subsection{NetCDF4 CLM Output}
Similar to ParFlow binary and silo, following keys can be used to write output CLM variables in a single NetCDF file containing multiple time steps.

\pfkey{integer}{NetCDF.CLMNumStepsPerFile}{None}
{This key sets number of time steps to be written to a single NetCDF file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.CLMNumStepsPerFile 24
\end{verbatim}\end{display}

\pfkey{string}{NetCDF.WriteCLM}{False}
{This key sets CLM variables to be written in a NetCDF file.} 
\begin{display}\begin{verbatim}
pfset NetCDF.WriteCLM         True
\end{verbatim}\end{display}
The output variables are:
\begin{description}
\item \file{eflx_lh_tot} for latent heat flux total $[W/m^2]$ using the silo variable {\em LatentHeat};
\item \file{eflx_lwrad_out} for outgoing long-wave radiation $[W/m^2]$ using the silo variable {\em LongWave};
\item \file{eflx_sh_tot} for sensible heat flux total $[W/m^2]$ using the silo variable {\em SensibleHeat};
\item \file{eflx_soil_grnd} for ground heat flux $[W/m^2]$ using the silo variable {\em GroundHeat};
\item \file{qflx_evap_tot} for total evaporation $[mm/s]$ using the silo variable {\em EvaporationTotal}; 
\item \file{qflx_evap_grnd} for ground evaporation without condensation $[mm/s]$ using the silo variable {\em EvaporationGroundNoSublimation};
\item \file{qflx_evap_soi} for soil evaporation $[mm/s]$ using the silo variable {\em EvaporationGround};
\item \file{qflx_evap_veg} for vegetation evaporation $[mm/s]$ using the silo variable {\em EvaporationCanopy};
\item \file{qflx_tran_veg} for vegetation transpiration $[mm/s]$ using the silo variable {\em Transpiration}; 
\item \file{qflx_infl} for soil infiltration $[mm/s]$ using the silo variable {\em Infiltration};
\item \file{swe_out} for snow water equivalent $[mm]$ using the silo variable {\em SWE};
\item \file{t_grnd} for ground surface temperature $[K]$ using the silo variable {\em TemperatureGround}; and 
\item \file{t_soil} for soil temperature over all layers $[K]$ using the silo variable {\em TemperatureSoil}.
\end{description}

\subsection{NetCDF4 CLM Input/Forcing}
NetCDF based meteorological forcing can be used with following TCL keys. It is built similar to 2D forcing case for CLM with parflow binary files. All the required forcing variables must be present in one single NetCDF file spanning entire length of simulation. If the simulation ends before number of time steps in NetCDF forcing file, next cycle of simulation can be restarted with same forcing file provided it covers the time span of this cycle.\\
e.g. If the NetCDF forcing file contains 100 time steps and simulation CLM-ParFlow simulation runs for 10 cycles containing 10 time steps in each cycle, the same forcing file can be reused. The user has to set correct value for the key \texttt{Solver.CLM.IstepStart}\\
The format of input file looks as follows. The variable names should match exactly as follows. The names of the dimensions are not important. The order of dimensions is important e.g. \textit{(time, lev, lat, lon) or (time,z, y, x)}
\begin{display}\begin{verbatim}
netcdf metForcing {
dimensions:
	lon = 41 ;
	lat = 41 ;
	time = UNLIMITED ; // (72 currently)
variables:
	double time(time) ;
	double APCP(time, lat, lon) ;
	double DLWR(time, lat, lon) ;
	double DSWR(time, lat, lon) ;
	double Press(time, lat, lon) ;
	double SPFH(time, lat, lon) ;
	double Temp(time, lat, lon) ;
	double UGRD(time, lat, lon) ;
	double VGRD(time, lat, lon) ;
\end{verbatim}\end{display}
\textbf{\textit{Note: While using NetCDF based CLM forcing, \texttt{Solver.CLM.MetFileNT} should be set to its default value of 1}}

\pfkey{string}{Solver.CLM.MetForcing}{no default}
{This key sets meteorological forcing to be read from NetCDF file.} 
\begin{display}\begin{verbatim}
pfset Solver.CLM.MetForcing     NC
\end{verbatim}\end{display}
Set the name of the input/forcing file as follows.
\begin{display}\begin{verbatim}
pfset Solver.CLM.MetFileName   "metForcing.nc"
\end{verbatim}\end{display}
This file should be present in experiment directory. User may create soft links in experiment directory in case where data can not be moved.

\subsection{NetCDF Testing Little Washita Test Case}
The basic NetCDF functionality of output (pressure and saturation) and initial conditions (pressure) can be tested with following tcl script. CLM input/output functionality can also be tested with this case.
\begin{display}\begin{verbatim}
parflow/test/washita/tcl_scripts/LW_NetCDF_Test.tcl
\end{verbatim}\end{display}
This test case will be initialized with following initial condition file, slopes and meteorological forcing.
\begin{display}\begin{verbatim}
parflow/test/washita/parflow_input/press.init.nc
parflow/test/washita/parflow_input/slopes.nc
parflow/test/washita/clm_input/metForcing.nc
\end{verbatim}\end{display}

\section{ParFlow Binary Files (.pfb)}
\label{ParFlow Binary Files (.pfb)}

The \file{.pfb} file format is a binary file format which is used
to store \parflow{} grid data.  It is written as BIG ENDIAN binary bit ordering \cite{endian}.
The format for the file is:

\begin{display}\begin{verbatim}
<double : X>    <double : Y>    <double : Z>
<integer : NX>  <integer : NY>  <integer : NZ>
<double : DX>   <double : DY>   <double : DZ>

<integer : num_subgrids>
FOR subgrid = 0 TO <num_subgrids> - 1
BEGIN
   <integer : ix>  <integer : iy>  <integer : iz>
   <integer : nx>  <integer : ny>  <integer : nz>
   <integer : rx>  <integer : ry>  <integer : rz>
   FOR k = iz TO iz + <nz> - 1
   BEGIN
      FOR j = iy TO iy + <ny> - 1
      BEGIN
         FOR i = ix TO ix + <nx> - 1
         BEGIN
            <double : data_ijk>
         END
      END
   END
END
\end{verbatim}\end{display}
%=============================================================================
%=============================================================================

\section{ParFlow CLM Single Output Binary Files (.c.pfb)}
\label{ParFlow Binary Files (.c.pfb)}

The \file{.pfb} file format is a binary file format which is used
to store \code{CLM} output data in a single file.  It is written as BIG ENDIAN binary bit ordering \cite{endian}.
The format for the file is:

\begin{display}\begin{verbatim}
<double : X>    <double : Y>    <double : Z>
<integer : NX>  <integer : NY>  <integer : NZ>
<double : DX>   <double : DY>   <double : DZ>

<integer : num_subgrids>
FOR subgrid = 0 TO <num_subgrids> - 1
BEGIN
   <integer : ix>  <integer : iy>  <integer : iz>
   <integer : nx>  <integer : ny>  <integer : nz>
   <integer : rx>  <integer : ry>  <integer : rz>
      FOR j = iy TO iy + <ny> - 1
      BEGIN
         FOR i = ix TO ix + <nx> - 1
         BEGIN
            eflx_lh_tot_ij
	    eflx_lwrad_out_ij
	    eflx_sh_tot_ij
	    eflx_soil_grnd_ij
	    qflx_evap_tot_ij
	    qflx_evap_grnd_ij
	    qflx_evap_soi_ij
	    qflx_evap_veg_ij
	    qflx_infl_ij
	    swe_out_ij
	    t_grnd_ij
     IF (clm_irr_type == 1)  qflx_qirr_ij 
ELSE IF (clm_irr_type == 3)  qflx_qirr_inst_ij
ELSE                         NULL
	    FOR k = 1 TO clm_nz
	    tsoil_ijk
	    END
         END
      END
END
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\section{ParFlow Scattered Binary Files (.pfsb)}
\label{ParFlow Scattered Binary Files (.pfsb)}

The \file{.pfsb} file format is a binary file format which is used
to store \parflow{} grid data.
This format is used when the grid data is ``scattered'', that is,
when most of the data is 0.
For data of this type, the \file{.pfsb} file format can reduce
storage requirements considerably.
The format for the file is:

\begin{display}\begin{verbatim}
<double : X>    <double : Y>    <double : Z>
<integer : NX>  <integer : NY>  <integer : NZ>
<double : DX>   <double : DY>   <double : DZ>

<integer : num_subgrids>
FOR subgrid = 0 TO <num_subgrids> - 1
BEGIN
   <integer : ix>  <integer : iy>  <integer : iz>
   <integer : nx>  <integer : ny>  <integer : nz>
   <integer : rx>  <integer : ry>  <integer : rz>
   <integer : num_nonzero_data>
   FOR k = iz TO iz + <nz> - 1
   BEGIN
      FOR j = iy TO iy + <ny> - 1
      BEGIN
         FOR i = ix TO ix + <nx> - 1
         BEGIN
            IF (<data_ijk> > tolerance)
            BEGIN
               <integer : i>  <integer : j>  <integer : k>
               <double : data_ijk>
            END
         END
      END
   END
END
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================

\section{ParFlow Solid Files (.pfsol)}
\label{ParFlow Solid Files (.pfsol)}

The \file{.pfsol} file format is an ASCII file format which is
used to define 3D solids.
The solids are represented by closed triangulated surfaces,
and surface ``patches'' may be associated with each solid.

Note that unlike the user input files, the solid file cannot contain comment
lines. 

The format for the file is:

\begin{display}\begin{verbatim}
<integer : file_version_number>

<integer : num_vertices>
# Vertices
FOR vertex = 0 TO <num_vertices> - 1
BEGIN
   <real : x>  <real : y>  <real : z>
END

# Solids
<integer : num_solids>
FOR solid = 0 TO <num_solids> - 1
BEGIN
   #Triangles
   <integer : num_triangles>
   FOR triangle = 0 TO <num_triangles> - 1
   BEGIN
      <integer : v0> <integer : v1> <integer : v2>
   END

   # Patches
   <integer : num_patches>
   FOR patch = 0 TO <num_patches> - 1
   BEGIN
      <integer : num_patch_triangles>
      FOR patch_triangle = 0 TO <num_patch_triangles> - 1
      BEGIN
         <integer : t>
      END
   END
END
\end{verbatim}\end{display}

\noindent
The field \code{<file_version_number>} is used to make file format changes
more manageable.
The field \code{<num_vertices>} specifies the number of vertices to follow.
The fields \code{<x>}, \code{<y>}, and \code{<z>} define the coordinate
of a triangle vertex.
The field \code{<num_solids>} specifies the number of solids to follow.
The field \code{<num_triangles>} specifies the number of triangles to follow.
The fields \code{<v0>}, \code{<v1>}, and \code{<v2>} are vertex indexes
that specify the 3 vertices of a triangle.  Note that the vertices for each
triangle MUST be specified in an order that makes the normal vector point
outward from the domain.
The field \code{<num_patches>} specifies the number of surface patches
to follow.
The field \code{num_patch_triangles} specifies the number of triangles
indices to follow (these triangles make up the surface patch).
The field \code{<t>} is an index of a triangle on the solid \code{solid}.

\parflow{} \file{.pfsol} files can be created from GMS \file{.sol}
files using the utility \kbd{gmssol2pfsol} located in the
\file{\$PARFLOW_DIR/bin} directory.
This conversion routine takes any number of GMS \file{.sol} files,
concatenates the vertices of the solids defined in the files, throws
away duplicate vertices, then prints out the \file{.pfsol} file.
Information relating the solid index in the resulting \file{.pfsol}
file with the GMS names and material IDs are printed to stdout.

%=============================================================================
%=============================================================================

\section{ParFlow Well Output File (.wells)}
\label{ParFlow Well Output File (.wells)}

A well output file is produced by \parflow{} when wells are defined.
The well output file contains information about the well data being
used in the internal computations and accumulated statistics about
the functioning of the wells.

\noindent
The header section has the following format:
\begin{display}\begin{verbatim}
LINE
BEGIN
   <real : BackgroundX>
   <real : BackgroundY>
   <real : BackgroundZ>
   <integer : BackgroundNX>
   <integer : BackgroundNY>
   <integer : BackgroundNZ>
   <real : BackgroundDX>
   <real : BackgroundDY>
   <real : BackgroundDZ>
END

LINE
BEGIN
   <integer : number_of_phases>
   <integer : number_of_components>
   <integer : number_of_wells>
END

FOR well = 0 TO <number_of_wells> - 1
BEGIN
   LINE
   BEGIN
      <integer : sequence_number>
   END

   LINE
   BEGIN
      <string : well_name>
   END

   LINE
   BEGIN
      <real : well_x_lower>
      <real : well_y_lower>
      <real : well_z_lower>
      <real : well_x_upper>
      <real : well_y_upper>
      <real : well_z_upper>
      <real : well_diameter>
   END

   LINE
   BEGIN
     <integer : well_type>
     <integer : well_action>
   END
END
\end{verbatim}\end{display}

\noindent
The data section has the following format:
\begin{display}\begin{verbatim}
FOR time = 1 TO <number_of_time_intervals>
BEGIN
   LINE
   BEGIN
      <real : time>
   END

   FOR well = 0 TO <number_of_wells> - 1
   BEGIN
      LINE
      BEGIN
         <integer : sequence_number>
      END

      LINE
      BEGIN
         <integer : SubgridIX>
         <integer : SubgridIY>
         <integer : SubgridIZ>
         <integer : SubgridNX>
         <integer : SubgridNY>
         <integer : SubgridNZ>
         <integer : SubgridRX>
         <integer : SubgridRY>
         <integer : SubgridRZ>
      END

      FOR well = 0 TO <number_of_wells> - 1
      BEGIN
         LINE
         BEGIN
            FOR phase = 0 TO <number_of_phases> - 1
            BEGIN
               <real : phase_value>
            END
         END

         IF injection well
         BEGIN
            LINE
            BEGIN
               FOR phase = 0 TO <number_of_phases> - 1
               BEGIN
                  <real : saturation_value>
               END
            END

            LINE
            BEGIN
               FOR phase = 0 TO <number_of_phases> - 1
               BEGIN
                  FOR component = 0 TO <number_of_components> - 1
                  BEGIN
                     <real : component_value>
                  END
               END
            END
         END

         LINE
         BEGIN
            FOR phase = 0 TO <number_of_phases> - 1
            BEGIN
               FOR component = 0 TO <number_of_components> - 1
               BEGIN
                  <real : component_fraction>
               END
            END
         END

         LINE
         BEGIN
            FOR phase = 0 TO <number_of_phases> - 1
            BEGIN
               <real : phase_statistic>
            END
         END

         LINE
         BEGIN
            FOR phase = 0 TO <number_of_phases> - 1
            BEGIN
               <real : saturation_statistic>
            END
         END

         LINE
         BEGIN
            FOR phase = 0 TO <number_of_phases> - 1
            BEGIN
               FOR component = 0 TO <number_of_components> - 1
               BEGIN
                  <real : component_statistic>
               END
            END
         END

         LINE
         BEGIN
            FOR phase = 0 TO <number_of_phases> - 1
            BEGIN
               FOR component = 0 TO <number_of_components> - 1
               BEGIN
                  <real : concentration_data>
               END
            END
         END
      END
   END
END
\end{verbatim}\end{display}
%=============================================================================
%=============================================================================

\section{ParFlow Simple ASCII  and Simple Binary Files (.sa and .sb)}
\label{ParFlow Simple ASCII Files (.sa and .sb)}

The simple binary,\file{.sa}, file format is an ASCII file format which is used by \code{pftools}
to write out \parflow{} grid data.  The simple binary,\code{.sb}, file format is exactly the same, just written as BIG ENDIAN binary bit ordering \cite{endian}.
The format for the file is:

\begin{display}\begin{verbatim}
<integer : NX>  <integer : NY>  <integer : NZ>

   FOR k = 0 TO  <nz> - 1
   BEGIN
      FOR j = 0 TO  <ny> - 1
      BEGIN
         FOR i = 0 TO  <nx> - 1
         BEGIN
            <double : data_ijk>
         END
      END
   END
\end{verbatim}\end{display}


