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

\chapter{The ParFlow System}
\label{The ParFlow System}

The \parflow{} system is still evolving, but here we discuss how to define the problem in
\S~\ref{Defining the Problem}, how to run \parflow{} in
\S~\ref{Running ParFlow}, and restart a simulation in \S~\ref{Restarting a Run}.
We also cover options for visualizing the results in
\S~\ref{Visualizing Output} and summarize the contents of a directory of test problems 
provided with \parflow{} in \S~\ref{Test Directory}.  Finally in \S~\ref{Tutorial}
we walk through two \parflow{} input scripts in detail. 


The reader is also referred to \S~\ref{Manipulating Data}
for a detailed listing of the of functions for manipulating \parflow{} data.  

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

\section{Defining the Problem}
\label{Defining the Problem}

There are many ways to define a problem in \parflow{}, here we summarize the general 
approach for defining a domain (\S~\ref{Defining a domain}) and simulating
 a real watershed (\S~\ref{Defining a Real domain}). 

In all cases the ``main" \parflow{} input file is the \file{.tcl} file.
This input file is a TCL script with some special routines to
create a database which is used as the input for \parflow{}.
See \S~\ref{Main Input File (.tcl)} for details on the format
of this file.
The input values into \parflow{} are defined by a key/value pair.  For
each key you provide the associated value using the \code{pfset} command
inside the input script. 

Since the input file is a TCL script you can use any feature of TCL to
define the problem.  This manual will make no effort to teach TCL so
refer to one of the available TCL manuals for more information
(``Practical Programming in TCL and TK'' by Brent Welch \cite{welch.95} 
is a good starting point).  
This is NOT required, you can get along fine without understanding TCL/TK.

Looking at the example programs in the \file{test} directory (\S~\ref{Test Directory}) and 
going through the annotated input scripts included in this manual ( \S~\ref{Tutorial}) is 
one of the best ways to understand what a \parflow{} input file looks like. 

\subsection{Basic Domain Definition}
\label{Defining a domain}

ParFlow can handle complex geometries and defining the problem may involve several steps.
Users can specify simple box domains directly in the \file{tcl} script. If a more complicated
domain is required, the user may convert geometries into the \file{.pfsol} file format 
(\S~\ref{ParFlow Solid Files (.pfsol)}) using the appropriate \pftools{} conversion utility 
(\S~\ref{Manipulating Data}).  Alternatively, the topography can be specified using 
\file{.pfb} files of the slopes in the x and y directions.

Regardless of the approach the user must set the computational grid within the \file{.pfb} 
script as follows:

\begin{display}\begin{verbatim}
#-----------------------------------------------------------------------------
# Computational Grid
#-----------------------------------------------------------------------------
pfset ComputationalGrid.Lower.X                -10.0
pfset ComputationalGrid.Lower.Y                 10.0
pfset ComputationalGrid.Lower.Z                  1.0
 
pfset ComputationalGrid.DX                       8.89
pfset ComputationalGrid.DY                      10.67
pfset ComputationalGrid.DZ                       1.0
 
pfset ComputationalGrid.NX                      18
pfset ComputationalGrid.NY                      15
pfset ComputationalGrid.NZ                       8

\end{verbatim}\end{display}

The value is normally a single string, double, or integer.  In some cases,
in particular for a list of names, you need to supply a space seperated
sequence.  This can be done using either a double quote or braces.

\begin{display}\begin{verbatim}
pfset Geom.domain.Patches "left right front back bottom top"

pfset Geom.domain.Patches {left right front back bottom top}
\end{verbatim}\end{display}

For commands longer than a single line, the TCL continuation character can be
used, 
\begin{display}\begin{verbatim}
pfset Geom.domain.Patches "very_long_name_1 very_long_name_2 very_long_name_3 \
                           very_long_name_4 very_long_name_5 very_long_name_6"
\end{verbatim}\end{display}

%=============================================================================
%=============================================================================
\subsection{Setting Up a Real Domain}
\label{Defining a Real domain}

This section provides a brief outline of a sample workflow for setup \parflow{} \code{CLM} 
simulation of a real domain. Of course there are many ways to accomplish this and users
are encouraged to develop a workflow that works for them. 

This example assumes that you are running with \parflow{} \code{CLM} and it uses 
slope files and an indicator file to define the topography and geologic units of the domain. 
An alternate approach would be to define geometries by building a \file{.pfsol} file 
(\S~\ref{ParFlow Solid Files (.pfsol)}) using the appropriate \pftools{} conversion 
utility(\S~\ref{Manipulating Data}).

The general approach is as follows:

\begin{enumerate}
   \item Gather input datasets to define the domain. First decide the resolution that you 
   would like to simulate at. Then gather the following datasets at the appropriate resolution 
   for your domain:
   \begin{enumerate}
		\item Elevation (DEM)
	    \item Soil data for the near surface layers 
	    \item Geologic maps for the deeper subsurface
	    \item Land Cover
	\end{enumerate}
	
   \item Create consistent gridded layers that are all clipped to your domain and have 
  the same number of grid cells

   \item Convert gridded files to \file{.pfb} (\S~\ref{ParFlow Binary Files (.pfb)}). One way 
   to accomplish this is by reformatting the gridded outputs to the correct \parflow{}
   \file{.sa} order (\S~\ref{ParFlow Simple ASCII Files (.sa and .sb)}) and to convert
   the \file{.sa} file to \file{.pfb} using the conversion tools (see \S~\ref{common_pftcl}
   Example 1). If you have an elevation file in \file{.pfb} format, you may wish to preserve
   it as provenance for the slopes and for use in post-processing tools. You may point
   ParFlow to the elevation:
   \begin{display}\begin{verbatim}
   pfset TopoSlopes.Elevation.FileName "elevation.pfb"
   \end{verbatim}\end{display}
  
   \item Calculate slopes in the x and y directions from the elevation dataset. This can be 
   done with the built in tools as shown in \S~\ref{common_pftcl} Example 5. 
   In most cases some additional processing of  the DEM will be required to ensure that 
   the drainage patterns are correct. To check this  you  can run a ``parking lot test" by 
   setting the permeability of surface to almost zero and adding a flux to the top surface.
   If the results from this test don't look right (i.e. your runoff patterns don't match 
   what you expect) you will need to go back and modify your DEM.  The built in \parflow{}
   tools pitfill and flatfill can be used to address some issues. 
   (These tools are also shown in \S~\ref{common_pftcl} Example 5). 

   \item Create an indicator file for the subsurface. The indicator file is a 3D \file{.pfb} 
   file with the same dimensions as your domain that has an integer for every cell 
   designating which unit it belongs to. The units you define will correspond to the soil 
   types and geologic units from your input datasets. 

   \item Determine the hydrologic properties for each of the subsurface units defined in 
   the indicator file. You will need: Permeability, specific storage, porosity and 
   vanGenuchten parameters. 

   \item At this point you are ready to run a \parflow{} model without \code{CLM} and if you 
   don't need to include the land surface model in your simulations you can ignore the 
   following steps. Either way, at this point it is advisable to run a ``spinup" simulation 
   to initialize the water table. There are several ways to approach this. One way is to 
   start with the water table at a constant depth and run for 
   a long time with a constant recharge forcing until the water table reaches a steady 
   state. There are some additional key for spinup runs that are provided in \S~\ref{Spinup Options}.

   \item Convert land cover classifications to the IGBP\footnote{http://www.igbp.net}
   land cover classes that are used in CLM. 
   
   \begin{itemize}
		 \item[] 1. Evergreen Needleleaf Forest
		 \item[] 2. Evergreen Broadleaf Forest
		 \item[] 3. Deciduous Needleleaf Forest
		 \item[] 4. Deciduous Broadleaf Forest
		 \item[] 5. Mixed Forests
		 \item[] 6. Closed Shrublands
		 \item[] 7. Open Shrublands
		 \item[] 8. Woody Savannas
		 \item[] 9. Savannas
     	 \item[] 10. Grasslands
		 \item[] 11. Permanent Wetlands
		 \item[] 12. Croplands
		 \item[] 13. Urban and Built-Up
		 \item[] 14. Cropland/Natural Vegetation Mosaic
		 \item[] 15. Snow and Ice
		 \item[] 16. Barren or Sparsely Vegetated
		 \item[] 17. Water
		 \item[] 18. Wooded Tundra
  \end{itemize}

   \item Create a \code{CLM} vegm file that designates the land cover fractions for every 
   cell (Refer to the \file{clm input} directory in the Washita Example for an sample of 
   what a \file{vegm} file should look like).

   \item Create a \code{CLM} driver file to set the parameters for the CLM model (Refer to the 
   \file{clm input} directory in the Washita Example for a sample of a \code{CLM} driver file).

   \item Assemble meteorological forcing data for your domain. CLM uses Greenwich Mean Time (GMT), not local time. The year, date and hour (in GMT) that the forcing begins should match the values in \code{drv_clmin.dat}.  CLM requires the following 
   variables (also described on p. \pageref{clm_forcing}): 
   
   \begin{itemize}
		\item DSWR: Visible or short-wave radiation $[W/m^2]$.
		\item DLWR: Long wave radiation $[W/m^2]$
		\item APCP: Precipitation $[mm/s]$
		\item Temp: Air Temperature $[K]$
		\item UGRD: East-west wind speed $[m/s]$
                \item VGRD: South-to-North wind speed $[m/s]$
		\item Press: Atmospheric pressure $[pa]$
		\item SPFH:	Specific humidity $[kg/kg]$
	\end{itemize}
	
	If you choose to do spatially heterogenous forcings you will need to generate separate 
	files for each variable. The files should be formatted in the standard ParFlow format 
	with the third (i.e. z dimension) as time. If you are doing hourly simulations it is 
	standard practice to put 24 hours in one file, but you can decide how many time steps per 
	file. For an example of heterogenous forcing files refer to the \file{NLDAS} directory
	in the Washita Example)

	Alternatively, if you would like to force the model with spatially homogenous forcings, 
	then a single file can be provided where each variable is a column and rows designate 
	time steps. 

   \item Run your simulation!
   
\end{enumerate}

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

\section{Running ParFlow}
\label{Running ParFlow}

