7. Manipulating Data: PFTools¶
7.1. Introduction to the ParFlow TCL commands (PFTCL)¶
Several tools for manipulating data are provided in PFTCL command set. Tools can be accessed directly from the TCL shell or within a ParFlow input script. In both cases you must first load the ParFlow package into the TCL shell as follows:
#
# To Import the ParFlow TCL package
#
lappend auto_path $env(PARFLOW_DIR)/bin
package require parflow
namespace import Parflow::*
In addition to these methods xpftools provides GUI access to most of
these features. However the simplest approach is generally to include
the tools commands within a tcl script. The following section lists all
of the available ParFlow TCL commands along with detailed instructions
for their use. §8.2 PFTCL Commands provides several examples of
pre and post processing using the tools. In addition, a list of tools
can be obtained by typing pfhelp
into a TCL shell after importing
ParFlow. Typing ¿pfhelp¿
followed by a command name will display a
detailed description of the command in question.
7.2. PFTCL Commands¶
The tables that follow 4.1, 4.2 and 4.3 provide a list of ParFlow commands with short descriptions grouped according to their function. The last two columns in this table indicate what examples from §4.3 Common examples using ParFlow TCL commands (PFTCL), if any, the command is used in and whether the command is compatible with a terrain following grid domain formulation.
Name 
Short Description 
Examples 
Compatible with TFG? 

pfhelp 
Get help for PF Tools 
X 

Mathematical Operations 

pfcellsum 
datasetx + datasety 
X 

pfcelldiff 
datasetx  datasety 
X 

pfcellmult 
datasetx * datasety 
X 

pfcelldiv 
datasetx / datasety 
X 

pfcellsumconst 
dataset + constant 
X 

p fcelldiffconst 
dataset  constant 
X 

p fcellmultconst 
dataset * constant 
X 

pfcelldivconst 
dataset / constant 
X 

pfsum 
Sum dataset 
7, 9 
X 
pfdiffelt 
Element difference 
X 

pfprintdiff 
Print difference 
X 

pfmdiff 
Calculate area where the difference between two datasets is less than a threshold 
X 

pfprintmdiff 
Print the locations with differences greater than a minimum threshold 
X 

pfsavediff 
Save the difference between two datasets 
X 

pfaxpy 
y=alpha*x+y 
X 

pfgetstats 
Calculate dataset statistics (min, max, mean, var, stdev) 
X 

pfprintstats 
Print formatted statistics 
X 

pfstats 
Calculate and print dataset statistics (min, max, mean, var, stdev) 
X 

Calculate physical parameters 

pfbfcvel 
Calculate block face centered velocity 

pfcvel 
Calculate Darcy velocity 

pfvvel 
Calculate Darcy velocity at cell vertices 

pfvmag 
Calculate velocity magnitude given components 

pfflux 
Calculate Darcy flux 

pfhhead 
Calculate hydraulic head 
2 

pfphead 
Calculate pressure head from hydraulic head 

pfsattrans 
calculate saturated transmissivity 
X 

pfupstreamarea 
Calculate upstream area 
X 

pfeff ectiverecharge 
Calculate effective recharge 
X 

pfw atertabledepth 
Calculate water table from saturation 
X 

pfhydrostatic 
Calculate hydrostatic pressure field 

pfsub surfacestorage 
Calculate total subsurface storage 
7 
X 
pfgwstorage 
Calculate saturated subsurface storage 
X 

p fsurfacerunoff 
Calculate total surface runoff 
9 
X 
pf surfacestorage 
Calculate total surface storage 
8 
X 
Name 
Short Description 
Examples 
Compatible with TFG? 

DEM Operations 

pfslopex 
Calculate slopes in the xdirection 
5 
X 
pfslopey 
Calculate slope in the ydirection 
5 
X 
pfchildD8 
Calculate D8 child 
X 

pfsegmentD8 
Calculate D8 segment lengths 
X 

pfslopeD8 
Calculate D8 slopes 
X 

pfslopexD4 
Calculate D4 slopes in the xdirection 
X 

pfslopeyD4 
Calculate D4 slopes in the ydirection 
X 

pffillflats 
Fill DEM flats 
5 
X 
pfmovingavgdem 
Fill dem sinks with moving average 
X 

pfpitfilldem 
Fill sinks in the dem using iterative pitfilling routine 
5 
X 
pfflintslawfit 
Calculate Flint’s Law parameters 
X 

pfflintslaw 
Smooth DEM using Flints Law 
X 

pffl intslawbybasin 
Smooth DEM using Flints Law by basin 
X 

Topmodel functions 

pftopodeficit 
Calculate TOPMODEL water deficit 
X 

pftopoindex 
Calculate topographic index 
X 

pftopowt 
Calculate watertable based on topographic index 
X 

pftoporecharge 
Calculate effective recharge 
X 

Domain Operations 

p fcomputedomain 
Compute domain mask 
3 
X 
pfcomputetop 
Compute domain top 
3, 6, 8, 9 
X 
pfextracttop 
Extract domain top 
6 
X 
p fcomputebottom 
Compute domain bottom 
3 
X 
pfsetgrid 
Set grid 
5 
X 
pfgridtype 
Set grid type 
X 

pfgetgrid 
Return grid information 
X 

pfgetelt 
Extract element from domain 
10 
X 
pfe xtract2Ddomain 
Build 2D domain 
X 

pfenlargebox 
Compute expanded dataset 
X 

pfgetsubbox 
Return subset of data 
X 

