Welcome to ctraj! The purpose of these codes is to perform 2-D, Lagrangian advection. Input to the codes can be any source of wind data or pre-generated velocity fields, since the velocity vectors are first converted to an intermediate file format for easy data-abstraction. As such, the codes require two subsidiary libraries, one for the data structures containing the wind fields and, since the semi-Lagrangian tracer simulations generate a matrix as intermediate output, a sparse matrix library. These libraries are available with the "libpetey" numerical library, available from the sourceforge project, "msci".
INSTALLATION
1. Download the files. They can be obtained either from the sourceforge download page as a gzipped tarball or more directly from the Subversion repository. If you want the very latest version or want to be a developer, take the latter option. Four external libraries are required: the Gnu Scientific Libary (GSL), NetCDF (available from UCAR), ARPACK (Arnoldi package for sparse eigensystems) and a smaller scientific library called libpetey. The latter can be downloaded from the Sourceforge page. For plotting you will also need the Generic Mapping Tools (GMT) as well as a postscript viewer, eg, ghostview.
2. Modify the makefiles for your system. Important macros include:
CC: C++ compiler
CFLAGS: options to pass to the compiler
LD_OPTIONS: options to pass to the linker
LIB_PATH: Directory to store library files
BIN_PATH: Directory for binary executables
INCLUDE_PATH: Directory for include files
GSL_LIB: Location of the GSL libraries
GSL_INCLUDE: Location of the GSL include files
NCDF_LIB: Location of the NetCDF libraries
NCDF_INCLUDE: Location of the NetCDF include files
ARPATH: Location of the ARPACK library
LIBARPACK: Name of the ARPACK library
FORTRAN_RUNTIME: Name of the Fortran runtime library--for linking in ARPACK eigen-system codes.
Take particular note of the OPT_VER macro. This macro defines the level of optimization for each library. In the libpetey library, the flag is appended to all the object files, including the libraries and executables. This allows you to store multiple versions with different levels of optimization, which is particularly handy for debugging purposes. This is not the case for the ctraj library which simply uses the flag to decide which version of the external libraries to link in.
I do not use a configure script for several reasons. The configure script for simple programs is frequently longer than the program itself. This is NOT a recipe for portability! Moreover, configure scripts fail almost as often as makefiles and they are much harder to debug. This software is aimed at the advanced user; as such, he or she should be intimately aware of his or her machine and will want to closely customize the software the given application. This can be done more directly and more thoroughly using only a makefile.
3. Compile the programs. Type "make" in the base directory. Executables are copied automatically to your executable directory.
QUICKSTART
All programs are run directly from the command line. For a command synopsis, including the format of the initialization file, type the command with no arguments. If the command takes no arguments, use the "-?" switch.
The following programs are included with the software bundle:
ctraj_model: General trajectory model. Currently supports 2-D trajectories
on a global surface, 2-D trajectories on a Cartesian surface
and analytic trajectories in a Cartesian space.
ctraj_tracer: General tracer model. Outputs a series of sparse matrices
representing the tracer dynamics.
ctraj_contour: General contour advection model.
contour_dimension: Calculates the fractal dimension of an advected contour.
contour_length: Calculates the length of an advected contour.
lla2aeb: Converts ASCII tracer fields to the native binary format.
extract_field: Converts binary tracer fields to ASCII.
xy2bev: Converts ASCII contours to binary.
bev2xy: Converts binary contours to ASCII.
bev2bev2: Converts binary ".bev" files to a link-chained, paging file
compatible with contour3.
zonally_symmetric_tracer: Generates a zonally symmetric tracer field.
C0: Generates a circularly symmetric contour along a single line of latitude.
nc2ds: Converts NCEP files to "dataset" files. Converts
a single vertical level.
vdata2ds: Converts velocity fields in ASCII format to dataset files.
ecmwf2ds.sh: Downloads ECMWF data from the data server and converts it to
a dataset. Must have wgrib installed and an account with the
ECMWF.
pc_proxy: For "principal-component proxy" tracer analysis.
PC-proxy is a new method for dynamical reconstruction of
long-lived tracers such as ozone. For an explanation, please
see the paper, "Principal component proxy tracer analysis,"
arxiv:1202.1999.
Please reference this paper if you use PC-proxy in
your work. A peer-reviewed article is currently in the works.
pc_proxy_predict: PC-proxy reconstruction of multiple fields.
proxy_tracer: More traditional proxy tracer reconstruction (doesn't work).
random_global: Generates a series of random points, evenly distributed around
the globe.
select_meas: Selects a subset of sparse measurements from a file.
tracer_interpolate: Interpolates points within a tracer.
Executables accept some subset of the following options. These can also be found in the file "ctraj_optargs.cc". Square brackets denote defaults.
-- restrict to Southern hemisphere
-+ restrict to Northern hemisphere
-0 initial time index [0]
-1 index of record in first file [0]
-2 index of record in second file [0]
-a use 'corners' in tracer simulation
-A number of Arnoldi vectors [20]
-B data file page size in bytes
-c maximum degrees of arc between pairs of points [1.]
-C count the number of measurements in PC proxy
-d width of date field [23]
-D calculate fractal dimension of type:
u=uncertainty exponent
0=box-counting dimension
l=box-counting dimension using line segments
(multiple types are supported, e.g. -D u0l)
-e number of values of dependent variable (epsilon) [20]
-E generate meridionally symmetric field
-f final date [last available]
-F final value for dependent variable (epsilon/z) (1000.)
-g use geometric progression
-h non-dimensional 'coarse' time step [1.]
-H print header information
-i initial date [first available]
-I initial value for dependent variable (epsilon/z) [1.]
-k number of R-K 'fine' time steps per coarse step [6]
-K apply constant term in PC proxy
-l lead time for the measurement window [default is the same as for the tracer]
-m minimum change in path in km [10.]
-M maximum number of points in advected contour
-n grid-points per side [120/100]
-N number of time-grids/time-steps [to end of file]
-o wrap off
-O write interval [1]
-p data path [./]
-P calculate Pearson correlation coefficient
-q sidelengths of tracer field in the form: d1xd2xd3...
-Q query time grids (write to stdout)
-r sidelength/2 [12000/10000]
-R (return) record size
-s maximum change in path in km [100/20]
-S sort measurement data by date
-T use theta (potential temperature) levels
-u epsilon-uncertain: minimum number [100]
-U epsilon-uncertain: maximum trials [10000]
-v number of singular vectors [5]
-V type of velocity field [0]
0=2-D global, azimuthal equidistant projection
1=2-D, date-based time grid
2=analytical
-w convert vertical velocities
-W output all (interpolated) tracer fields
-x number of x/longitude grids [360]
-X length of tracer field in x-direction [20000]
-y number of y/latitude grids [181]
-Y length of tracer field in y-direction [20000]
-z number of contours [20]
-Z run jobs in background
-? print help screen
SCRIPTS
There are several directories under scripts/. All of these contain makefiles
and shell scripts for higher-level processing of the compiled command-line
executables. Most of the shell scripts accept the same options as the
program they are calling or wrapping. To access the help screen, use -H
instead of -? since the latter is not supported by the default command option
processor for bash. Here is a breakdown of what's in the different diretories:
scripts/test/
Contains a makefile for running and plotting both an advected contour and a
passive tracer simulation. This can either be studied to pick up how to use
the codes or can simply be modified for use in your own research.
"make contour" generates an animation of an advected contour;
"make tracer" generates an animation of a passive tracer;
"make both" (the default) animates both in the same gif.
To produce the animations you will need both the Generic Mapping Tools (GMT)
and Imagemagick (convert).
scripts/plotting/
plot_frame.sh: plots a single frame from either a contour advection simulation
or a tracer simulation
animate_ctraj.sh: Animates an advected contour or a tracer simulation.
animate_both.sh: Animates a contour advection simulation AND a tracer
simulation on the same set of frames.
gen_zgrid: Generates a set of contour levels for plotting a tracer.
scripts/pc_proxy/
Tests the PC Proxy algorithm on both a simulated dataset and on actual ozone
measurements from the Polar Ozone and Aerosol Measurement (POAM) III instrument.
Compares results from the latter with a series of sonde measurements.
See the paper, "Principal component proxy tracer analysis." arxiv:1202.1999v1.
scripts/space_filling/
Fills up a lot of space on your hard-drive. Eats a lot of CPU cycles as well.
TECHNICAL NOTES
Most of what the codes do is integrate the differential equation:
dx/dt = v(x, t)
where x is position (a vector), t is unitless time (normalized to the time grid-spacing) and v is assumed to be known at all desired positions and is contained as a gridded dataset in the specified input file (note: analytical wind fields are now supported). In order to make this work for a curved Earth, both the gridding of the velocity fields and the velocity fields themselves must be suitably transformed. All codes designed for global simulations take two velocity fields as input: one for the Northern hemisphere and one for the Southern hemisphere, both of which must be gridded in an azimuthal equidistant coordinate system with the velocities suitably transformed.
The utility, nc2ds, converts NCEP data in NetCDF format to a format recognizable by the trajectory codes. Both pressure levels and potential temperature (theta) levels may be extracted. You will need files containing one year of either daily or 4-times daily wind data: uwnd.yyyy.nc and vwnd.yyyy.nc, where yyyy is the year and also temperature data, air.yyyy.nc, for converting theta levels. The utilities take as arguments the start and end dates, the desired level and the name of the two output files: one for the Southern hemisphere and one for the Southern. These are typically assigned a ".ds" extension. Optional arguments are the side-length and the number of grid points per side.
Since trajectories are integrated on one of two hemispheres and the hemisphere is corrected only once per large time-step, there must be some overlap between the two hemispheres. A side-length of 20000 will produce a grid whose boundaries only just reach the equator at the middle of the side edge, while extending well across the equator at the corners.
A note on the time-stepping: except for the trajectory model which works differently, output is done once per large time step. There can be any number of Runge-Kutta steps in between. In the contour advection simulation, points are "fixed" (both for the hemisphere and to patch the advected contour) once per large time step. In the tracer simulation, back-trajectories are performed and interpolation coefficients calculated for each large time step.
The tracer simulations output a series of sparse matrices that define a mapping from the initial tracer distribution to all subsequent time steps. There are two versions: one excludes all redundant points, while the other does not. The version that excludes the redundant points will output a mapping from each grid point to the position in the vector--values of -1 are excluded. Once the initial tracer state has been converted into a vector, it can be multiplied through with one of the utilities in the sparse matrix libary. Vectors are stored as simple binary files with a 4-byte header denoting their size--this is the same format as that used for the libagf library.
STEP-BY-STEP
The codes can be run by hand or for an automated version, see the Makefile under scripts/.
Step 1: Download wind data
The trajectory codes are driven by wind fields from GCMs or assimilation models. These data need to be converted to the native binary format. Currently, conversion utilities only exist to convert National Center for Environmental Prediction (NCEP) reanalyis data which can be obtained free of charge from the National Center for Atmospheric Research (NCAR) website. You will need daily or 4x daily air temperature data ("air.YYYY.nc"), zonal wind data ("uwind.YYYY.nc") and meridional wind data ("vwind.YYYY.nc") on pressure levels for every year for which you plan on integrating the models. Store all years in the same directory.
Step 2: Convert the wind data
Suppose you have stored the data in the directory, "/home/joe/data", and you want to integrate the models between November 13, 2001 and March 3, 2002 on a 700 K isentropic level. The following command will convert data for the Northern hemisphere and store it in the files, "vfield_S.ds" and "vfield_N.ds":
> nc2ds -T -p /home/joe/data -i 2001/11/13 -f 2002/3/3 700 vfield_S.ds vfield_N.ds
The "-T" option specifies isentropic (potential temperature) levels (default is pressure). These files can be used both for contour advection and two-dimensional tracer simulations. Multiple files on different levels will also be used for three-dimensional tracer simulations, however this has not been implemented yet.
- FOR 2-D PASSIVE TRACER MODELS
Step 3: Generate the tracer mapping
We want to integrate the velocity data created above over their entire date range (assuming 6hr NCEP data) with a resolution of 100 grid points per side per hemisphere and place the data in the file "mat.dat".
This file can be used for fast integration of any tracer field within the specified time span.
Note that the new programs only support command line options instead of an initialization file.
> ctraj_tracer -i 2001/11/13 -f 2002/3/3 -h 1 -k 4 -n 100 vfield_N.ds vfield_S.ds mat.dat
Step 4: Create the initial tracer field
The command "lla2aeb" transforms tracer fields in ASCII format, longitude-latitude coordinates to a binary vector in azimuthal equidistant coordinates. The command takes takes a list of each tracer value from STDIN. Input tracer values are arranged to cycle over longitude grids first, then latitude grids. Use the "-x" and "-y" options to specify the resolution (default is 360 x 181).
> lla2aeb -n 100 output.dat < tracer_in.txt
To start with a zonally symmetric tracer:
> zonally_symmetric_tracer | lla2aeb -n 100 output.dat
Note that the tracer resolution for the azimuthally equidistant coordinate system (the working coordinates) must match that of the tracer mapping. Here we have placed the initial tracer field in the file, "output.dat".
Step 5: Perform the matrix multiplication
Use the executable from the "libsparse" library which multiplies a vector with an array of sparse matrices:
> sparse_vect_mult mat.dat output.dat
The results are appended to the file "output.dat". You can also use the sparse matrix calculator, "sparse_calc".
Step 6: Extract the desired field
To extract the 100th field and print it to STDOUT:
> extract_field output.dat 100
To get the number of fields, omit the index:
> extract_field output.dat
For a list of dates, use "tracer4_matout" with the "-Q" option:
> tracer4_matout -Q tracer.init
- FOR 2-D CONTOUR ADVECTION
Step 3: Create the initial contour
To convert lon-lat pairs to the native binary format, use the command xy2bev:
> xy2bev contour_in.bev < contour_in.txt
For a single line of latitude, use the "C0" command:
> C0 45 | xy2bev
which generates a contour along 45 degrees North.
Step 4: Run the contour advection simulation
We need to specify the maximum fraction of arc traced out by each pair of points, as well as the minimum and maximum spacing using the -c, -m and -s options, respectively.
Once again, we can only use command line switches:
> cp contour_in.bev output.dat
> ctraj_contour -i 2001/11/13 -N 90 -h 1 -k 4 -c 1 -m 10 -s 100 vfield_N.ds vfield_S.ds output.dat
Because it lacks surgery, the ctraj_contour program is limited in how long it can advect the contours. To run longer contours, have a look at the "space_filling" directory under scripts--here you will find scripts that recursively split the contour to help it run longer.
Step 5: Extract the desired contour
To print out the 50th contour:
> bev2xy output.dat 50
To print out a list of records, omit the index:
> bev2xy output.dat
Step 6/7: Plotting
To plot the 50th field:
> plot_frame.sh output.dat 50
To animate the entire tracer integration:
> animate_ctraj.sh output.dat output.gif
Plotting scripts are currently contained in the "scripts" directory.
--------------
Chief developer: Peter Mills (petey@peteysoft.org)
Download Latest Version
libpetey_r289.tgz (141.5 kB)