Once the problem input is defined, you need to add a few things to 
the script to make it execute \parflow{}.  First you need to add
the TCL commands to load the \parflow{} command package.

\begin{display}\begin{verbatim}
#
# Import the ParFlow TCL package
#
lappend auto_path $env(PARFLOW_DIR)/bin 
package require parflow
namespace import Parflow::*
\end{verbatim}\end{display}

This loads the \code{pfset} and other \parflow{} commands into the
TCL shell.

Since this is a script you need to actually run \parflow{}. These are
normally the last lines of the input script.

\begin{display}\begin{verbatim}
#-----------------------------------------------------------------------------
# Run and Unload the ParFlow output files
#-----------------------------------------------------------------------------
pfrun default_single
pfundist default_single
\end{verbatim}\end{display}


The \code{pfrun} command runs \parflow{} with the database as it
exists at that point in the file.  The argument is the name to give to
the output files (which will normally be the same as the name of the
script).  Advanced users can set up multiple problems within the input
script by using different output names.  

The \code{pfundist} command takes the output files from the \parflow{}
run and undistributes them.  \parflow{} uses a virtual file system
which allows files to be distributed across the processors.  The
\code{pfundist} takes these files and collapses them into a single
file.  On some machines if you don't do the \code{pfundist} you will
see many files after the run.  Each of these contains the output from
a single node; before attempting using them you should undistribute them.

Since the input file is a TCL script run it using TCL:

\begin{display}\begin{verbatim}
tclsh runname.tcl
\end{verbatim}\end{display}

NOTE: Make sure you are using TCL 8.0 or later.  The script will not
work with earlier releases.

One output file of particular interest is the \file{<run
name>.out.log} file.  This file contains information about the run
such as number of processes used, convergence history of algorithms,
timings and MFLOP rates.  For Richards' equation problems (including overland flow) the \file{<run
name>.out.kinsol.log} file contains the nonlinear convergence information for each timestep.  Additionally, the \file{<run
name>.out.tx} contains all information routed to \file{standard out} of the machine you are running on and often contains error messages and other control information.

%=============================================================================
%=============================================================================
\section{Restarting a Run}
\label{Restarting a Run}

A \parflow{} run may need to be restarted because either a
system time limit has been reached, \parflow{} has been prematurely
terminated or the user specifically sets up a problem to run in segments.
In order to restart a run the user needs to know the conditions under which
\parflow{} stopped.  If \parflow{} was prematurely terminated then the
user must examine the output files from the last ``timed dump" to see if
they are complete.  If not then those data files should be discarded and the
output files from the next to last ``timed dump" will be used in the
restarting procedure.  As an important note, if any set of ``timed dump"
files are deleted remember to also delete corresponding lines in the well
output file or recombining the well output files from the individual segments
afterwards will be difficult.  It is not necessary to delete lines from
the log file as you will only be noting information from it.  To summarize,
make sure all the important output data files are complete, accurate and
consistent with each other.

Given a set of complete, consistent output files - to restart a run follow
this procedure :

\begin{enumerate}
   \item Note the important information for restarting :
   \begin{itemize}
      \item Write down the dump sequence number for the last collection of
           ``timed dump'' data.
      \item Examine the log file to find out what real time that ``timed dump"
            data was written out at and write it down.
   \end{itemize}
   \item Prepare input data files from output data files :
   \begin{itemize}
      \item Take the last pressure output file before the restart with the sequence number
            from above and format them for regular input using the keys detailed in \S~\ref{Initial Conditions: Pressure} and possibly the \code{pfdist} utility in the input script.
   \end{itemize}
   \item Change the Main Input File \S~\ref{Main Input File (.tcl)} :
   \begin{itemize}
      \item Edit the .tcl file (you may want to save the old one) and
            utilize the pressure initial condition input file option (as referenced above) to specify
             the input files you created above as initial conditions
            for concentrations.
   \end{itemize}
   \item Restart the run :
   \begin{itemize}
      \item Utilizing an editor recreate all the input parameters used
            in the run except for the following two items :
            \begin{itemize}
               \item Use the dump sequence number from step 1
                     as the start\_count.
               \item Use the real time that the dump occured at from step 1
                     as the start\_time.
	     \item To restart with \code{CLM}, use the \file{Solver.CLM.IstepStart} key described in \S~\ref{CLM Solver Parameters} with a value equal to the dump sequence plus one. Make sure this corresponds to changes to \code{drv_clmin.dat}. 
            \end{itemize}
   \end{itemize}
\end{enumerate}

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

\section{Visualizing Output}
\label{Visualizing Output}