pfprintdomain 
Print domain 
3 
X 
pfbuilddomain 
Build a subgrid array from a ParFlow database 
X 

Dataset operations 

pflistdata 
Return dataset names and labels 
X 

pfgetlist 
Return dataset descriptions 
X 

pfprintlist 
Print list of datasets and their labels 
X 

pfnewlabel 
Change dataset label 
X 

pfnewdata 
Create new dataset 
X 

pfprintgrid 
Print grid 
X 

pfnewgrid 
Set grid for new dataset 
X 

pfdelete 
Delete dataset 
X 

pfreload 
Reload dataset 
X 

pfreloadall 
Reload all current datasets 
X 

pfprintdata 
Print all elements of a dataset 
X 

pfprintelt 
Print a single element 
X 
Name 
Short Description 
Examples 
Compatible with TFG? 

File Operations 

pfload 
Load file 
All 
X 
pfloadsds 
Load Scientific Data Set from HDF file 
X 

pfdist 
Distribute files based on processor topology 
4 
X 
pfdistondomain 
Distribute files based on domain 
X 

pfundist 
Undistribute files 
X 

pfsave 
Save dataset 
1,2,5,6 
X 
pfsavesds 
Save dataset in an HDF format 
X 

pfvtksave 
Save dataset in VTK format using DEM 
X 
X 
pfwritedb 
Write the settings for a PF run to a database 
X 

Solid file operations 

pfpatchysolid 
Build a solid file between two complex surfaces and assign userdefined patches around the edges 
X 

pfs olidfmtconvert 
Converts back and forth between ascii and binary formats for solid files 
X 
Detailed descriptions of every command are included below in alphabetical order. Note that the required inputs are listed following each command. Commands that perform operations on data sets will require an identifier for each data set it takes as input. Inputs listed in square brackets are optional and do not need to be provided.
pfaxpy alpha x y
This command computes y = alpha*x+y where alpha is a scalar and x and y are identifiers representing data sets. No data set identifier is returned upon successful completion since data set y is overwritten.
pfbfcvel conductivity phead
This command computes the block face centered flow velocity at every grid cell. Conductivity and pressure head data sets are given as arguments. The output includes x, y, and z velocity components that are appended to the Tcl result.
pfbuilddomain database
This command builds a subgrid array given a ParFlow database that contains the domain parameters and the processor topology.
pfcelldiff datasetx datasety mask
This command computes cellwise differences of two datasets (diff=datasetxdatasety). This is the difference at each individual cell, not over the domain. Datasets must have the same dimensions.
pfcelldiffconst dataset constant mask
This command subtracts a constant value from each (active) cell of dataset (dif=dataset  constant).
pfcelldiv datasetx datasety mask
This command computes the cellwise quotient of datasetx and datasety (div = datasetx/datasety). This is the quotient at each individual cell. Datasets must have the same dimensions.
pfcelldivconst dataset constant mask
This command divides each (active) cell of dataset by a constant (div=dataset/constant).
pfcellmult datasetx datasety mask
This command computes the cellwise product of datasetx and datasety (mult = datasetx * datasety). This is the product at each individual cell. Datasets must have the same dimensions.
pfcellmultconst dataset constant mask
This command multiplies each (active) cell of dataset by a constant (mult=dataset * constant).
pfcellsum datasetp datasetq mask
This command computes the cellwise sum of two datasets (i.e., the sum at each individual cell, not the sum over the domain). Datasets must have the same dimensions.
pfcellsumconst dataset constant mask
This command adds the value of constant to each (active) cell of dataset.
pfchildD8 dem
This command computes the unique D8 child for all cells. Child[i,j] is the elevation of the cell to which [i,j] drains (i.e. the elevation of [i,j]’s child). If [i,j] is a local minima the child elevation set the elevation of [i,j].
pfcomputebottom mask
This command computes the bottom of the domain based on the mask of active and inactive zones. The identifier of the data set created by this operation is returned upon successful completion.
pfcomputedomain top bottom
This command computes a domain based on the top and bottom data sets. The domain built will have a single subgrid per processor that covers the active data as defined by the top and botttom. This domain will more closely follow the topology of the terrain than the default single computation domain.
A typical usage pattern for this is to start with a mask file (zeros in inactive cells and nonzero in active cells), create the top and bottom from the mask, compute the domain and then write out the domain. Refer to example number 3 in the following section.
pfcomputetop mask
This command computes the top of the domain based on the mask of
active and inactive zones. This is the landsurface in clm
or overland flow simulations. The identifier of the data set created
by this operation is returned upon successful completion.
pfcvel conductivity phead
This command computes the Darcy velocity in cells for the conductivity data set represented by the identifier ‘conductivity’ and the pressure head data set represented by the identifier ‘phead’. (note: This “cell” is not the same as the grid cells; its corners are defined by the grid vertices.) The identifier of the data set created by this operation is returned upon successful completion.
pfdelete dataset
This command deletes the data set represented by the identifier ‘dataset’. This command can be useful when working with multiple datasets / time series, such as those created when many timesteps of a file are loaded and processed. Deleting these datasets in between reads can help with tcl memory management.
pfdiffelt datasetp datasetq i j k digits [zero]
This command returns the difference of two corresponding coordinates from ‘datasetp’ and ‘datasetq’ if the number of digits in agreement (significant digits) differs by more than ‘digits’ significant digits and the difference is greater than the absolute zero given by ‘zero’.
pfdist [options] filename
Distribute the file onto the virtual file system. This utility must be used to create files which ParFlow can use as input. ParFlow uses a virtual file system which allows each node of the parallel machine to read from the input file independently. The utility does the inverse of the pfundist command. If you are using a ParFlow binary file for input you should do a pfdist just before you do the pfrun. This command requires that the processor topology and computational grid be set in the input file so that it knows how to distribute the data. Note that the old syntax for pfdist required the NZ key be set to 1 to indicate a two dimensional file but this can now be specified manually when pfdist is called by using the optional argument nz followed by the number of layers in the file to be distributed, then the filename. If the nz argument is absent the NZ key is used by default for the processor topology.
For example,
pfdist nz 1 slopex.pfb
pfdistondomain filename domain
Distribute the file onto the virtual file system based on the domain provided rather than the processor topology as used by pfdist. This is used by the SAMRAI version of which allows for a more complicated computation domain specification with different sized subgrids on each processor and allows for more than one subgrid per processor. Frequently this will be used with a domain created by the pfcomputedomain command.
pfeffectiverecharge precip et slopex slopey dem
This command computes the effective recharge at every grid cell based on total precipitation minus evapotranspiration (PET)in the upstream area. Effective recharge is consistent with TOPMODEL definition, NOT local PET. Inputs are total annual (or average annual) precipitation (precip) at each point, total annual (or average annual) evapotranspiration (ET) at each point, slope in the x direction, slope in the y direction and elevation.
pfenlargebox dataset sx sy sz
This command returns a new dataset which is enlarged to be of the new size indicated by sx, sy and sz. Expansion is done first in the z plane, then y plane and x plane.
pfextract2Ddomain domain
This command builds a 2D domain based off a 3D domain. This can be used for a pfdistondomain command for Parflow 2D data (such as slopes and soil indices).
pfextracttop top data
This command computes the top of the domain based on the top of the domain and another dataset. The identifier of the data set created by this operation is returned upon successful completion.
pffillflats dem
This command finds the flat regions in the DEM and eliminates them by bilinearly interpolating elevations across flat region.
pfflintslaw dem c p
This command smooths the digital elevation model dem according to Flints Law, with Flints Law parameters specified by c and p, respectively. Flints Law relates the slope magnitude at a given cell to its upstream contributing area: S = c*A**p. In this routine, elevations at local minima retain the same value as in the original dem. Elevations at all other cells are computed by applying Flints Law recursively up each drainage path, starting at its terminus (a local minimum) until a drainage divide is reached. Elevations are computed as:
dem[i,j] = dem[child] + c*(A[i,j]**p)*ds[i,j]
where child is the D8 child of [i,j] (i.e., the cell to which [i,j] drains according to the D8 method); ds[i,j] is the segment length between the [i,j] and its child; A[i,j] is the upstream contributing area of [i,j]; and c and p are constants.
pfflintslawbybasin dem c0 p0 maxiter
This command smooths the digital elevation model (dem) using the same approach as “pfflints law”. However here the c and p parameters are fit for each basin separately. The Flint’s Law parameters are calculated for the provided digital elevation model dem using the iterative LevenbergMarquardt method of nonlinear least squares minimization, as in “pfflintslawfit”. The user must provide initial estimates of c0 and p0; results are not sensitive to these initial values. The user must also specify the maximum number of iterations as maxiter.
pfflintslawfit dem c0 p0 maxiter
This command fits Flint’s Law parameters c and p for the provided digital elevation model dem using the iterative LevenbergMarquardt method of nonlinear least squares minimization. The user must provide initial estimates of c0 and p0; results are not sensitive to these initial values. The user must also specify the maximum number of iterations as maxiter. Final values of c and p are printed to the screen, and a dataset containing smoothed elevation values is returned. Smoothed elevations are identical to running pfflintslaw for the final values of c and p. Note that dem must be a ParFlow dataset and must have the correct grid information – dx, dy, nx, and ny are used in parameter estimation and Flint’s Law calculations. If gridded elevation values are read in from a text file (e.g., using pfload’s simple ascii format), grid information must be specified using the pfsetgrid command.
pfflux conductivity hhead
This command computes the net Darcy flux at vertices for the conductivity data set ‘conductivity’ and the hydraulic head data set given by ‘hhead’. An identifier representing the flux computed will be returned upon successful completion.
pfgetelt dataset i j k
This command returns the value at element (i,j,k) in data set ‘dataset’. The i, j, and k above must range from 0 to (nx  1), 0 to (ny  1), and 0 to (nz  1) respectively. The values nx, ny, and nz are the number of grid points along the x, y, and z axes respectively. The string ‘dataset’ is an identifier representing the data set whose element is to be retrieved.
pfgetgrid dataset
This command returns a description of the grid which serves as the domain of data set ‘dataset’. The format of the description is given below.
(nx, ny, nz)
The number of coordinates in each direction.
(x, y, z)
The origin of the grid.
(dx, dy, dz)
The distance between each coordinate in each direction.
The above information is returned in the following Tcl list format: nx ny nz x y z dx dy dz
pfgetlist dataset
This command returns the name and description of a dataset if an argument is provided. If no argument is given, then all of the data set names followed by their descriptions is returned to the TCL interpreter. If an argument (dataset) is given, it should be the it should be the name of a loaded dataset.
pfgetstats dataset
This command calculates the following statistics for the data set represented by the identifier dataset: minimum, maximum, mean, sum, variance, and standard deviation.
pfgetsubbox dataset il jl kl iu ju ku
This command computes a new dataset with the subbox starting at il, jl, kl and going to iu, ju, ku.
pfgridtype gridtype
This command sets the grid type to either cell centered if ‘gridtype’ is set to ‘cell’ or vetex centered if ‘gridtype’ is set to ‘vertex’. If no new value for ‘gridtype’ is given, then the current value of ‘gridtype’ is returned. The value of ‘gridtype’ will be returned upon successful completion of this command.
pfgwstorage mask porosity pressure saturation specific_storage
This command computes the subsurface water storage (compressible and incompressible components) based on mask, porosity, saturation, storativity and pressure fields, similar to pfsubsurfacestorage, but only for the saturated cells.
pfhelp [command]
This command returns a list of pftools commands. If a command is provided it gives a detailed description of the command and the necessary inputs.
pfhhead phead
This command computes the hydraulic head from the pressure head represented by the identifier ‘phead’. An identifier for the hydraulic head computed is returned upon successful completion.
pfhydrostatic wtdepth top mask
Compute hydrostatic pressure field from water table depth
pflistdata dataset
This command returns a list of pairs if no argument is given. The first item in each pair will be an identifier representing the data set and the second item will be that data set’s label. If a data set’s identifier is given as an argument, then just that data set’s name and label will be returned.
pfload [file format] filename
Loads a dataset into memory so it can be manipulated using the other utilities. A file format may preceed the filename in order to indicate the file’s format. If no file type option is given, then the extension of the filename is used to determine the default file type. An identifier used to represent the data set will be returned upon successful completion.
File type options include:
pfb
ParFlow binary format. Default file type for files with a ‘.pfb’ extension.
pfsb
ParFlow scattered binary format. Default file type for files with a ‘.pfsb’ extension.
sa
ParFlow simple ASCII format. Default file type for files with a ‘.sa’ extension.
sb
ParFlow simple binary format. Default file type for files with a ‘.sb’ extension.
silo
Silo binary format. Default file type for files with a ‘.silo’ extension.
rsa
ParFlow real scattered ASCII format. Default file type for files with a ‘.rsa’ extension
pfloadsds filename dsnum
This command is used to load Scientific Data Sets from HDF files. The SDS number ‘dsnum’ will be used to find the SDS you wish to load from the HDF file ‘filename’. The data set loaded into memory will be assigned an identifier which will be used to refer to the data set until it is deleted. This identifier will be returned upon successful completion of the command.
pfmdiff datasetp datasetq digits [zero]
If ‘digits’ is greater than or equal to zero, then this command computes the grid point at which the number of digits in agreement (significant digits) is fewest and differs by more than ‘digits’ significant digits. If ‘digits’ is less than zero, then the point at which the number of digits in agreement (significant digits) is minimum is computed. Finally, the maximum absolute difference is computed. The above information is returned in a Tcl list of the following form: mi mj mk sd adiff
Given the search criteria, (mi, mj, mk) is the coordinate where the minimum number of significant digits ‘sd’ was found and ‘adiff’ is the maximum absolute difference.
pfmovingaveragedem dem wsize maxiter
This command fills sinks in the digital elevation model dem by a standard iterative movingaverage routine. Sinks are identified as cells with zero slope in both x and ydirections, or as local minima in elevation (i.e., all adjacent cells have higher elevations). At each iteration, a moving average is taken over a window of width wsize around each remaining sink; sinks are thus filled by averaging over neighboring cells. The procedure continues iteratively until all sinks are filled or the number of iterations reaches maxiter. For most applications, sinks should be filled prior to computing slopes (i.e., prior to executing pfslopex and pfslopey).
pfnewdata {nx ny nz} {x y z} {dx dy dz} label
This command creates a new data set whose dimension is described by the lists nx ny nz, x y z, and dx dy dz. The first list, describes the dimensions, the second indicates the origin, and the third gives the length intervals between each coordinate along each axis. The ‘label’ argument will be the label of the data set that gets created. This new data set that is created will have all of its data points set to zero automatically. An identifier for the new data set will be returned upon successful completion.
pfnewgrid {nx ny nz} {x y z} {dx dy dz} label
Create a new data set whose grid is described by passing three lists and a label as arguments. The first list will be the number of coordinates in the x, y, and z directions. The second list will describe the origin. The third contains the intervals between coordinates along each axis. The identifier of the data set created by this operation is returned upon successful completion.
pfnewlabel dataset newlabel
This command changes the label of the data set ‘dataset’ to ‘newlabel’.
pfphead hhead
This command computes the pressure head from the hydraulic head represented by the identifier ‘hhead’. An identifier for the pressure head is returned upon successful completion.
pfpatchysolid top topdata bot botdata msk emaskdata [optional args]
Creates a solid file with complex upper and lower surfaces from a top surface elevation dataset (topdata), a bottom elevation dataset (botdata), and an enhanced mask dataset (emaskdata) all of which must be passed as handles to 2d datasets that share a common size and origin. The solid is built as the volume between the top and bottom surfaces using the mask to deactivate other regions. The “enhanced mask” used here is a gridded dataset containing integers where all active cells have values of one but inactive cells may be given a positive integer value that identifies a patch along the model edge or the values may be zero. Any mask cell with value 0 is omitted from the active domain and is not written to a patch. If an active cell is adjacent to a nonzero mask cell, the face between the active and inactive cell is assigned to the patch with the integer value of the adjacent inactive cell. Bottom and Top patches are always written for every active cell and the West, East, South, and North edges are written automatically anytime active cells touch the edges of the input dataset(s). Up to 30 user defined patches can be specified using arbitrary integer values that are greater than 1. Note that the msk flag may be omitted and doing so will make every cell active.
The top and bot flags, and msk if it is used, MUST each be followed by the handle for the relevant dataset. Optional argument flagname pairs include:
pfsol <file name>.pfsol (or pfsolb <file name>.pfsolb)
vtk <file name>.vtk
sub
where <file name> is replaced by the desired text string. The pfsolb option creates a compact binary solid file; pfsolb cannot currently be read directly by ParFlow but it can be converted with pfsolidfmtconvert and full support is under development. If pfsol (or pfsolb) is not specified the default name “SolidFile.pfsol” will be used. If vtk is omitted, no vtk file will be created. The vtk attributes will contain mean patch elevations and patch IDs from the enhanced mask. Edge patch IDs are shown as negative values in the vtk. The patchysolid tool also outputs the list of the patch names in the order they are written, which can be directly copied into a ParFlow TCL script for the list of patch names. The sub option writes separate patches for each face (left,right,front,back), which are indicated in the output patch write order list.
Assuming $Msk, $Top, and $Bot are valid dataset handles from pfload, two valid examples are:
pfpatchysolid msk $Msk top $Top bot $Bot pfsol "MySolid.pfsol" vtk "MySolid.vtk"
pfpatchysolid bot $Bot top $Top vtk "MySolid.vtk" sub
Note that all flagname pairs may be specified in any order for this tool as long as the required argument immediately follows the flag. To use with a terrain following grid, you will need to subtract the surface elevations from the top and bottom datasets (this makes the top flat) then add back in the total thickness of your grid, which can be done using “pfcelldiff” and “pfcellsumconst”.
pfpitfilldem dem dpit maxiter
This command fills sinks in the digital elevation model dem by a standard iterative pitfilling routine. Sinks are identified as cells with zero slope in both x and ydirections, or as local minima in elevation (i.e., all adjacent neighbors have higher elevations). At each iteration, the value dpit is added to all remaining sinks. The procedure continues iteratively until all sinks are filled or the number of iterations reaches maxiter. For most applications, sinks should be filled prior to computing slopes (i.e., prior to executing pfslopex and pfslopey).
pfprintdata dataset
This command executes ‘pfgetgrid’ and ‘pfgetelt’ in order to display all the elements in the data set represented by the identifier ‘dataset’.
pfprintdiff datasetp datasetq digits [zero]
This command executes ‘pfdiffelt’ and ‘pfmdiff’ to print differences to standard output. The differences are printed one per line along with the coordinates where they occur. The last two lines displayed will show the point at which there is a minimum number of significant digits in the difference as well as the maximum absolute difference.
pfprintdomain domain
This command creates a set of TCL commands that setup a domain as specified by the provided domain input which can be then be written to a file for inclusion in a Parflow input script. Note that this kind of domain is only supported by the SAMRAI version of Parflow.
pfprintelt i j k dataset
This command prints a single element from the provided dataset given an i, j, k location.
pfprintgrid dataset
This command executes pfgetgrid and formats its output before printing it on the screen. The triples (nx, ny, nz), (x, y, z), and (dx, dy, dz) are all printed on seperate lines along with labels describing each.
pfprintlist [dataset]
This command executes pflistdata and formats the output of that command. The formatted output is then printed on the screen. The output consists of a list of data sets and their labels one per line if no argument was given or just one data set if an identifier was given.
pfprintmdiff datasetp datasetq digits [zero]
This command executes ‘pfmdiff’ and formats that command’s output before displaying it on the screen. Given the search criteria, a line displaying the point at which the difference has the least number of significant digits will be displayed. Another line displaying the maximum absolute difference will also be displayed.
printstats dataset
This command executes ‘pfstats’ and formats that command’s output before printing it on the screen. Each of the values mentioned in the description of ‘pfstats’ will be displayed along with a label.
pfreload dataset
This argument reloads a dataset. Only one arguments is required, the name of the dataset to reload.
pfreloadall
This command reloads all of the current datasets.
pfsattrans mask perm
Compute saturated transmissivity for all [i,j] as the sum of the permeability[i,j,k]*dz within a column [i,j]. Currently this routine uses dz from the input permeability so the dz in permeability must be correct. Also, it is assumed that dz is constant, so this command is not compatible with variable dz.
pfsave dataset filetype filename
This command is used to save the data set given by the identifier ‘dataset’ to a file ‘filename’ of type ‘filetype’ in one of the ParFlow formats below.
File type options include:
pfb ParFlow binary format.
sa ParFlow simple ASCII format.
sb ParFlow simple binary format.
silo Silo binary format.
vis Vizamrai binary format.
pfsavediff datasetp datasetq digits [zero] file filename
This command saves to a file the differences between the values of the data sets represented by ‘datasetp’ and ‘datasetq’ to file ‘filename’. The data points whose values differ in more than ‘digits’ significant digits and whose differences are greater than ‘zero’ will be saved. Also, given the above criteria, the minimum number of digits in agreement (significant digits) will be saved.
If ‘digits’ is less than zero, then only the minimum number of significant digits and the coordinate where the minimum was computed will be saved.
In each of the above cases, the maximum absolute difference given the criteria will also be saved.
pfsavesds dataset filetype filename
This command is used to save the data set represented by the identifier ‘dataset’ to the file ‘filename’ in the format given by ‘filetype’.
The possible HDF formats are:
float32
float64
int8
uint8
int16
uint16
int32
uint32
pfsegmentD8 dem
This command computes the distance between the cell centers of every parent cell [i,j] and its child cell. Child cells are determined using the eightpoint pour method (commonly referred to as the D8 method) based on the digital elevation model dem. If [i,j] is a local minima the segment length is set to zero.
pfsetgrid {nx ny nz} {x0 y0 z0} {dx dy dz} dataset
This command replaces the grid information of dataset with the values provided.
pfslopeD8 dem
This command computes slopes according to the eightpoint pour method (commonly referred to as the D8 method) based on the digital elevation model dem. Slopes are computed as the maximum downward gradient between a given cell and it’s lowest neighbor (adjacent or diagonal). Local minima are set to zero; where local minima occur on the edge of the domain, the 1st order upwind slope is used (i.e., the cell is assumed to drain out of the domain). Note that dem must be a ParFlow dataset and must have the correct grid information – dx and dy both used in slope calculations. If gridded elevation values are read in from a text file (e.g., using pfload’s simple ascii format), grid information must be specified using the pfsetgrid command. It should be noted that ParFlow uses slopex and slopey (NOT D8 slopes!) in runoff calculations.
pfslopex dem
This command computes slopes in the xdirection using 1st order upwind finite differences based on the digital elevation model dem. Slopes at local maxima (in xdirection) are calculated as the maximum downward gradient to an adjacent neighbor. Slopes at local minima (in xdirection) do not drain in the xdirection and are therefore set to zero. Note that dem must be a ParFlow dataset and must have the correct grid information – dx in particular is used in slope calculations. If gridded elevation values are read from a text file (e.g., using pfload’s simple ascii format), grid inforamtion must be specified using the pfsetgrid command.
pfslopexD4 dem
This command computes the slope in the xdirection for all [i,j] using a four point (D4) method. The slope is set to the maximum downward slope to the lowest adjacent neighbor. If [i,j] is a local minima the slope is set to zero (i.e. no drainage).
pfslopey dem
This command computes slopes in the ydirection using 1st order upwind finite differences based on the digital elevation model dem. Slopes at local maxima (in ydirection) are calculated as the maximum downward gradient to an adjacent neighbor. Slopes at local minima (in ydirection) do not drain in the ydirection and are therefore set to zero. Note that dem must be a ParFlow dataset and must have the correct grid information  dy in particular is used in slope calculations. If gridded elevation values are read in from a text file (e.g., using pfload’s simple ascii format), grid information must be specified using the pfsetgrid command.
pfslopeyD4 dem
This command computes the slope in the ydirection for all [i,j] using a four point (D4) method. The slope is set to the maximum downward slope to the lowest adjacent neighbor. If [i,j] is a local minima the slope is set to zero (i.e. no drainage).
pfsolidfmtconvert filename1 filename2
This command converts solid files back and forth between the ascii .pfsol format and the binary .pfsolb format. The tool automatically detects the conversion mode based on the extensions of the input file names. The filename1 is the name of source file and filename2 is the target output file to be created or overwritten. Support to directly use a binary solid (.pfsolb) is under development but this allows a significant reduction in file sizes.
For example, to convert from ascii to binary, then back to ascii:
pfsolidfmtconvert "MySolid.pfsol" "MySolid.pfsolb"
pfsolidfmtconvert "MySolid.pfsolb" "NewSolid.pfsol"
pfstats dataset
This command prints various statistics for the data set represented by the identifier ‘dataset’. The minimum, maximum, mean, sum, variance, and standard deviation are all computed. The above values are returned in a list of the following form: min max mean sum variance (standard deviation)
pfsubsurfacestorage mask porosity pressure saturation specific_storage
This command computes the subsurface water storage (compressible and incompressible components) based on mask, porosity, saturation, storativity and pressure fields. The equations used to calculate this quantity are given in §5.9 Water Balance. The identifier of the data set created by this operation is returned upon successful completion.
pfsum dataset
This command computes the sum over the domain of the dataset.
pfsurfacerunoff top slope_x slope_y mannings pressure
This command computes the surface water runoff (out of the domain) based on a computed top, pressure field, slopes and mannings roughness values. This is integrated along all domain boundaries and is calculated at any location that slopes at the edge of the domain point outward. This data is in units of \([L^3 T^{1}]\) and the equations used to calculate this quantity are given in §5.9 Water Balance. The identifier of the data set created by this operation is returned upon successful completion.
pfsurfacestorage top pressure
This command computes the surface water storage (ponded water on top of the domain) based on a computed top and pressure field. The equations used to calculate this quantity are given in §5.9 Water Balance. The identifier of the data set created by this operation is returned upon successful completion.
pftopodeficit profile m trans dem slopex slopey recharge ssat sres porosity mask
Compute water deficit for all [i,j] based on TOPMODEL/topographic index. For more details on methods and assumptions refer to toposlopes.c in pftools.
pftopoindex dem sx sy
Compute topographic index for all [i,j]. Here topographic index is defined as the total upstream area divided by the contour length, divided by the local slope. For more details on methods and assumptions refer to toposlopes.c in pftools.
pftoporecharge riverfile nriver trans dem sx sy
Compute effective recharge at all [i,j] over upstream area based on topmodel assumptions and given list of river points. Notes: See detailed notes in toposlopes.c regarding assumptions, methods, etc. Input Notes: nriver is an integer (number of river points) river is an array of integers [nriver][2] (list of river indices, ordered from outlet to headwaters) is a Databox of saturated transmissivity dem is a Databox of elevations at each cell sx is a Databox of slopes (xdir) – lets you use processed slopes! sy is a Databox of slopes (ydir) – lets you use processed slopes!
pftopowt deficit porosity ssat sres mask top wtdepth
Compute water depth from column water deficit for all [i,j] based on TOPMODEL/topographic index.
pfundist filename, pfundist runname
The command undistributes a ParFlow output file. ParFlow uses a distributed file system where each node can write to its own file. The pfundist command takes all of these individual files and collapses them into a single file.
The arguments can be a runname or a filename. If a runname is given then all of the output files associated with that run are undistributed.
Normally this is done after every pfrun command.
pfupstreamarea slope_x slope_y
This command computes the upstream area contributing to surface
runoff at each cell based on the x and y slope values provided in
datasets slope_x
and slope_y
, respectively. Contributing
area is computed recursively for each cell; areas are not weighted
by slope direction. Areas are returned as the number of upstream
(contributing) cells; to compute actual area, simply multiply by
the cell area (dx*dy).
pfvmag datasetx datasety datasetz
This command computes the velocity magnitude when given three velocity components. The three parameters are identifiers which represent the x, y, and z components respectively. The identifier of the data set created by this operation is returned upon successful completion.
pfvtksave dataset filetype filename [options]
This command loads PFB or SILO output, reads a DEM from a file and generates a 3D VTK output field from that ParFlow output.
The options: Any combination of these can be used and they can be specified in any order as long as the required elements immediately follow each option.
var specifies what the variable written to the dataset will be called. This is followed by a text string, like “Pressure” or “Saturation” to define the name of the data that will be written to the VTK. If this isn’t specified, you’ll get a property written to the file creatively called “Variable”. This option is ignored if you are using clmvtk since all its variables are predefined.
dem specifies that a DEM is to be used. The argument following dem MUST be the handle of the dataset containing the elevations. If it cannot be found, the tool ignores it and reverts to nondem mode. If the nx and ny dimensions of the grids don’t match, the tool will error out. This option shifts the layers so that the top of the domain coincides with the land surface defined by the DEM. Regardless of the actual number of layers in the DEM file, the tool only uses the elevations in the top layer of this dataset, meaning a 1layer PFB can be used.
flt tells the tool to write the data as type float instead of double. Since the VTKs are really only used for visualization, this reduces the file size and speeds up plotting.
tfg causes the tool to override the specified dz in the dataset PFB and uses a user specified list of layer thicknesses instead. This is designed for terrain following grids and can only be used in conjunction with a DEM. The argument following the flag is a text string containing the number of layers and the dz list of actual layer thicknesses (not dz multipliers) for each layer from the bottom up such as: tfg “5 200.0 1.0 0.7 0.2 0.1” Note that the quotation marks around the list are necessary.
Example:
file copy force CLM_dem.cpfb CLM_dem.pfb
set CLMdat [pfload pfb clm.out.clm_output.00005.C.pfb]
set Pdat [pfload pfb clm.out.press.00005.pfb]
set Perm [pfload pfb clm.out.perm_x.pfb]
set DEMdat [pfload pfb CLM_dem.pfb]
set dzlist "10 6.0 5.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5"
pfvtksave $Pdat vtk "CLM.out.Press.00005a.vtk" var "Press"
pfvtksave $Pdat vtk "CLM.out.Press.00005b.vtk" var "Press" flt
pfvtksave $Pdat vtk "CLM.out.Press.00005c.vtk" var "Press" dem $DEMdat
pfvtksave $Pdat vtk "CLM.out.Press.00005d.vtk" var "Press" dem $DEMdat flt
pfvtksave $Pdat vtk "CLM.out.Press.00005e.vtk" var "Press" dem $DEMdat flt tfg $dzlist
pfvtksave $Perm vtk "CLM.out.Perm.00005.vtk" var "Perm" flt dem $DEMdat tfg $dzlist
pfvtksave $CLMdat clmvtk "CLM.out.CLM.00005.vtk" flt
pfvtksave $CLMdat clmvtk "CLM.out.CLM.00005.vtk" flt dem $DEMdat
pfvtksave $DEMdat vtk "CLM.out.Elev.00000.vtk" flt var "Elevation" dem $DEMdat
pfvvel conductivity phead
This command computes the Darcy velocity in cells for the conductivity data set represented by the identifier ‘conductivity’ and the pressure head data set represented by the identifier ‘phead’. The identifier of the data set created by this operation is returned upon successful completion.
pfwatertabledepth top saturation
This command computes the water table depth (distance from top to first cell with saturation = 1). The identifier of the data set created by this operation is returned upon successful completion.
pfwritedb runname
This command writes the settings of parflow run to a pfidb database that can be used to run the model at a later time. In general this command is used in lieu of the pfrun command.
7.3. Common examples using ParFlow TCL commands (PFTCL)¶
This section contains some brief examples of how to use the pftools commands (along with standard TCL commands) to postprocess data.
Load a file as one format and write as another format.
set press [pfload harvey_flow.out.press.pfb]
pfsave $press sa harvey_flow.out.sa
#####################################################################
# Also note that PFTCL automatically assigns
#identifiers to each data set it stores. In this
# example we load the pressure file and assign
#it the identifier press. However if you
#read in a file called foo.pfb into a TCL shell
#with assigning your own identifier, you get
#the following:
#parflow> pfload foo.pfb
#dataset0
# In this example, the first line is typed in by the
#user and the second line is printed out
#by PFTCL. It indicates that the data read
#from file foo.pfb is associated with the
#identifier dataset0.
Load pressurehead output from a file, convert to headpotential and write out as a new file.
set press [pfload harvey_flow.out.press.pfb]
set head [pfhhead $press]
pfsave $head pfb harvey_flow.head.pfb
Build a SAMARI compatible domain decomposition based off of a mask file.
#
# This example script takes 3 command line arguments
# for P,Q,R and then builds a SAMRAI compatible
# domain decomposition based off of a mask file.
#
# Processor Topology
set P [lindex $argv 0]
set Q [lindex $argv 1]
set R [lindex $argv 2]
pfset Process.Topology.P $P
pfset Process.Topology.Q $Q
pfset Process.Topology.R $R
# Computational Grid
pfset ComputationalGrid.Lower.X 10.0
pfset ComputationalGrid.Lower.Y 10.0
pfset ComputationalGrid.Lower.Z 1.0
pfset ComputationalGrid.DX 8.8888888888888893
pfset ComputationalGrid.DY 10.666666666666666
pfset ComputationalGrid.DZ 1.0
pfset ComputationalGrid.NX 10
pfset ComputationalGrid.NY 10
pfset ComputationalGrid.NZ 8
# Calculate top and bottom and build domain
set mask [pfload samrai.out.mask.pfb]
set top [pfcomputetop $mask]
set bottom [pfcomputebottom $mask]
set domain [pfcomputedomain $top $bottom]
set out [pfprintdomain $domain]
set grid\_file [open samrai_grid.tcl w]
puts $grid_file $out
close $grid_file
#
# The resulting TCL file samrai_grid.tcl may be read into
# a Parflow input file using ¿¿source samrai_grid.tcl¿¿.
#
Distributing input files before running [dist example]
#
# A common problem for new ParFlow users is to
# distribute slope files using
# the 3D computational grid that is
# set at the begging of a run script.
# This results in errors because slope
# files are 2D.
# To avoid this problem the computational
# grid should be reset before and after
# distributing slope files. As follows:
#
#First set NZ to 1 and distribute the 2D slope files
pfset ComputationalGrid.NX 40
pfset ComputationalGrid.NY 40
pfset ComputationalGrid.NZ 1
pfdist slopex.pfb
pfdist slopey.pfb
#Reset NZ to the correct value and distribute any 3D inputs
pfset ComputationalGrid.NX 40
pfset ComputationalGrid.NY 40
pfset ComputationalGrid.NZ 50
pfdist IndicatorFile.pfb
Calculate slopes from an elevation file
#Read in DEM
set dem [pfload sa dem.txt]
pfsetgrid {209 268 1} {0.0 0.0 0.0} {100 100 1.0} $dem
# Fill flat areas (if any)
set flatfill [pffillflats $dem]
# Fill pits (if any)
set pitfill [pfpitfilldem $flatfill 0.01 10000]
# Calculate Slopes
set slope_x [pfslopex $pitfill]
set slope_y [pfslopey $pitfill]
# Write to output...
pfsave $flatfill silo klam.flatfill.silo
pfsave $pitfill silo klam.pitfill.silo
pfsave $slope_x pfb klam.slope_x.pfb
pfsave $slope_y pfb klam.slope_y.pfb
Calculate and output the subsurface storage in the domain at a point in time.
set saturation [pfload runname.out.satur.00001.silo]
set pressure [pfload runname.out.press.00001.silo]
set specific_storage [pfload runname.out.specific_storage.silo]
set porosity [pfload runname.out.porosity.silo]
set mask [pfload runname.out.mask.silo]
set subsurface_storage [pfsubsurfacestorage $mask $porosity \
$pressure $saturation $specific_storage]
set total_subsurface_storage [pfsum $subsurface_storage]
puts [format "Subsurface storage\t\t\t\t : %.16e" $total_subsurface_storage]
Calculate and output the surface storage in the domain at a point in time.
set pressure [pfload runname.out.press.00001.silo]
set mask [pfload runname.out.mask.silo]
set top [pfcomputetop $mask]
set surface_storage [pfsurfacestorage $top $pressure]
set total_surface_storage [pfsum $surface_storage]
puts [format "Surface storage\t\t\t\t : %.16e" $total_surface_storage]
Calculate and output the runoff out of the entire domain over a timestep.
set pressure [pfload runname.out.press.00001.silo]
set slope_x [pfload runname.out.slope_x.silo]
set slope_y [pfload runname.out.slope_y.silo]
set mannings [pfload runname.out.mannings.silo]
set mask [pfload runname.out.mask.silo]
set top [pfcomputetop $mask]
set surface_runoff [pfsurfacerunoff $top $slope_x $slope_y $mannings $pressure]
set total_surface_runoff [expr [pfsum $surface_runoff] * [pfget TimeStep.Value]]
puts [format "Surface runoff from pftools\t\t\t : %.16e" $total_surface_runoff]
Calculate overland flow at a point using Manning’s equation
#Set the location
set Xloc 2
set Yloc 2
set Zloc 50 #This should be a z location on the surface of your domain
#Set the grid dimension and Mannings roughness coefficient
set dx 1000.0
set n 0.000005
#Get the slope at the point
set slopex [pfload runname.out.slope_x.pfb]
set slopey [pfload runname.out.slope_y.pfb]
set sx1 [pfgetelt $slopex $Xloc $Yloc 0]
set sy1 [pfgetelt $slopey $Xloc $Yloc 0]
set S [expr ($sx**2+$sy**2)**0.5]
#Get the pressure at the point
set press [pfload runname.out.press.00001.pfb]
set P [pfgetelt $press $Xloc $Yloc $Zloc]
#If the pressure is less than zero set to zero
if {$P < 0} { set P 0 }
set QT [expr ($dx/$n)*($S**0.5)*($P**(5./3.))]
puts $QT