While \parflow{} does not have any visualization capabilities built-in, there are a number
flexible, free options.  Probably the best option is to use \emph{VisIt}.
\emph{VisIt} is a powerful, free, open-source, rendering environment.  It is 
multiplatform and may be downloaded directly from:
\code{https://visit.llnl.gov/}  \\
The most flexible option for using \emph{VisIt} to view \parflow{} output is to write files using the \emph{SILO} format, which is available either as a direct output option (described in \S~\ref{Code Parameters}) or a conversion option using pftools.  Many other output conversion options exist as described in \S~\ref{Manipulating Data} and this allows \parflow{} output to be converted into formats used by almost all visualization software.
%=============================================================================

\section{Directory of Test Cases}
\label{Test Directory}

\parflow{} comes with a directory containing a few simple input files for use
as templates in making new files and for use in testing the code.  These files sit in the \file{/test} directory described earlier. This section gives a brief description of the problems in this directory.

\begin{description}

\item{\file{crater2D.tcl} An example of a two-dimensional, variably-saturated crater infiltration problem with time-varying boundary conditions.  It uses the solid file \file{crater2D.pfsol.}}

\item{\file{default_richards.tcl} The default variably-saturated Richards' Equation simulation test script.}

\item{\file{default_single.tcl} The default parflow, single-processor, fully-saturated test script.}

\item{\file{forsyth2.tcl} An example two-dimensional, variably-saturated infiltration problem with layers of different hydraulic properties. It runs problem 2 in \cite{FWP95} and uses the solid file \file{fors2_hf.pfsol}.}

\item{\file{harvey.flow.tcl} An example from \cite{MWH07} for the Cape Cod bacterial injection site.  This example is a three-dimensional, fully-saturated flow problem with spatially heterogeneous media (using a correlated, random field approach).  It also provides examples of how tcl/tk scripts may be used in conjunction with ParFlow to loop iteratively or to run other scripts or programs.  It uses the input text file \file{stats4.txt}. This input script is fully detailed in \S~\ref{Tutorial}}

\item{\file{default_overland.tcl} An overland flow boundary condition test and example script based loosely on the V-catchment problem in \cite{KM06}.  There are options provided to expand this problem into other overland flow-type, transient boundary-type problems included in the file as well.}

\item{\file{LW_var_dz_spinup.tcl} An example that uses the Little Washita domain to demonstrate a steady-state spinup initialization using P-E forcing. It also demonstrates the variable dz keys.}

\item{\file {LW_var_dz.tcl} An example that uses the Little Washita domain to demonstrate surface flow network development.  It also uses the variable dz keys.}

\item{\file {Evap_Trans_test.tcl} An example that modifies the \file{default_overland.tcl} to demonstrate steady-state external flux \file{.pfb} files.}

\item{\file {overland_flux.tcl} An example that modifies the \file{default_overland.tcl} to demonstrate transient external flux \file{.pfb} files.}

\item{\file{/clm/clm.tcl} An example of how to use \parflow{} coupled to {\file{clm}}.  This directory also includes {\file{clm}}-specific input. Note: this problem will only run if {\file{--with-clm}} flag is used during the configure and build process.}

\item{\file{water_balance_x.tcl} and {\file{water_balance_y.tcl}}. An overland flow example script that uses the water-balance routines integrated into {\file{pftools}}.  These two problems are based on simple overland flow conditions with slopes primarily in the x or y-directions.  Note: this problem only will run if the Silo file capability is used, that is a {\file{--with-silo=PATH}} flag is used during the configure and build process.}

\item{\file{pfmg.tcl} and \file{pfmg_octree.tcl}. Tests of the external \emph{Hypre} preconditioner options.  Note: this problem only will run if the \emph{Hypre} capability is used, that is a {\file{--with-hypre=PATH}} flag is used during the configure and build process.}

\item{\file{test_x.tcl} A test problem for the Richards' solver that compares output to an analytical solution.}

\item{\file{/washita/tcl_scripts/LW_Test.tcl} A three day simulation of the Little Washita domain using \parflow{} \code{CLM} with 3D forcings. }

\end{description}

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

\section{Annotated Input Scripts}
\label{Tutorial}

This section contains two annotated input scripts:

\begin{itemize}
	\item \S~\ref{Harvey Flow Example} contains the harvey flow example (\file{harvey.flow.tcl})
	which is an idealized domain with a heterogenous subsurface. The example also demonstrates how to generate 
	multiple realizations of the subsurface and add pumping wells. 
	 
	\item \S~\ref{Little Washita Example} contains the Little Washita example (\file{LW_Test.tcl})
	which simulates a moderately sized (41km by 41km) real domain using \parflow{} \code{CLM}
	with 3D meteorological forcings.
	 
\end{itemize}

To run \parflow{}, you use a script written in Tcl/TK.  This script has a lot of flexibility, as it is somewhere in between a program 
and a user interface.  The tcl script gives \parflow{} the data it requires (or tells \parflow{} where to find or read in that data) 
and also tells \parflow{} to run.

To run the simulation:
\begin{enumerate}
\item{make any modifications to the tcl input script (and give a new name, if you want to)}
\item{save the tcl script}
\item{For Linux/Unix/OSX: invoke the script from the command line using the tcl-shell, this looks like: \file{>tclsh filename.tcl}}
\item{Wait patiently for the command prompt to return (Linux/Unix/OSX) indicating that \parflow{} has finished.  Intermediate files are written as the simulation runs, however there is no other indication that \parflow{} is running.}
\end{enumerate}

To modify a tcl script, you right-click and select edit from the menu.  If you select open, you will run the script.

{\bf Note:}  The units for {\bf K} (\i{m/d}, usually) are critical to the entire construction.  
These length and time units for {\bf K} set the units for all other variables (input or generated, 
throughout the entire simulation) in the simulation.  \parflow{} can set to solve using hydraulic 
conductivity by literally setting density, viscosity and gravity to one (as is done in the script below).  
This means the pressure units are in length (meters), so pressure is now so-called pressure-head.

%________________________________
\subsection{Harvey Flow Example}
\label{Harvey Flow Example}

This tutorial matches the \file{harvey_flow.tcl} file found in the \file{/test} directory.  This example is directly from \cite{MWH07}.  
This example demonstrates how to set up and run a fully saturated flow problem with heterogeneous hydraulic conductivity using 
the turning bands approach \cite{TAG89}. Given statistical parameters describing the geology of your site, 
this script can be easily modified to make as many realizations 
of the subsurface as you like, each different and yet having the same statistical parameters, useful for a Monte Carlo simulation.  This example
is the basis for several fully-saturated \parflow{} applications \cite{Siirila12a,Siirila12b,SNSMM10,Atchley13a,Atchley13b,Cui14}.

When the script runs, it creates a new directory named \file{/flow} right in the directory where the tcl script is stored.  \parflow{} then puts all its output in \file{/flow}.  Of course, you can change the name and location of this output directory by modifying the tcl script that runs \parflow{}.


Now for the tcl script:

\begin{verbatim}
#
# Import the ParFlow TCL package
#
\end{verbatim}

These first three lines are what link \parflow{} and the tcl script, thus allowing you to use a set of commands seen later, such as \code{pfset}, etc.

\begin{verbatim}
lappend auto_path $env(PARFLOW_DIR)/bin 
package require parflow
namespace import Parflow::*

#-----------------------------------------------------------------------------
# File input version number
#-----------------------------------------------------------------------------
pfset FileVersion 4
\end{verbatim}

These next lines set the parallel process topology.  The domain is divided in \emph{x},\emph{y} and \emph{z} by \code{P}, \code{Q} and \code{R}.
The total number of processors is \code{P*Q*R}  (\emph{see} \S~\ref{Computing Topology}).

\begin{verbatim}
#----------------------------------------------------------------------------
# Process Topology
#----------------------------------------------------------------------------

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

Next we set up the computational grid (\emph{see} \S~\ref{Defining the Problem} and  \S~\ref{Computational Grid}).

\begin{verbatim}
#----------------------------------------------------------------------------
# Computational Grid
#----------------------------------------------------------------------------
\end{verbatim}

Locate the origin in the domain.

\begin{verbatim}
pfset ComputationalGrid.Lower.X    0.0
pfset ComputationalGrid.Lower.Y    0.0
pfset ComputationalGrid.Lower.Z    0.0
\end{verbatim}

Define the size of the domain grid block.  Length units, same as those on hydraulic conductivity.

\begin{verbatim}
pfset ComputationalGrid.DX     0.34
pfset ComputationalGrid.DY     0.34
pfset ComputationalGrid.DZ     0.038
\end{verbatim}

Define the number of grid blocks in the domain.

\begin{verbatim}
pfset ComputationalGrid.NX      50
pfset ComputationalGrid.NY      30
pfset ComputationalGrid.NZ      100
\end{verbatim}

This next piece is comparable to a pre-declaration of variables.  These will be areas in our domain geometry. 
The regions themselves will be defined later.  You must always have one that is the name of your entire domain.  
If you want subsections within your domain, you may declare these as well.  For Cape Cod, we have the entire domain,
and also the 2 (upper and lower) permeability zones in the aquifer.

\begin{verbatim}
#----------------------------------------------------------------------------
# The Names of the GeomInputs
#----------------------------------------------------------------------------
pfset GeomInput.Names "domain_input upper_aquifer_input lower_aquifer_input"
\end{verbatim}

Now you characterize your domain that you just pre-declared to be a \code{box} (\emph{see} \S~\ref{Geometries}),
and you also give it a name, \code{domain}.

\begin{verbatim}
#----------------------------------------------------------------------------
# Domain Geometry Input
#----------------------------------------------------------------------------
pfset GeomInput.domain_input.InputType       Box
pfset GeomInput.domain_input.GeomName   domain
\end{verbatim}

Here, you set the limits in space for your entire domain.  The span from \file{Lower.X} to \file{Upper.X} 
will be equal to the product of \file{ComputationalGrid.DX} times \file{ComputationalGrid.NX}.  
Same for {\em Y} and {\em Z} (i.e. the number of grid elements times size of the grid element has to
equal the size of the grid in each dimension).   The \code{Patches} key assigns names to the outside edges, 
because the domain is the limit of the problem in space.

\begin{verbatim}
#----------------------------------------------------------------------------
# Domain Geometry
#----------------------------------------------------------------------------
pfset Geom.domain.Lower.X       0.0 
pfset Geom.domain.Lower.Y       0.0
pfset Geom.domain.Lower.Z       0.0

pfset Geom.domain.Upper.X       17.0
pfset Geom.domain.Upper.Y       10.2
pfset Geom.domain.Upper.Z       3.8

pfset Geom.domain.Patches "left right front back bottom top"
\end{verbatim}

Just like domain geometry, you also set the limits in space for the individual components 
(upper and lower, as defined in the Names of GeomInputs pre-declaration).  There are no patches 
for these geometries as they are internal to the domain.

\begin{verbatim}
#----------------------------------------------------------------------------
# Upper Aquifer Geometry Input
#----------------------------------------------------------------------------
pfset GeomInput.upper_aquifer_input.InputType        Box
pfset GeomInput.upper_aquifer_input.GeomName   upper_aquifer

#----------------------------------------------------------------------------
# Upper Aquifer Geometry
#----------------------------------------------------------------------------
pfset Geom.upper_aquifer.Lower.X                        0.0 
pfset Geom.upper_aquifer.Lower.Y                        0.0
pfset Geom.upper_aquifer.Lower.Z                        1.5

pfset Geom.upper_aquifer.Upper.X                        17.0
pfset Geom.upper_aquifer.Upper.Y                        10.2
pfset Geom.upper_aquifer.Upper.Z                        1.5

#----------------------------------------------------------------------------
# Lower Aquifer Geometry Input
#----------------------------------------------------------------------------
pfset GeomInput.lower_aquifer_input.InputType       Box
pfset GeomInput.lower_aquifer_input.GeomName   lower_aquifer

#----------------------------------------------------------------------------
# Lower Aquifer Geometry
#----------------------------------------------------------------------------
pfset Geom.lower_aquifer.Lower.X      0.0 
pfset Geom.lower_aquifer.Lower.Y      0.0
pfset Geom.lower_aquifer.Lower.Z      0.0

pfset Geom.lower_aquifer.Upper.X     17.0
pfset Geom.lower_aquifer.Upper.Y     10.2
pfset Geom.lower_aquifer.Upper.Z      1.5
\end{verbatim}

Now you add permeability data to the domain sections defined above (\S~\ref{Permeability}).  
You can reassign values simply by re-stating them -- there is no need to comment out or 
delete the previous version -- the final statement is the only one that counts.

\begin{verbatim}
#----------------------------------------------------------------------------
# Perm
#----------------------------------------------------------------------------
\end{verbatim}

Name the permeability regions you will describe.

\begin{verbatim}
pfset Geom.Perm.Names "upper_aquifer lower_aquifer"
\end{verbatim}

You can set, for example homogeneous, constant permeability, or you can generate a random 
field that meets your statistical requirements.  To define a constant permeability for the entire domain:

\begin{verbatim}
#pfset Geom.domain.Perm.Type     Constant
#pfset Geom.domain.Perm.Value    4.0
\end{verbatim}

However, for Cape Cod, we did not want a constant permeability field, so we instead generated a 
random permeability field meeting our statistical parameters for each the upper and lower zones.  
Third from the bottom is the \code{Seed}.  This is a random starting point to generate the K field.  
Pick any large {\em ODD} number.  First we do something tricky with { \file Tcl/TK}. 
We use the native commands within tcl to open a text file and read in locally set variables. 
Note we use set here and not pfset.  One is a native tcl command, the other a \parflow{}-specific command.  
For this problem, we are linking the parameter estimation code, PEST to \parflow{}. 
PEST writes out the ascii file \file{stats4.txt} (also located in the \file{/test} directory) 
as the result of a calibration run.  Since we are not coupled to PEST in this example, 
we just read in the file and use the values to assign statistical properties.

\begin{verbatim}
# we open a file, in this case from PEST to set upper and lower # kg and sigma
#
set fileId [open stats4.txt r 0600]
set kgu [gets $fileId]
set varu [gets $fileId]
set kgl [gets $fileId]
set varl [gets $fileId]
close $fileId
\end{verbatim}

Now we set the heterogeneous parameters for the Upper and Lower aquifers (\emph{see} \S~\ref{Permeability}).  
Note the special section at the very end of this block where we reset the geometric mean and standard 
deviation to our values we read in from a file.  {\bf Note:} \parflow{} uses {\em Standard Deviation} not {\em Variance}.

\begin{verbatim}
pfset Geom.upper_aquifer.Perm.Type "TurnBands"
pfset Geom.upper_aquifer.Perm.LambdaX  3.60
pfset Geom.upper_aquifer.Perm.LambdaY  3.60
pfset Geom.upper_aquifer.Perm.LambdaZ  0.19
pfset Geom.upper_aquifer.Perm.GeomMean  112.00

pfset Geom.upper_aquifer.Perm.Sigma   1.0
pfset Geom.upper_aquifer.Perm.Sigma   0.48989794
pfset Geom.upper_aquifer.Perm.NumLines 150
pfset Geom.upper_aquifer.Perm.RZeta  5.0
pfset Geom.upper_aquifer.Perm.KMax  100.0
pfset Geom.upper_aquifer.Perm.DelK  0.2
pfset Geom.upper_aquifer.Perm.Seed  33333
pfset Geom.upper_aquifer.Perm.LogNormal Log
pfset Geom.upper_aquifer.Perm.StratType Bottom
pfset Geom.lower_aquifer.Perm.Type "TurnBands"
pfset Geom.lower_aquifer.Perm.LambdaX  3.60
pfset Geom.lower_aquifer.Perm.LambdaY  3.60
pfset Geom.lower_aquifer.Perm.LambdaZ  0.19

pfset Geom.lower_aquifer.Perm.GeomMean  77.0
pfset Geom.lower_aquifer.Perm.Sigma   1.0
pfset Geom.lower_aquifer.Perm.Sigma   0.48989794
pfset Geom.lower_aquifer.Perm.NumLines 150
pfset Geom.lower_aquifer.Perm.RZeta  5.0
pfset Geom.lower_aquifer.Perm.KMax  100.0
pfset Geom.lower_aquifer.Perm.DelK  0.2
pfset Geom.lower_aquifer.Perm.Seed  33333
pfset Geom.lower_aquifer.Perm.LogNormal Log
pfset Geom.lower_aquifer.Perm.StratType Bottom

#pfset lower aqu and upper aq stats to pest/read in values

pfset Geom.upper_aquifer.Perm.GeomMean  $kgu
pfset Geom.upper_aquifer.Perm.Sigma  $varu

pfset Geom.lower_aquifer.Perm.GeomMean  $kgl
pfset Geom.lower_aquifer.Perm.Sigma  $varl
\end{verbatim}

The following section allows you to specify the permeability tensor.  
In the case below, permeability is symmetric in all directions (x, y, and z) and therefore each is set to 1.0.

\begin{verbatim}
pfset Perm.TensorType               TensorByGeom

pfset Geom.Perm.TensorByGeom.Names  "domain"

pfset Geom.domain.Perm.TensorValX  1.0
pfset Geom.domain.Perm.TensorValY  1.0
pfset Geom.domain.Perm.TensorValZ  1.0
\end{verbatim}

Next we set the specific storage, though this is not used in the IMPES/steady-state calculation.

\begin{verbatim}
#----------------------------------------------------------------------------
# Specific Storage
#----------------------------------------------------------------------------
# specific storage does not figure into the impes (fully sat) 
# case but we still need a key for it

pfset SpecificStorage.Type            Constant
pfset SpecificStorage.GeomNames       ""
pfset Geom.domain.SpecificStorage.Value 1.0e-4
\end{verbatim}

\parflow{} has the capability to deal with a multiphase system, but we only have one (water) at Cape Cod. 
As we stated earlier, we set density and viscosity artificially (and later gravity) both to 1.0.  
Again, this is merely a trick to solve for hydraulic conductivity and pressure head.  
If you were to set density and viscosity to their true values, the code would calculate {\bf k} (permeability). 
By using the \emph{normalized} values instead, you effectively embed the conversion of {\bf k} to {\bf K} (hydraulic conductivity). 
So this way, we get hydraulic conductivity, which is what we want for this problem.

\begin{verbatim}
#----------------------------------------------------------------------------
# Phases
#----------------------------------------------------------------------------

pfset Phase.Names "water"

pfset Phase.water.Density.Type	Constant
pfset Phase.water.Density.Value	1.0

pfset Phase.water.Viscosity.Type	Constant
pfset Phase.water.Viscosity.Value	1.0
\end{verbatim}

We will not use the \parflow{} grid based transport scheme.  We will then leave 
contaminants blank because we will use a different code to model (virus, tracer) contamination.

\begin{verbatim}
#----------------------------------------------------------------------------
# Contaminants
#----------------------------------------------------------------------------
pfset Contaminants.Names			""
\end{verbatim}

As with density and viscosity, gravity is normalized here.  If we used the true value
(in the \emph{[L]} and \emph{[T]} units of hydraulic conductivity) the code would be calculating permeability.  
Instead, we normalize so that the code calculates hydraulic conductivity.

\begin{verbatim}
#----------------------------------------------------------------------------
# Gravity
#----------------------------------------------------------------------------

pfset Gravity				1.0

#----------------------------------------------------------------------------
# Setup timing info
#----------------------------------------------------------------------------
\end{verbatim}

This basic time unit of 1.0 is used for transient boundary and well conditions. 
We are not using those features in this example.

\begin{verbatim} 
pfset TimingInfo.BaseUnit		1.0
\end{verbatim}

Cape Cod is a steady state problem, so these timing features are again unused, but need to be included.

\begin{verbatim}
pfset TimingInfo.StartCount   -1
pfset TimingInfo.StartTime     0.0
pfset TimingInfo.StopTime     0.0
\end{verbatim}

Set the \code{dump interval} to -1 to report info at the end of every calculation, 
which in this case is only when steady state has been reached.

\begin{verbatim}
pfset TimingInfo.DumpInterval	       -1
\end{verbatim}

Next, we assign the porosity (\emph{see} \S~\ref{Porosity}).  For the Cape Cod, the porosity is 0.39.  

\begin{verbatim}
#----------------------------------------------------------------------------
# Porosity
#----------------------------------------------------------------------------

pfset Geom.Porosity.GeomNames          domain

pfset Geom.domain.Porosity.Type    Constant
pfset Geom.domain.Porosity.Value   0.390
\end{verbatim}

Having defined the geometry of our problem before and named it \code{domain}, we are
now ready to report/upload that problem, which we do here.

\begin{verbatim}
#----------------------------------------------------------------------------
# Domain
#----------------------------------------------------------------------------
pfset Domain.GeomName domain
\end{verbatim}

Mobility between phases is set to 1.0 because we only have one phase (water).

\begin{verbatim}
#----------------------------------------------------------------------------
# Mobility
#----------------------------------------------------------------------------
pfset Phase.water.Mobility.Type        Constant
pfset Phase.water.Mobility.Value       1.0
\end{verbatim}

Again, \parflow{} has more capabilities than we are using here in the Cape Cod example. 
For this example, we handle monitoring wells in a separate code as we assume they
do not remove a significant amount of water from the domain.  Note that since 
there are no well names listed here, \parflow{} assumes we have no wells. 
If we had pumping wells, we would have to include them here, because they 
would affect the head distribution throughout our domain.  See below for an example of how to 
include pumping wells in this script.

\begin{verbatim}
#----------------------------------------------------------------------------
# Wells
#----------------------------------------------------------------------------
pfset Wells.Names ""
\end{verbatim}

You can give certain periods of time names if you want to (ie. Pre-injection, post-injection, etc).  
Here, however we do not have multiple time intervals and are simulating in steady state, so time cycle 
keys are simple.  We have only one time cycle and it's constant for the duration of the simulation. 
We accomplish this by giving it a repeat value of -1, which repeats indefinitely.  The length of 
the cycle is the length specified below (an integer) multiplied by the base unit value we specified earlier.

\begin{verbatim}
#----------------------------------------------------------------------------
# Time Cycles
#----------------------------------------------------------------------------
pfset Cycle.Names constant
pfset Cycle.constant.Names		"alltime"
pfset Cycle.constant.alltime.Length	 1
pfset Cycle.constant.Repeat		-1
\end{verbatim}

Now, we assign Boundary Conditions for each face (each of the Patches in the domain defined before).  
Recall the previously stated Patches and associate them with the boundary conditions that follow.

\begin{verbatim}
pfset BCPressure.PatchNames "left right front back bottom top"
\end{verbatim}

These are Dirichlet BCs (i.e. constant head over cell so the pressure head 
is set to hydrostatic-- \emph{see} \S~\ref{Boundary Conditions: Pressure}).  
There is no time dependence, so use the \code{constant} time cycle we defined previously.  
\code{RefGeom} links this to the established domain geometry and tells \parflow{} 
what to use for a datum when calculating hydrostatic head conditions.

\begin{verbatim}
pfset Patch.left.BCPressure.Type          DirEquilRefPatch
pfset Patch.left.BCPressure.Cycle        "constant"
pfset Patch.left.BCPressure.RefGeom domain
\end{verbatim}

Reference the current (left) patch to the bottom to define the line of intersection between the two.

\begin{verbatim}
pfset Patch.left.BCPressure.RefPatch  bottom
\end{verbatim}

Set the head permanently to 10.0m.  Pressure-head will of course vary top to bottom because of hydrostatics, 
but head potential will be constant.

\begin{verbatim}
pfset Patch.left.BCPressure.alltime.Value  10.0
\end{verbatim}

Repeat the declarations for the rest of the faces of the domain.  The left to right (\emph{X}) 
dimension is aligned with the hydraulic gradient.  The difference between the values assigned to 
right and left divided by the length of the domain corresponds to the correct hydraulic gradient.

\begin{verbatim}
pfset Patch.right.BCPressure.Type               DirEquilRefPatch
pfset Patch.right.BCPressure.Cycle              "constant"
pfset Patch.right.BCPressure.RefGeom       domain
pfset Patch.right.BCPressure.RefPatch       bottom
pfset Patch.right.BCPressure.alltime.Value 9.97501

pfset Patch.front.BCPressure.Type                FluxConst
pfset Patch.front.BCPressure.Cycle              "constant"
pfset Patch.front.BCPressure.alltime.Value 0.0

pfset Patch.back.BCPressure.Type               FluxConst
pfset Patch.back.BCPressure.Cycle              "constant"
pfset Patch.back.BCPressure.alltime.Value 0.0

pfset Patch.bottom.BCPressure.Type               FluxConst
pfset Patch.bottom.BCPressure.Cycle              "constant"
pfset Patch.bottom.BCPressure.alltime.Value 0.0

pfset Patch.top.BCPressure.Type			FluxConst
pfset Patch.top.BCPressure.Cycle			"constant"
pfset Patch.top.BCPressure.alltime.Value		0.0
\end{verbatim}

Next we define topographic slopes and Mannings \emph{n} values. 
These are not used, since we do not solve for overland flow.  
However, the keys still need to appear in the input script.

\begin{verbatim}
#---------------------------------------------------------
# Topo slopes in x-direction
#---------------------------------------------------------
# topo slopes do not figure into the impes (fully sat) case but we still
# need keys for them

pfset TopoSlopesX.Type "Constant"
pfset TopoSlopesX.GeomNames ""

pfset TopoSlopesX.Geom.domain.Value 0.0

#---------------------------------------------------------
# Topo slopes in y-direction
#---------------------------------------------------------

pfset TopoSlopesY.Type "Constant"
pfset TopoSlopesY.GeomNames ""

pfset TopoSlopesY.Geom.domain.Value 0.0

# You may also indicate an elevation file used to derive the slopes.
# This is optional but can be useful when post-processing terrain-
# following grids:
pfset TopoSlopes.Elevation.FileName "elevation.pfb"

#---------------------------------------------------------
# Mannings coefficient 
#---------------------------------------------------------
# mannings roughnesses do not figure into the impes (fully sat) case but we still
# need a key for them

pfset Mannings.Type "Constant"
pfset Mannings.GeomNames ""
pfset Mannings.Geom.domain.Value 0.
\end{verbatim}

Phase sources allows you to add sources other than wells and boundaries, but 
we do not have any so this key is constant, 0.0 over entire domain.

\begin{verbatim}
#----------------------------------------------------------------------------
# Phase sources:
#----------------------------------------------------------------------------

pfset PhaseSources.water.Type                         Constant
pfset PhaseSources.water.GeomNames                    domain
pfset PhaseSources.water.Geom.domain.Value        0.0
\end{verbatim}

Next we define solver parameters for {\bf IMPES}.  Since this is the default solver, we do not need a solver key.

\begin{verbatim}
#---------------------------------------------------------
#  Solver Impes  
#---------------------------------------------------------
\end{verbatim}

We allow up to 50 iterations of the linear solver before it quits or converges.

\begin{verbatim}
pfset Solver.MaxIter 50
\end{verbatim}

The solution must be accurate to this level 

\begin{verbatim}
pfset Solver.AbsTol  1E-10
\end{verbatim}

We drop significant digits beyond E-15

\begin{verbatim}
pfset Solver.Drop   1E-15

#--------------------------------------------------------
# Run and Unload the ParFlow output files
#---------------------------------------------------------
\end{verbatim}

Here you set the number of realizations again using a local tcl variable.  
We have set only one run but by setting the \file{n_runs} variable to something 
else we can run more than one realization of hydraulic conductivity. 

\begin{verbatim}
# this script is setup to run 100 realizations, for testing we just run one
###set n_runs 100
set n_runs 1
\end{verbatim}

Here is where you tell \parflow{} where to put the output.  
In this case, it is a directory called flow.  
Then you cd (change directory) into that new directory.  
If you wanted to put an entire path rather than just a name, 
you would have more control over where your output file goes.  
For example, you would put \file{file mkdir ``/cape_cod/revised_statistics/flow"} 
and then change into that directory. 

\begin{verbatim}
file mkdir "flow"
cd "flow"
\end{verbatim}

Now we loop through the realizations, again using tcl.  \file{k} is the
integer counter that is incremented for each realization.  When you use a variable
(rather than define it), you precede it with\verb+$+.  The hanging character \verb+{+ 
opens the do loop for \file{k}.

\begin{verbatim}
#
#  Loop through runs
#
for {set k 1} {$k <= $n_runs} {incr k 1} {
\end{verbatim}

The following expressions sets the variable \code{seed} equal to the expression in brackets,
which increments with each turn of the do loop and each seed will produce a different random field of {\bf K}.  
You set upper and lower aquifer, because in the Cape Cod site, these are the two subsets of the domain. 
Note the seed starts at a different point to allow for different random field generation for the upper 
and lower zones.

\begin{verbatim}
#
# set the random seed to be different for every run
#
pfset Geom.upper_aquifer.Perm.Seed  [ expr 33333+2*$k ] 
pfset Geom.lower_aquifer.Perm.Seed  [ expr 31313+2*$k ]
\end{verbatim}

The following command runs \parflow{} and gives you a suite of output files for each realization.  
The file names will begin \file{harvey_flow.1.xxxxx}, \file{harvey_flow.2.xxxx}, etc up to as many realizations as you run.  
The .xxxxx part will designate x, y, and z permeability, etc. Recall that in this case, since we normalized gravity, viscosity, 
and density, remember that we are really getting hydraulic conductivity.

\begin{verbatim}
pfrun harvey_flow.$k
\end{verbatim}

This command removes a large number of superfluous dummy files or un-distributes parallel files back into a single file. If you
compile with the \code{-- with-amps-sequential-io} option then a single \parflow{} file is written with corresponding \file{XXXX.dist} files and the 
\code{pfundist} command just removes these \file{.dist} files (though you don't really need to remove them if you don't want to).

\begin{verbatim}
pfundist harvey_flow.$k
\end{verbatim}

The following commands take advantage of PFTools (\emph{see} \S~\ref{PFTCL Commands}) and load pressure head output of the
/parflow{} model into a pressure matrix.

\begin{verbatim}
# we use pf tools to convert from pressure to head
# we could do a number of other things here like copy files to different
# format
set press [pfload harvey_flow.$k.out.press.pfb]
\end{verbatim}

The next command takes the pressures that were just loaded and converts it to head and loads them into a head matrix tcl variable.

\begin{verbatim}
set head [pfhhead $press]
\end{verbatim}

Finally, the head matrix is saved as a \parflow{} binary file (.pfb) and the k do loop is closed by the \verb+}+ character.  
Then we move up to the root directory when we are finished

\begin{verbatim}
 pfsave $head -pfb harvey_flow.$k.head.pfb
}

cd ".."
\end{verbatim}

Once you have modified the tcl input script (if necessary) and run \parflow{}, you will have as
many realizations of your subsurface as you specified.  Each of these realizations will be used 
as input for a particle or streamline calculation in the future.  We can see below, that since 
we have a tcl script as input, we can do a lot of different operations, for example, we might 
run a particle tracking transport code simulation using the results of the \parflow{} runs.  
This actually corresponds to the example presented in the \code{SLIM} user's manual.

\begin{verbatim}
# this could run other tcl scripts now an example is below
#puts stdout "running SLIM"
#source bromide_trans.sm.tcl
\end{verbatim}

We can add options to this script.  For example if we wanted to add a pumping well these additions
are described below.
%==============================
\subsection*{Adding a Pumping Well}

Let us change the input problem by adding a pumping well:
\begin{enumerate}


\item 
Add the following lines to the input file near where the existing
well information is in the input file.  You need to replace
the ``Wells.Names'' line with the one included here to get both
wells activated (this value lists the names of the wells):

\begin{display}\begin{verbatim}
pfset Wells.Names {new_well}

pfset Wells.new_well.InputType                Recirc

pfset Wells.new_well.Cycle		    constant

pfset Wells.new_well.ExtractionType	    Flux
pfset Wells.new_well.InjectionType            Flux

pfset Wells.new_well.X			    10.0 
pfset Wells.new_well.Y			    10.0
pfset Wells.new_well.ExtractionZLower	     0.5
pfset Wells.new_well.ExtractionZUpper	     0.5
pfset Wells.new_well.InjectionZLower	     0.2
pfset Wells.new_well.InjectionZUpper	     0.2

pfset Wells.new_well.ExtractionMethod	    Standard
pfset Wells.new_well.InjectionMethod          Standard

pfset Wells.new_well.alltime.Extraction.Flux.water.Value        	     0.50
pfset Wells.new_well.alltime.Injection.Flux.water.Value		     0.75

\end{verbatim}\end{display}

\end{enumerate}
For more information on defining the problem,
see \S~\ref{Defining the Problem}.

We could also visualize the results of the \parflow{} simulations, using \emph{VisIt}. 
For example, we can turn on \emph{SILO} file output which allows these files to be directly
read and visualized.  We would do this by adding the following \code{pfset} commands, I usually 
add them to the solver section:
\begin{display}\begin{verbatim}
pfset Solver.WriteSiloSubsurfData True
pfset Solver.WriteSiloPressure True
pfset Solver.WriteSiloSaturation True
\end{verbatim}\end{display}


You can then directly open the file \file{harvey_flow.#.out.perm_x.silo} (where \file{#} is the realization number).
The resulting image will be the hydraulic conductivity field of your domain, showing the variation in x-permeability in 3-D space.
You can also generate representations of head or pressure (or y or z permeability) throughout your domain using \parflow{} output files. 
See the section on visualization for more details.

%----------------------------------
\subsection{Little Washita Example}
\label{Little Washita Example}
This tutorial matches the \file{LW_Test.tcl} file found in the \file{/test/washita/tcl_scripts} 
directory and corresponds to \cite{Condon14a, Condon14b}.  This script runs the Little Washita domain for three days using \parflow{} \code{CLM} with 3D forcings. 
The domain is setup using terrain following grid (\S~\ref{TFG}) and subsurface geologes are specified using a \file{.pfb} indicator file. 
Input files were generated using the workflow detailed in \S~\ref{Defining a Real domain}. 

Now for the tcl script:

\begin{verbatim}
#
# Import the ParFlow TCL package
#
\end{verbatim}

These first three lines are what link \parflow{} and the tcl script, thus allowing you to use a set of commands seen later, such as \code{pfset}, etc.

\begin{verbatim}
lappend auto_path $env(PARFLOW_DIR)/bin 
package require parflow
namespace import Parflow::*

#-----------------------------------------------------------------------------
# File input version number
#-----------------------------------------------------------------------------
pfset FileVersion 4
\end{verbatim}

These next lines set the parallel process topology.  The domain is divided in \emph{x},\emph{y} and \emph{z} by \code{P}, \code{Q} and \code{R}.
The total number of processors is \code{P*Q*R}  (\emph{see} \S~\ref{Computing Topology}).

\begin{verbatim}
#----------------------------------------------------------------------------
# Process Topology
#----------------------------------------------------------------------------

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

Before we really get started make a directory for our outputs and copy all of the required input files into the run directory. 
These files will be described in detail later as they get used.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Make a directory for the simulation and copy inputs into it
#-----------------------------------------------------------------------------
exec mkdir "Outputs"
cd "./Outputs"

# ParFlow Inputs
file copy -force "../../parflow_input/LW.slopex.pfb" .
file copy -force "../../parflow_input/LW.slopey.pfb" .
file copy -force "../../parflow_input/IndicatorFile_Gleeson.50z.pfb"   .
file copy -force "../../parflow_input/press.init.pfb"  .

#CLM Inputs
file copy -force "../../clm_input/drv_clmin.dat" .
file copy -force "../../clm_input/drv_vegp.dat"  .
file copy -force "../../clm_input/drv_vegm.alluv.dat"  . 

puts "Files Copied"
\end{verbatim}

Next we set up the computational grid (\emph{see} \S~\ref{Defining the Problem} and  \S~\ref{Computational Grid}).

\begin{verbatim}
#----------------------------------------------------------------------------
# Computational Grid
#----------------------------------------------------------------------------
\end{verbatim}

Locate the origin in the domain.

\begin{verbatim}
pfset ComputationalGrid.Lower.X    0.0
pfset ComputationalGrid.Lower.Y    0.0
pfset ComputationalGrid.Lower.Z    0.0
\end{verbatim}

Define the size of the domain grid block.  Length units, same as those on hydraulic conductivity.

\begin{verbatim}
pfset ComputationalGrid.DX     1000.0
pfset ComputationalGrid.DY     1000.0
pfset ComputationalGrid.DZ     2.0
\end{verbatim}

Define the number of grid blocks in the domain.

\begin{verbatim}
pfset ComputationalGrid.NX      41
pfset ComputationalGrid.NY      41
pfset ComputationalGrid.NZ      50
\end{verbatim}

This next piece is comparable to a pre-declaration of variables.  These will be areas in our domain geometry. 
The regions themselves will be defined later.  You must always have one that is the name of your entire domain.  
If you want subsections within your domain, you may declare these as well. Here we define two geometries one is the
domain and one is for the indicator file (which will also span the entire domain).

\begin{verbatim}
#-----------------------------------------------------------------------------
# The Names of the GeomInputs
#-----------------------------------------------------------------------------
pfset GeomInput.Names                     "box_input indi_input"
\end{verbatim}

Now you characterize the domain that you just pre-declared to be a \code{box} (\emph{see} \S~\ref{Geometries}),
and you also give it a name, \code{domain}.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Domain Geometry Input
#-----------------------------------------------------------------------------
pfset GeomInput.box_input.InputType      Box
pfset GeomInput.box_input.GeomName      domain
\end{verbatim}

Here, you set the limits in space for your entire domain.  The span from \file{Lower.X} to \file{Upper.X} 
will be equal to the product of \file{ComputationalGrid.DX} times \file{ComputationalGrid.NX}.  
Same for {\em Y} and {\em Z} (i.e. the number of grid elements times size of the grid element has to
equal the size of the grid in each dimension).   The \code{Patches} key assigns names to the outside edges, 
because the domain is the limit of the problem in space.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Domain Geometry 
#-----------------------------------------------------------------------------
pfset Geom.domain.Lower.X                        0.0
pfset Geom.domain.Lower.Y                        0.0
pfset Geom.domain.Lower.Z                        0.0
 
pfset Geom.domain.Upper.X                        41000.0
pfset Geom.domain.Upper.Y                        41000.0
pfset Geom.domain.Upper.Z                          100.0

pfset Geom.domain.Patches             "x-lower x-upper y-lower y-upper z-lower z-upper"
\end{verbatim}

Now we setup the indicator file. As noted above, the indicator  file has integer values for every grid
cell in the domain designating what geologic unit it belongs to. The \file{GeomNames} list
should include a name for every unit in your indicator file. In this example we have thirteen soil
units and eight geologic units. The \file{FileName} points to the indicator file that \parflow{}
will read. Recall that this file into the run directory at the start of the script. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Indicator Geometry Input
#-----------------------------------------------------------------------------
pfset GeomInput.indi_input.InputType      IndicatorField
pfset GeomInput.indi_input.GeomNames      "s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 g1 g2 g3 g4 g5 g6 g7 g8"
pfset Geom.indi_input.FileName            "IndicatorFile_Gleeson.50z.pfb"
\end{verbatim}

For every name in the \file{GeomNames} list we define the corresponding value in the indicator file. 
For example, here we are saying that our first soil unit (\file{s1}) is represented by the number ``1"
in the indicator file, while the first geologic unit (\file{g1}) is represented by the number ``21".
Note that the integers used in the indicator file do not need to be consecutive.

\begin{verbatim}
pfset GeomInput.s1.Value                1
pfset GeomInput.s2.Value                2
pfset GeomInput.s3.Value                3
pfset GeomInput.s4.Value                4
pfset GeomInput.s5.Value                5
pfset GeomInput.s6.Value                6
pfset GeomInput.s7.Value                7
pfset GeomInput.s8.Value                8
pfset GeomInput.s9.Value                9
pfset GeomInput.s10.Value               10
pfset GeomInput.s11.Value               11
pfset GeomInput.s12.Value               12
pfset GeomInput.s13.Value               13
pfset GeomInput.g1.Value                21
pfset GeomInput.g2.Value                22
pfset GeomInput.g3.Value                23
pfset GeomInput.g4.Value                24
pfset GeomInput.g5.Value                25
pfset GeomInput.g6.Value                26
pfset GeomInput.g7.Value                27
pfset GeomInput.g8.Value                28
\end{verbatim}

Now you add permeability data to the domain sections defined above (\S~\ref{Permeability}).  
You can reassign values simply by re-stating them -- there is no need to comment out or 
delete the previous version -- the final statement is the only one that counts. Also, note
that you do not need to assign permeability values to all of the geometries names. 
Any geometry that is not assigned its own permeability value will take the \file{domain}
value. However, every geometry listed in \file{Porosity.GeomNames} must have 
values assigned. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Permeability (values in m/hr)
#-----------------------------------------------------------------------------
pfset Geom.Perm.Names                     "domain s1 s2 s3 s4 s5 s6 s7 s8 s9 g2 g3 g6 g8"

pfset Geom.domain.Perm.Type           Constant
pfset Geom.domain.Perm.Value          0.2

pfset Geom.s1.Perm.Type               Constant
pfset Geom.s1.Perm.Value              0.269022595

pfset Geom.s2.Perm.Type               Constant
pfset Geom.s2.Perm.Value              0.043630356

pfset Geom.s3.Perm.Type               Constant
pfset Geom.s3.Perm.Value              0.015841225

pfset Geom.s4.Perm.Type               Constant
pfset Geom.s4.Perm.Value              0.007582087

pfset Geom.s5.Perm.Type               Constant
pfset Geom.s5.Perm.Value              0.01818816

pfset Geom.s6.Perm.Type               Constant
pfset Geom.s6.Perm.Value              0.005009435

pfset Geom.s7.Perm.Type               Constant
pfset Geom.s7.Perm.Value              0.005492736

pfset Geom.s8.Perm.Type               Constant
pfset Geom.s8.Perm.Value              0.004675077

pfset Geom.s9.Perm.Type               Constant
pfset Geom.s9.Perm.Value              0.003386794

pfset Geom.g2.Perm.Type               Constant
pfset Geom.g2.Perm.Value              0.025

pfset Geom.g3.Perm.Type               Constant
pfset Geom.g3.Perm.Value              0.059

pfset Geom.g6.Perm.Type               Constant
pfset Geom.g6.Perm.Value              0.2

pfset Geom.g8.Perm.Type              Constant
pfset Geom.g8.Perm.Value             0.68
\end{verbatim}

The following section allows you to specify the permeability tensor.  
In the case below, permeability is symmetric in all directions (x, y, and z) and therefore each is set to 1.0.
Also note that we just specify this once for the whole domain because we want isotropic
permeability everywhere. You can specify different tensors for different units by repeating 
these lines with different \file{Geom.Names}.

\begin{verbatim}
pfset Perm.TensorType                     TensorByGeom
pfset Geom.Perm.TensorByGeom.Names        "domain"
pfset Geom.domain.Perm.TensorValX         1.0d0
pfset Geom.domain.Perm.TensorValY         1.0d0
pfset Geom.domain.Perm.TensorValZ         1.0d0
\end{verbatim}

Next we set the specific storage. Here again we specify one value for the whole domain but
these lines can be easily repeated to set different values for different units. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Specific Storage
#-----------------------------------------------------------------------------
pfset SpecificStorage.Type                Constant
pfset SpecificStorage.GeomNames           "domain"
pfset Geom.domain.SpecificStorage.Value   1.0e-5
\end{verbatim}

\parflow{} has the capability to deal with a multiphase system, but we only have one (water) in this example. 
As we stated earlier, we set density and viscosity artificially (and later gravity) both to 1.0.  
Again, this is merely a trick to solve for hydraulic conductivity and pressure head.  
If you were to set density and viscosity to their true values, the code would calculate {\bf k} (permeability). 
By using the \emph{normalized} values instead, you effectively embed the conversion of {\bf k} to {\bf K} (hydraulic conductivity). 
So this way, we get hydraulic conductivity, which is what we want for this problem.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Phases
#-----------------------------------------------------------------------------
pfset Phase.Names                         "water"

pfset Phase.water.Density.Type            Constant
pfset Phase.water.Density.Value           1.0

pfset Phase.water.Viscosity.Type          Constant
pfset Phase.water.Viscosity.Value         1.0
\end{verbatim}

This example does not include the \parflow{} grid based transport scheme.  Therefore we leave 
contaminants blank.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Contaminants
#-----------------------------------------------------------------------------
pfset Contaminants.Names                  ""
\end{verbatim}

As with density and viscosity, gravity is normalized here.  If we used the true value
(in the \emph{[L]} and \emph{[T]} units of hydraulic conductivity) the code would be calculating permeability.  
Instead, we normalize so that the code calculates hydraulic conductivity.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Gravity
#-----------------------------------------------------------------------------
pfset Gravity                             1.0
\end{verbatim}

Next we set up the timing for our simulation.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Timing (time units is set by units of permeability)
#-----------------------------------------------------------------------------
\end{verbatim}

This specifies the base unit of time for all time values entered. All time should be expressed
as multiples of this value. To keep things simple here we set it to 1. Because we expressed
our permeability in units of m/hr in this example this means that our basin unit of time is 1hr. 

\begin{verbatim}
pfset TimingInfo.BaseUnit                 1.0
\end{verbatim}

This key specifies the time step number that will be associated with the first advection cycle
of the transient problem. Because we are starting from scratch we set this to 0. If we were 
restarting a run we would set this to the last time step of your previous simulation. Refer to 
\S~\ref{Restarting a Run} for additional instructions on restarting a run.

\begin{verbatim}
pfset TimingInfo.StartCount               0.0
\end{verbatim}

\file{StartTime} and \file{StopTime} specify the start and stop times for the simulation. These 
values should correspond with the forcing files you are using. 

\begin{verbatim}
pfset TimingInfo.StartTime                0.0
pfset TimingInfo.StopTime                 72.0
\end{verbatim}

This key specifies the timing interval at which \parflow{} time dependent outputs will be written. 
Here we have a base unit of 1hr so a dump interval of 24 means that we are writing daily outputs. 
Note that this key only controls the \parflow{} output interval and not the interval that \code{CLM}
outputs will be written out at. 

\begin{verbatim}
pfset TimingInfo.DumpInterval             24.0
\end{verbatim}

Here we set the time step value. For this example we use a constant time step of 1hr. 

\begin{verbatim}
pfset TimeStep.Type                       Constant
pfset TimeStep.Value                      1.0
\end{verbatim}

Next, we assign the porosity (\emph{see} \S~\ref{Porosity}).  As with the permeability
we assign different values for different indicator geometries. Here we assign 
values for all of our soil units but not for the geologic units, they will default to the
domain value of 0.4. Note that every geometry listed in \file{Porosity.GeomNames} must have 
values assigned. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Porosity
#-----------------------------------------------------------------------------
pfset Geom.Porosity.GeomNames             "domain s1 s2 s3 s4 s5 s6 s7 s8 s9"

pfset Geom.domain.Porosity.Type          Constant
pfset Geom.domain.Porosity.Value         0.4

pfset Geom.s1.Porosity.Type    Constant
pfset Geom.s1.Porosity.Value   0.375

pfset Geom.s2.Porosity.Type    Constant
pfset Geom.s2.Porosity.Value   0.39

pfset Geom.s3.Porosity.Type    Constant
pfset Geom.s3.Porosity.Value   0.387

pfset Geom.s4.Porosity.Type    Constant
pfset Geom.s4.Porosity.Value   0.439

pfset Geom.s5.Porosity.Type    Constant
pfset Geom.s5.Porosity.Value   0.489

pfset Geom.s6.Porosity.Type    Constant
pfset Geom.s6.Porosity.Value   0.399

pfset Geom.s7.Porosity.Type    Constant
pfset Geom.s7.Porosity.Value   0.384

pfset Geom.s8.Porosity.Type            Constant
pfset Geom.s8.Porosity.Value           0.482

pfset Geom.s9.Porosity.Type            Constant
pfset Geom.s9.Porosity.Value           0.442
\end{verbatim}

Having defined the geometry of our problem before and named it \code{domain}, we are
now ready to report/upload that problem, which we do here.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Domain
#-----------------------------------------------------------------------------
pfset Domain.GeomName                     "domain"
\end{verbatim}

Mobility between phases is set to 1.0 because we only have one phase (water).

\begin{verbatim}
#----------------------------------------------------------------------------
# Mobility
#----------------------------------------------------------------------------
pfset Phase.water.Mobility.Type        Constant
pfset Phase.water.Mobility.Value       1.0
\end{verbatim}

Again, \parflow{} has more capabilities than we are using here in this example. 
Note that since there are no well names listed here, \parflow{} assumes we have no wells. 
If we had pumping wells, we would have to include them here, because they 
would affect the head distribution throughout our domain.  See \S~\ref{Harvey Flow Example} 
for an example of how to include pumping wells in this script.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Wells
#-----------------------------------------------------------------------------
pfset Wells.Names                         ""
\end{verbatim}

You can give certain periods of time names if you want. For example if you aren't running with \code{CLM} 
and you would like to have periods with rain and periods without.   
Here, however we have only one time cycle because \code{CLM} will handle the variable forcings.
Therefore, we specify one time cycle and it's constant for the duration of the simulation. 
We accomplish this by giving it a repeat value of -1, which repeats indefinitely.  The length of 
the cycle is the length specified below (an integer) multiplied by the base unit value we specified earlier.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Time Cycles
#-----------------------------------------------------------------------------
pfset Cycle.Names                         "constant"
pfset Cycle.constant.Names                "alltime"
pfset Cycle.constant.alltime.Length        1
pfset Cycle.constant.Repeat               -1
\end{verbatim}

Now, we assign Boundary Conditions for each face (each of the Patches in the domain defined before).  
Recall the previously stated Patches and associate them with the boundary conditions that follow.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Boundary Conditions
#-----------------------------------------------------------------------------
pfset BCPressure.PatchNames                   [pfget Geom.domain.Patches]
\end{verbatim}

The bottom and sides of our domain are all set to no-flow (i.e. constant flux of 0) boundaries.

\begin{verbatim}
pfset Patch.x-lower.BCPressure.Type		      FluxConst
pfset Patch.x-lower.BCPressure.Cycle		      "constant"
pfset Patch.x-lower.BCPressure.alltime.Value	      0.0

pfset Patch.y-lower.BCPressure.Type		      FluxConst
pfset Patch.y-lower.BCPressure.Cycle		      "constant"
pfset Patch.y-lower.BCPressure.alltime.Value	      0.0

pfset Patch.z-lower.BCPressure.Type		      FluxConst
pfset Patch.z-lower.BCPressure.Cycle		      "constant"
pfset Patch.z-lower.BCPressure.alltime.Value	      0.0

pfset Patch.x-upper.BCPressure.Type		      FluxConst
pfset Patch.x-upper.BCPressure.Cycle		      "constant"
pfset Patch.x-upper.BCPressure.alltime.Value	      0.0

pfset Patch.y-upper.BCPressure.Type		      FluxConst
pfset Patch.y-upper.BCPressure.Cycle		      "constant"
pfset Patch.y-upper.BCPressure.alltime.Value	      0.0
\end{verbatim}

The top is set to an \code{OverlandFLow} boundary to turn on the fully-coupled overland flow
routing. 

\begin{verbatim}
pfset Patch.z-upper.BCPressure.Type		      OverlandFlow
pfset Patch.z-upper.BCPressure.Cycle		      "constant"
pfset Patch.z-upper.BCPressure.alltime.Value	      0.0
\end{verbatim}

Next we define topographic slopes and values. These slope values were derived from a digital
elevation model of the domain following the workflow outlined in \S~\ref{Defining a Real domain}.
In this example we read the slope files in from \file{.pfb} files that were copied into the run directory
at the start of this script. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Topo slopes in x-direction
#-----------------------------------------------------------------------------
pfset TopoSlopesX.Type                                "PFBFile"
pfset TopoSlopesX.GeomNames                           "domain"
pfset TopoSlopesX.FileName                            "LW.slopex.pfb"

#-----------------------------------------------------------------------------
# Topo slopes in y-direction
#-----------------------------------------------------------------------------
pfset TopoSlopesY.Type                                "PFBFile"
pfset TopoSlopesY.GeomNames                           "domain"
pfset TopoSlopesY.FileName                            "LW.slopey.pfb"
\end{verbatim}

And now we define the Mannings \emph{n}, again just one value for the whole domain in this example. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Mannings coefficient
#-----------------------------------------------------------------------------
pfset Mannings.Type                                   "Constant"
pfset Mannings.GeomNames                              "domain"
pfset Mannings.Geom.domain.Value                      5.52e-6
\end{verbatim}

Following the same approach as we did for \file{Porosity} we define the relative permeability
inputs that will be used for Richards' equation implementation (\S~\ref{Richards RelPerm}).
Here we use \file{VanGenuchten} parameters. Note that every geometry listed in \file{Porosity.GeomNames} must have 
values assigned. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Relative Permeability
#-----------------------------------------------------------------------------
pfset Phase.RelPerm.Type                  VanGenuchten
pfset Phase.RelPerm.GeomNames             "domain s1 s2 s3 s4 s5 s6 s7 s8 s9 "

pfset Geom.domain.RelPerm.Alpha           3.5
pfset Geom.domain.RelPerm.N               2.0

pfset Geom.s1.RelPerm.Alpha        3.548
pfset Geom.s1.RelPerm.N            4.162

pfset Geom.s2.RelPerm.Alpha        3.467
pfset Geom.s2.RelPerm.N            2.738

pfset Geom.s3.RelPerm.Alpha        2.692
pfset Geom.s3.RelPerm.N            2.445

pfset Geom.s4.RelPerm.Alpha        0.501
pfset Geom.s4.RelPerm.N            2.659

pfset Geom.s5.RelPerm.Alpha        0.661
pfset Geom.s5.RelPerm.N            2.659

pfset Geom.s6.RelPerm.Alpha        1.122
pfset Geom.s6.RelPerm.N            2.479

pfset Geom.s7.RelPerm.Alpha        2.089
pfset Geom.s7.RelPerm.N            2.318

pfset Geom.s8.RelPerm.Alpha        0.832
pfset Geom.s8.RelPerm.N            2.514

pfset Geom.s9.RelPerm.Alpha        1.585
pfset Geom.s9.RelPerm.N            2.413
\end{verbatim}

Next we do the same thing for saturation (\S~\ref{Saturation}) again using the \file{VanGenuchten} parameters
Note that every geometry listed in \file{Porosity.GeomNames} must have 
values assigned. 

\begin{verbatim}
#-----------------------------------------------------------------------------
# Saturation
#-----------------------------------------------------------------------------
pfset Phase.Saturation.Type               VanGenuchten
pfset Phase.Saturation.GeomNames          "domain s1 s2 s3 s4 s5 s6 s7 s8 s9 "

pfset Geom.domain.Saturation.Alpha        3.5
pfset Geom.domain.Saturation.N            2.
pfset Geom.domain.Saturation.SRes         0.2
pfset Geom.domain.Saturation.SSat         1.0

pfset Geom.s1.Saturation.Alpha        3.548
pfset Geom.s1.Saturation.N            4.162
pfset Geom.s1.Saturation.SRes         0.000001
pfset Geom.s1.Saturation.SSat         1.0

pfset Geom.s2.Saturation.Alpha        3.467
pfset Geom.s2.Saturation.N            2.738
pfset Geom.s2.Saturation.SRes         0.000001
pfset Geom.s2.Saturation.SSat         1.0

pfset Geom.s3.Saturation.Alpha        2.692
pfset Geom.s3.Saturation.N            2.445
pfset Geom.s3.Saturation.SRes         0.000001
pfset Geom.s3.Saturation.SSat         1.0

pfset Geom.s4.Saturation.Alpha        0.501
pfset Geom.s4.Saturation.N            2.659
pfset Geom.s4.Saturation.SRes         0.000001
pfset Geom.s4.Saturation.SSat         1.0

pfset Geom.s5.Saturation.Alpha        0.661
pfset Geom.s5.Saturation.N            2.659
pfset Geom.s5.Saturation.SRes         0.000001
pfset Geom.s5.Saturation.SSat         1.0

pfset Geom.s6.Saturation.Alpha        1.122
pfset Geom.s6.Saturation.N            2.479
pfset Geom.s6.Saturation.SRes         0.000001
pfset Geom.s6.Saturation.SSat         1.0

pfset Geom.s7.Saturation.Alpha        2.089
pfset Geom.s7.Saturation.N            2.318
pfset Geom.s7.Saturation.SRes         0.000001
pfset Geom.s7.Saturation.SSat         1.0

pfset Geom.s8.Saturation.Alpha        0.832
pfset Geom.s8.Saturation.N            2.514
pfset Geom.s8.Saturation.SRes         0.000001
pfset Geom.s8.Saturation.SSat         1.0

pfset Geom.s9.Saturation.Alpha        1.585
pfset Geom.s9.Saturation.N            2.413
pfset Geom.s9.Saturation.SRes         0.000001
pfset Geom.s9.Saturation.SSat         1.0
\end{verbatim}

Phase sources allows you to add sources other than wells and boundaries, but 
we do not have any so this key is constant, 0.0 over entire domain.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Phase sources:
#-----------------------------------------------------------------------------
pfset PhaseSources.water.Type                         "Constant"
pfset PhaseSources.water.GeomNames                    "domain"
pfset PhaseSources.water.Geom.domain.Value            0.0
\end{verbatim}

In this example we are using \parflow{} \file{CLM} so we must provide some parameters for \file{CLM} (\S~\ref{CLM Solver Parameters}).
Note that \file{CLM} will also require some additional inputs outside of the tcl script. Refer to 
\file{/washita/clm_input/} for examples of the \code{CLM} \file{vegm} and \file{driver} files. These inputs
are also discussed briefly in \S~\ref{Defining a Real domain}.

\begin{verbatim}
#----------------------------------------------------------------
# CLM Settings:
# ------------------------------------------------------------
\end{verbatim}

First we specify that we will be using \code{CLM} as the land surface model and provide the name of a directory
that outputs will be written to. For this example we do not need outputs for each processor or a binary output directory. 
Finally we set the dump interval to 1, indicating that we will be writing outputs for every time step. Note that
this does not have to match the dump interval for \parflow{} outputs. Recall that earlier we set the \parflow{} dump
interval to 24. 

\begin{verbatim}
pfset Solver.LSM                                      CLM
pfset Solver.CLM.CLMFileDir                           "clm_output/"
pfset Solver.CLM.Print1dOut                           False
pfset Solver.BinaryOutDir                             False
pfset Solver.CLM.CLMDumpInterval                      1
\end{verbatim}

Next we specify the details of the meteorological forcing files that \code{clm} will read. First we provide the 
name of the files and the directory they can be found in. Next we specify that we are using
\file{3D} forcing files meaning that we have spatially distributed forcing with multiple time steps in every 
file. Therefore we must also specify the number of times steps (\code{MetFileNT}) in every file, in this case 24.  
Finally, we specify the initial value for the CLM counter. 

\begin{verbatim}
pfset Solver.CLM.MetFileName                          "NLDAS"
pfset Solver.CLM.MetFilePath                          "../../NLDAS/"
pfset Solver.CLM.MetForcing                           3D
pfset Solver.CLM.MetFileNT                            24
pfset Solver.CLM.IstepStart                           1
\end{verbatim}

This last set of \code{CLM} parameters refers to the physical properties of the system. Refer to \S~\ref{CLM Solver Parameters}
for details. 
\begin{verbatim}
pfset Solver.CLM.EvapBeta                             Linear
pfset Solver.CLM.VegWaterStress                       Saturation
pfset Solver.CLM.ResSat                               0.1
pfset Solver.CLM.WiltingPoint                         0.12
pfset Solver.CLM.FieldCapacity                        0.98
pfset Solver.CLM.IrrigationType                       none
\end{verbatim}

Next we set the initial conditions for the domain. In this example we are using a pressure \file{.pfb} 
file that was obtained by spinning up the model in the workflow outlined in \S~\ref{Defining a Real domain}.
Alternatively, the water table can be set to a constant value by changing the \code{ICPressure.Type}. Again, 
the input file that is referenced here was was copied into the run directory at the top of this script. 

\begin{verbatim}
#---------------------------------------------------------
# Initial conditions: water pressure
#---------------------------------------------------------
pfset ICPressure.Type                                 PFBFile
pfset ICPressure.GeomNames                            domain
pfset Geom.domain.ICPressure.RefPatch                   z-upper
pfset Geom.domain.ICPressure.FileName                   press.init.pfb
\end{verbatim}

Now we specify what outputs we would like written. In this example we specify that we would like to write
out \code{CLM} variables as well as \file{Pressure} and \file{Saturation}. However, there are many options 
for this and you should change these options according to what type of analysis you will be performing on your 
results. A complete list of print options is provided in \S~\ref{Code Parameters}. 

\begin{verbatim}
#----------------------------------------------------------------
# Outputs
# ------------------------------------------------------------
#Writing output (all pfb):
pfset Solver.PrintSubsurfData                         False
pfset Solver.PrintPressure                            True
pfset Solver.PrintSaturation                          True
pfset Solver.PrintMask                                True

pfset Solver.WriteCLMBinary                           False
pfset Solver.PrintCLM                                 True
pfset Solver.WriteSiloSpecificStorage                 False
pfset Solver.WriteSiloMannings                        False
pfset Solver.WriteSiloMask                            False
pfset Solver.WriteSiloSlopes                          False
pfset Solver.WriteSiloSubsurfData                     False
pfset Solver.WriteSiloPressure                        False
pfset Solver.WriteSiloSaturation                      False
pfset Solver.WriteSiloEvapTrans                       False
pfset Solver.WriteSiloEvapTransSum                    False
pfset Solver.WriteSiloOverlandSum                     False
pfset Solver.WriteSiloCLM                             False
\end{verbatim}

Next we specify the solver settings for the ParFlow (\S~\ref{RE Solver Parameters}).
First we turn on solver Richards and the terrain following grid.  We turn off 
variable dz.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Set solver parameters
#-----------------------------------------------------------------------------
# ParFlow Solution
pfset Solver                                          Richards
pfset Solver.TerrainFollowingGrid                     True
pfset Solver.Nonlinear.VariableDz                     False
\end{verbatim}

We then set the max solver settings and linear and nonlinear convergence
tolerance settings. The linear system will be solved to a norm of $10^{-8}$ and 
the nonlinear system will be solved to less than $10^{-6}$. Of note in latter key block 
is the EtaChoice and that we use the analytical Jacobian (\emph{UseJacobian=\bf{True}}).
We are using the \emph{FullJacobian} preconditioner, which is a more robust approach
but is more expensive.
\begin{verbatim}
pfset Solver.MaxIter                                  25000
pfset Solver.Drop                                     1E-20
pfset Solver.AbsTol                                   1E-8
pfset Solver.MaxConvergenceFailures                   8
pfset Solver.Nonlinear.MaxIter                        80
pfset Solver.Nonlinear.ResidualTol                    1e-6

pfset Solver.Nonlinear.EtaChoice                         EtaConstant
pfset Solver.Nonlinear.EtaValue                          0.001
pfset Solver.Nonlinear.UseJacobian                       True 
pfset Solver.Nonlinear.DerivativeEpsilon                 1e-16
pfset Solver.Nonlinear.StepTol				 			1e-30
pfset Solver.Nonlinear.Globalization                     LineSearch
pfset Solver.Linear.KrylovDimension                      70
pfset Solver.Linear.MaxRestarts                           2

pfset Solver.Linear.Preconditioner                       PFMG
pfset Solver.Linear.Preconditioner.PCMatrixType     FullJacobian
\end{verbatim}

This key is just for testing the Richards' formulation, so we are not using it.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Exact solution specification for error calculations
#-----------------------------------------------------------------------------
pfset KnownSolution                                   NoKnownSolution
\end{verbatim}

Next we distribute all the inputs as described by the keys in \S~\ref{PFTCL Commands}.  Note the slopes 
are 2D files, while the rest of the \parflow{} inputs are 3D so we need to alter the NZ accordingly
following example \ref{dist example} in \S~\ref{common_pftcl}.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Distribute inputs
#-----------------------------------------------------------------------------
pfset ComputationalGrid.NX                41 
pfset ComputationalGrid.NY                41 
pfset ComputationalGrid.NZ                1
pfdist LW.slopex.pfb
pfdist LW.slopey.pfb

pfset ComputationalGrid.NX                41 
pfset ComputationalGrid.NY                41 
pfset ComputationalGrid.NZ                50 
pfdist IndicatorFile_Gleeson.50z.pfb
pfdist press.init.pfb
\end{verbatim}

Now we run the simulation.  Note that we use a tcl variable to set the run name.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Run Simulation
#-----------------------------------------------------------------------------
set runname "LW"
puts $runname
pfrun    $runname
\end{verbatim}

All that is left is to undistribute files.

\begin{verbatim}
#-----------------------------------------------------------------------------
# Undistribute Files
#-----------------------------------------------------------------------------
pfundist $runname
pfundist press.init.pfb
pfundist LW.slopex.pfb
pfundist LW.slopey.pfb
pfundist IndicatorFile_Gleeson.50z.pfb

puts "ParFlow run Complete"
\end{verbatim}



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