-----------------------------------------------------------------
   |    E3D: 2D/3D ELASTIC FINITE-DIFFERENCE WAVE PROPAGATION CODE   |
    -----------------------------------------------------------------

    =================================================================
   | (c) Copyright 1992 The Regents of the University of California. |
   |     All Rights Reserved.                                        |
   |                                                                 |
   |     This work was produced under the sponsorship of the         |
   |     United States Department of Energy.                         |
   |     The Government retains certain rights therein.              |
    =================================================================


               *******************************************
               * Shawn Larsen                            *
               * L-203                                   *
               * Lawrence Livermore National Laboratory  * 
               * Livermore, CA  94551                    * 
               * (925) 423-9617 (Office)                 * 
               * (925) 423-2163 (FAX)                    * 
               * larsen8@llnl.gov                        * 
               *******************************************



                           INTRODUCTION
                           ------------


E3D is an explicit 2D/3D elastic finite-difference wave propagation code 
used for the modeling of seismic waves. It is 4th-order accurate in space 
and 2nd-order accurate in time. It is based on the elastodynamic formulation 
of the wave equation on a staggered grid [Madariaga, 1976, BSSA; Virieux, 
1986, Geophysics; Levander, 1988, Geophysics; Larsen and Harris, 1993, UCRL]. 
The grid is staggered in both space and time. It is regularly spaced, except 
for the static grid refinement option. The computed variables at each node 
are the velocities and the components of the stress tensor.

Input consists of a single parameter file containing information about 
run-time options, grid structure, time stepping parameters, source functions, 
velocity model, and output options. Various data also can be input from 
separate files. Output consists of seismograms (SAC format), images (floating 
point), and run-time visualization.

Features include: Absorbing (non-reflecting) boundary conditions, sponge
boundary conditions, stress-free surface conditions, multiple sources, 
attenuation, topography (2D), 1D static grid refinement, hybridization,
parallelization (shared memory/message passing), run-time visualization, 
SAC output, image output, pure acoustic modeling option (for efficiency).

Major features to be included: 3D topography, 3D static grid refinement,
anisotropy. 

E3D has been implemented on workstations (e.g., Sun, SGI, HP's, DEC Alpha's), 
workstation networks (PVM), shared memory multiprocessors (Sun Enterprise,
SGI Origin), vector processors (Cray YMP, Cray C-90), and massively parallel 
processors using mpsc or mpi (Meiko CS-2, Cray T3D/T3E, nCUBE-2, IBM SP-2,
Compaq/DEC).



                          PROGRAM EXECUTION
                          -----------------

E3D takes an input file containing input parameters (which can include the 
velocity model), run-time options, and output options. To run type:

                          e3d Inputfile

To get a brief usage listing containing InputFile parameters simply
type "e3d" without any arguments. To enter this documentation from
the vi editor type "e3d -v". 



                            INPUT FILE
                            ----------

The input file contains a series of line commands. Each line command is 
composed of a line-type identifier followed by one or more attributes and
with attribute values. Each line (and related attributes) defines some
functional operation of the E3D code. 

The general format of a line command is:

 <line-type>  <attr1=attr1_value> <attr2=attr2_value> ... <attrn=attrn_value>

Any number of white-space characters can precede the line-type definition.
The line-type definition must precede all attributes and attribute values.
The order of the attribute values is not important. There must not
be any space between the attribute definition, the equals sign, and the
attribute value. If there are too many attributes to fit on a single line,
a "\" character as the (absolute) last character in the line designates
a wrap-around or continuation line. That is, the remaining attributes
can be listed on the immediately succeeding line(s).

Any line with an unrecognized line type will be ignored. Any unrecognized
attribute definition also will be ignored.  Any line with a "#" as the first
non-white space character automatically will be considered a comment line. 
[Unimportant Detail: The only difference between an ignored unrecognized line 
and a comment line is that if a line begins with a "#" as the first non-white
space character and is designated a continuation line with a "\" as the 
last character, then the following line(s) that are part of the continuation
line must begin with a "#" (otherwise, it may be processed as a legitimate 
line). Any unrecognized attribute definition will be ignored.]

For example, a format to define the dimensions of a 3D finite-difference
grid is:

  grid    x=50. y=40. z=20. dh=0.15 
                       
In this example, "grid" is the line-type identifier and x, y, z, and dh,
are attributes followed by an attribute value. In this case x, y, and z
are the grid dimensions in kilometers and dh represents the grid spacing 
in kilometers.

There are several line-types (with corresponding attributes) and these are 
discussed in greater detail below. Each line-type defines a different 
function. Some define the operation of the code, while others define
desired input and output options. A few line-types are required, most
are optional. Correspondingly, some of the attributes for any given line-type
are required and others are optional. In some cases, default values will
be used if a given parameter is not specified.

For some line-types, there should be only 1 reference to that line-type in
the input file. These are known as single reference line types. For other 
line types, there can be multiple references (usually with different 
argument values). These are known as multiple reference line-types.
For some multiple reference line types, such as the "block" line-type, 
a succeeding line-type can "overwrite" some or all of the values defined
by a preceding reference to that same line-type.

The order of the input file is sequential. That is, the 1st line of the
file is processed first, followed by the 2nd line, and so on. If a
single reference line-type is encountered more than once in the input file,
the last one encountered will be the one used by the E3D code.

There are a few options for inputing additional information into the code.
The multiple reference line-type "file" allows a separate file containing
additional line commands to be read and processed (in place) into the input
file. This can be used to make the input file well organization. For some
line-types, there are attributes which allow additional data or parameters
to be read.



                            LINE TYPES
                            ----------

The following line types are recognized by the E3D code.

grid:     defines the grid dimensions, parameters, and some run time options.
time:     defines the time stepping parameters.
source:   defines the source function(s).
visual:   defines run-time visualization parameters.
image:    defines snap shots (or movies) of output images or data files.
volume:   defines snap shots (or movies) of output 3-D volumes of data.
sac:      defines location and style of seismogram output (in SAC format).
traces:   defines location and style of seismogram output (in multi-traces).
block:    defines the grid medium (velocity) values.
poly3d:   defines a 3D polygon velocity file
vfile:    defines a grid-based velocity file 
topo:     defines topography parameters
restart:  defines simulation dump and restart parameters [not supported]
parallel: defines parallelization parameters 
file:     defines in-line reading of additional line commands.

The line-types "grid", "time", "visual", "topo", "restart", and "parallel" are
single reference line-types.  The others are multiple reference line-types. 
The "block" line-type can overwrite some or all of the information supplied 
by a previous "block" line-type.

The "grid" and "time" line-types are required. It would be all but pointless
to run the code without at least 1 "source" line-type. It would be all
but pointless to run the code without at least 1 of the output
line-types ("visual", "image", or "sac"). The velocity grid is set with
the "block" line-type. A default velocity model will be used if no
"block" line-type is specified.

In addition, the line-type "medium" is available but has been superseded
by line-types "block" and "vfile". It shouldn't be used except for the
input of 2D polygon objects.



                  INPUT PARAMETERS (GENERAL NOTES)
                  --------------------------------

The input units tend to be commonly used values in solid earth seismology.
Time is in seconds. Distances are in km. Seismic velocities (velocity model)
are in km/s. Density is in g/cm**3. Moment is in dyne-cm.

MKS units are used in the code.

Output units are in MKS. For example, ground velocities are in m/s.

A left-handed coordinate system is used. In 2D, the origin is in the
upper left corner of the grid. The x coordinate is positive to the right,
and the z coordinate is positive down. In 3D, the origin is in the upper
left front of the grid. The x coordinate is positive to the right, the
z coordinate is positive down, and the y coordinate is positive into the
plane.

Many parameters accept input coordinates as either x, y, z, or as l, m, n.
The x, y, z coordinates are the distances (km) from the origin. The l, m, n
coordinates are the associated grid nodes, starting at 0 (l == y, m == z, 
n == x).



                           STAGGERED GRID
                           --------------

A staggered grid is used for the finite-difference grid, where the simulated 
state variables (the velocities Vx, Vy, and Vz, and the stress tensor 
components Txx, Tyy, Tzz, Txy, Txz, Tyz) are spatially staggered from one 
another. In the E3D implementation, Vx is considered the reference point for 
each virtual node (the other variables are staggered by 1/2 grid point from 
this reference point).

For example, a given virtual node has the following configuration:


            Txy ---- Vy               y
           /        / |              /
          /        /  |             /
         Vx ---- Tnn  |            o-------> x
         |        |  Tyz           |
         |        | /              |
         |        |/               |
        Txz ---- Vz                z
                     

Tnn corresponds to the normal stresses (Txx, Tyy, Tzz). All other variables
are spaced 1/2 grid point from the others. The Back Face is spaced 1/2 grid 
point in the positive y direction from the Front Face. The variable Vx is
the reference. That is, if it's at grid point (l,m,n), then variable Tyz
is at location (l+1/2,m+1/2,n+1/2). For the purposes of discussion, however,
Tyz is said to be at grid point (l,m,n).

The finite-difference technique is also staggered in time, but this 
should not be a significant user issue.



                        LINE TYPE ATTRIBUTES
                        --------------------

The following represent the attributes available for each line-type.
The default value (when applicable) is enclosed in brackets [default].

grid   
----
 x:        the grid dimensions in the x direction [required]
 y:        the grid dimensions in the y direction [0. == 2D]
 z:        the grid dimensions in the z direction [required]
 l:        the number of nodes in the y direction [1  == 2D]
 m:        the number of nodes in the z direction [required] (also v grid)
 n:        the number of nodes in the x direction [required]
 dh:       the grid spacing in all directions
 dx:       the grid spacing in kilometers in the x direction (use dh instead)
 dy:       the grid spacing in kilometers in the y direction (use dh instead)
 dz:       the grid spacing in kilometers in the z direction (use dh instead)
 active:   implement active grid algorithm (0=off;1=h;2=h&v;3=red,4=box) [1]
 vred:     reducing velocity [maximum P velocity]
 tlen:     length of time window [0.75*(time of calculation)]
 damp:     number of grid points along boundary to apply damping [0]
 adamp:    damping factor along boundary [0.0]
 q:        flag to implement attenuation (0=off; 1=on) [0]
 model:    modeling mode (1=elastic; 2=pure_acoustic) [1]
 b:        grid boundary conditions [3]
           1=reflecting          
           2=absorbing
           3=surface
           4=surface (modified for acoustic instabilities)

time
----
 dt:       time-stepping increment [required]
 t:        number of time-steps [required]

source   
------
 type:     source type [required] 
           1=P  
           2=S  
           4=Moment Tensor  
           5=Force  
           6=Fault (point source) 
           7=Fault (finite length)  
           8=Distributed Slip (from file)
           9=Distributed Slip (from file with rise-times)
 amp:      amplitude of source (in dyne-cm) [required]
 t0:       the shift it "start" time of the source [0.0]
 freq:     central frequency of source in Hz [1.0]
 x:        x position of source (kilometers) [center of grid]
 y:        y position of source (kilometers) [center of grid]
 z:        z position of source (kilometers) [center of grid]
 l:        the y grid point of source [center of grid]
 m:        the z grid point of source [center of grid]
 n:        the x grid point of source [center of grid]
 strike:   strike of moment source (degrees, aki and richards)
 dip:      dip    of moment source (degrees, aki and richards)
 rake:     rake   of moment source (degrees, aki and richards)
 length:   length of fault (km, along strike) [undefined]
 width:    width  of fault (km, down dip)     [undefined]
 depth:    depth  of fault (km, from top)     [undefined]
 v:        rupture velocity (km/s)            [undefined]
 s0:       distance of hypocenter from fault center in strike direction [   0.]
 d0:       distance of hypocenter from fault fault top (down dip)   [undefined]
 Mxx:      the xx moment tensor component  (mxx also works) [0.0]
 Myy:      the yy moment tensor component  (myy also works) [0.0]
 Mzz:      the zz moment tensor component  (mzz also works) [0.0]
 Mxy:      the xy moment tensor component  (mxy also works) [0.0]
 Mxz:      the xz moment tensor component  (mxy also works) [0.0]
 Myz:      the yz moment tensor component  (myz also works) [0.0]
 Fx:       forcing function in x direction (fx  also works) [0.0]
 Fy:       forcing function in y direction (fy  also works) [0.0]
 Fz:       forcing function in z direction (fz  also works) [0.0]
 join:     flag for hybrid grid option
 xjoin:    offset parameter for hybrid grid option
 tjoin:    number of time steps for hybrid grid option
 file:     name of input file containing sac file of source [none]
 dfile:    name of distributed slip input file [none]
 p:        number of points in distributed slip file [undefined]

block
-----
 p:        P wave velocity (P and alpha also work) (km/sec)  [none]
 s:        S wave velocity (S and beta  also work) (km/sec)  [none]
 r:        density (R and rho also work)           (g/cm**3) [none]
 ps:       p/s ratio                                         [none]
 gradient: vertical gradient (units per km)                  [  0.]
 Q:        Q value for P and S wave attenuation              [none]
 Qp:       Q value for P wave attenuation                    [none]
 Qs:       Q value for S wave attenuation                    [none]
 Qf:       central frequency of attenuation curve            [none]
 x1:       starting x position of block element       [         0.]
 x2:       ending   x position of block element       [end of grid]
 x :       sets x1 and x2 to same value
 y1:       starting y position of block element       [         0.]
 y2:       ending   y position of block element       [end of grid]
 y :       sets y1 and y2 to same value
 z1:       starting z position of block element       [         0.]
 z2:       ending   z position of block element       [end of grid]
 z :       sets z1 and z2 to same value
 l1:       starting y grid point of block element     [          0]
 l2:       ending   y grid point of block element     [end of grid]
 l:        sets both l1 and l2 to same value
 m1:       starting z grid point of block element     [          0]
 m2:       ending   z grid point of block element     [end of grid]
 m:        sets both m1 and m2 to same value
 n1:       starting x grid point of block element     [          0]
 n2:       ending   x grid point of block element     [end of grid]
 n:        sets both n1 and n2 to same value
 
poly3d
------
 p:        P wave velocity (P and alpha also work) (km/sec)  [none]
 s:        S wave velocity (S and beta  also work) (km/sec)  [none]
 r:        density (R and rho also work)           (g/cm**3) [none]
 ps:       p/s ratio                                         [none]
 gradient: vertical gradient (units per km)                  [  0.]
 Q:        Q value for P and S wave attenuation              [none]
 Qp:       Q value for P wave attenuation                    [none]
 Qs:       Q value for S wave attenuation                    [none]
 Qf:       central frequency of attenuation curve            [none]
 x1:       starting x position of block element       [         0.]
 x2:       ending   x position of block element       [end of grid]
 x :       sets x1 and x2 to same value
 y1:       starting y position of block element       [         0.]
 y2:       ending   y position of block element       [end of grid]
 y :       sets y1 and y2 to same value
 z1:       starting z position of block element       [         0.]
 z2:       ending   z position of block element       [end of grid]
 z :       sets z1 and z2 to same value
 l1:       starting y grid point of block element               [0]
 l2:       ending   y grid point of block element     [end of grid]
 l:        sets both l1 and l2 to same value
 m1:       starting z grid point of block element               [0]
 m2:       ending   z grid point of block element     [end of grid]
 m:        sets both m1 and m2 to same value
 n1:       starting x grid point of block element               [0]
 n2:       ending   x grid point of block element     [end of grid]
 n:        sets both n1 and n2 to same value
 lon:      grid origin longitude (longitude also works)  [required]
 lat:      grid origin latitude  (latitude also work)    [required]
 angle:    rotation angle of grid (ccw from east)              [0.]
 file:     name of input 3D polygon file                 [required]

vfile
-----
 p:        P wave velocity (P and alpha also work) (km/sec)  [none]
 s:        S wave velocity (S and beta  also work) (km/sec)  [none]
 r:        density (R and rho also work)           (g/cm**3) [none]
 ps:       p/s ratio                                         [none]
 gradient: vertical gradient (units per km)                  [  0.]
 Q:        Q value for P and S wave attenuation              [none]
 Qp:       Q value for P wave attenuation                    [none]
 Qs:       Q value for S wave attenuation                    [none]
 Qf:       central frequency of attenuation curve            [none]
 l1:       starting y grid point of block element        [required]
 l2:       ending   y grid point of block element        [required]
 l:        sets both l1 and l2 to same value
 m1:       starting z grid point of block element        [required]   
 m2:       ending   z grid point of block element        [required]
 m:        sets both m1 and m2 to same value
 n1:       starting x grid point of block element        [required]
 n2:       ending   x grid point of block element        [required]
 n:        sets both n1 and n2 to same value
 type:     input file type (p,P,s,S,r,R,Q,Qp,Qs,Qf)     [undefined]
 file:     name of input velocity file                   [required]
 c:        rle file compression flag (0=no, 1=yes)              [0] 
 seg:      undocumented seg filling flag (0=no, 1=yes)          [0]
 
medium  (this option has been superseded by "block", "vfile", and "poly3d")
------
 alpha:    P wave velocity in km/sec [6.0]
 beta:     S wave velocity in km/sec [alpha/PSratio]
 rho:      density in g/cm**3        [2.52 + .3788*alpha]
 PSratio:  P to S ratio [sqrt(3.0)]
 gradient: vertical gradient for P and S velocities [0.0]
 Q:        Q value for P and S wave attenuation [10,000]
 Qp:       Q value for P-wave attenuation [10,000]
 Qs:       Q value for S-wave attenuation [10,000]
 Qf:       central frequency of attenuation curve [1.0 or first source]
 x1:       starting x position of medium element [0.]
 x2:       ending   x position of medium element [end of grid]
 x :       sets x1 and x2 to same value
 y1:       starting x position of medium element [0.]
 y2:       ending   x position of medium element [end of grid]
 y :       sets x1 and x2 to same value
 z1:       starting x position of medium element [0.]
 z2:       ending   x position of medium element [end of grid]
 z :       sets x1 and x2 to same value
 l1:       starting y grid point of medium element [0]
 l2:       ending   y grid point of velocity element [end of grid]
 l:        sets both l1 and l2 to same value
 m1:       starting z grid point of velocity element [0]
 m2:       ending   z grid point of velocity element [end of grid]
 m:        sets both m1 and m2 to same value
 n1:       starting x grid point of velocity element [0]
 n2:       ending   x grid point of velocity element [end of grid]
 n:        sets both n1 and n2 to same value
 type:     type of velocity element (1=block; 3=polygon; 10=vfile) [1 or p]
 p:        number of points in polygon (which follow) (resets type)
 mode:     velocity file format (1=asciiP,2=binaryP,3=asciiALL,4=binaryALL) [1]
 file:     velocity file [none]

topo   
----
 file:     name of topography file (in sac format) [none]
 depth:    depth of flat topographic interface [0.]

visual
------
 movie:    time-step increment for run-time visualization [0]
 scale:    scale factor for run-time visualization image [1.0]
 sample:   sampling increment for run-time visualization image [1]
 mode:     mode of run-time visualization (1=normal) [1]
 x,y,z:    3D coordinate of visual plane  [y=source or middle of grid]
 n,l,m:    3D node number of visual plane [l=source or middle of grid]

image   
-----
 t:        time-step to output image [0]
 movie:    time-step increment to output images [0]
 file:     file name header of image ["image"]
 scale:    scale factor of image [1.0]
 sample:   sampling factor of image [1]
 mode:     mode of image [1]
           1  P,S 
           2  X,Y,Z 
           3  8 bit Sun Raster (don't use)
           4  maximum horizontal velocity
           5  energy array
           6  p velocity array
           7  s velocity array
           8  density array
           11 Vx  - x velocity  over grid (single file)
           12 Vy  - y velocity  over grid (single file)
           13 Vz  - z velocity  over grid (single file)
           14 Txx - Txx stress  over grid (single file)
           15 Tyy - Tyy stress  over grid (single file)
           16 Tzz - Tzz stress  over grid (single file)
           17 Txy - Txy stress  over grid (single file)
           18 Txz - Txz stress  over grid (single file)
           19 Tyz - Tyz stress  over grid (single file)
           20 P   - P potential over grid (single file)
           21 S   - S potential over grid (single file)
 x,y,z:    3D coordinate of image plane  [y=source or middle of grid] 
 n,l,m:    3D node number of image plane [l=source or middle of grid]

volume
------
 t:        time-step to output 3-D data volume [0]
 movie:    time-step increment to output 3-D data volume [0]
 file:     file name header of volume file ["volume"]
 scale:    scale factor of image [1.0] (doesn't work)
 sample:   sampling factor of image [1]
 mode:     mode of image [1]
           1  P potentials
           2  S potentials

sac   
---
 x:        x position of seismogram (kilometers) 
 y:        y position of seismogram (kilometers)
 z:        z position of seismogram (kilometers) [0.]
 l:        y grid point of SAC seismogram 
 m:        z grid point of SAC seismogram [0]
 n:        x grid point of SAC seismogram
 scale:    scale factor for seismogram [1.0] (doesn't work)
 sample:   sample factor for seismogram [1]
 mode:     type of seismogram output (1=regular; 2=moment) [1]
 file:     file name header of seismograms ["sac"]

traces
------
 mode:     bitwise component(s) to output (1=Vx;2=Vy;4=Vz;8=P;16=S) [31 (all)]
 sample:   time sampling factor for traces [1]
 scale:    scale factor for traces [1.0] (doesn't work)
 file:     file name header of traces file ["traces"]
 tfile:    input file containing trace coordinates [none]
 x1:       starting x coordinate for traces (km) [0.]
 x2:       ending   x coordinate for traces (km) [0.]
 y1:       starting y coordinate for traces (km) [0.]
 y2:       ending   y coordinate for traces (km) [0.]
 z1:       starting z coordinate for traces (km) [0.]
 z2:       ending   z coordinate for traces (km) [0.]
  p:       number of traces along coordinate line (inclusive) [0]

restart
-------
 restart:  time-step to find dump files and restart simulation
 dump:     time-step increment to dump simulation data for restart
 nolast:   flag to delete last (and all) dump/restart files at final timestep
 file:     filename header for dump and restart files

parallel    
--------
 nx:       number of parallel nodes in x direction [1]
 ny:       number of parallel nodes in y direction [1]
 nz:       number of parallel nodes in z direction [1]
 file:     input file for optional mpp grid decomposition [none]

file   
----
 file: file containing additional command lines for inline insertion [none]



                         EXAMPLE INPUT FILES
                         -------------------
Example 1: 
---------
   # This is a simple example input file. To run type "e3d inputFile".
   #
     grid    m=100 n=200 dh=0.15
     time    t=200 dt=0.045
     block   p=2.00 s=1.55 r=3.0
     source  type=1 amp=5.e23 freq=0.5 m=50 n=100
     visual  movie=10

This input file describes a finite-difference grid that has 200 grid points
in the horizontal direction and 100 grid points in the vertical direction.
The grid spacing is 0.15 km. Hence, the grid represents a region 30 km long 
and 15 km deep. The grid is considered 2D because "l" (or "y") is not 
indicated.  The time step is 0.045 s, which just satisfies the Courant 
condition. The simulation lasts for 200 time steps, or 9 seconds.  In this 
case, a homogeneous velocity model is used with a compressional velocity 
of 6.0 km/s, a shear wave velocity of 4.0 km/s, and a density of 
3.0 gm/cm**3. A compressional point source with a Ricker wavelet time 
history is used for the source term. The amplitude isn't really defined, 
since this is a 2D problem (in fact, the source mimics a 3D line source). 
The frequency is 0.5 hz, and the source is located in the middle of the 
grid.  There is no output, except for the run-time visualization which is
update every 10 time steps. Note that there are over 20 grid points per
wavelength (at the lowest non-zero grid velocity of 1.55 km/s), so this
grid is oversampled by about 2 times (10 grid points per wavelength is
a good rule of thumb).

Example 2:
---------
   # This is a slightly more complicated input file describing additional
   # input methodology. Note that unrecognized lines or lines beginning with
   # a "#" are treated as comments. Also note that "wrap around" lines are
   # allowed with a "\" symbol before the newline character, and that blank
   # lines are allowed. Furthermore, unrecognized lines and attributes will
   # be ignored.

   # We use the same sized grid as Example 1 but now we use km instead of 
   # grid nodes. The line "ggrid" is not recognized so it is ignored. The
   # attributes "mm" and "nn" are not recognized so they are ignored.
   # Since active=1 is the default, it can be included or ignored.
     ggrid   m=100 n=200 dh=0.15
     grid    x=30. z=15. dh=0.15 mm=100 nn=200 active=1

     time    t=200 dt=0.045

   # This velocity model represents a layered model, with a "box" of water
   # in the center right (obviously not a realistic model).
   # Note how succeeding "block" line-types overwrite previous information
   # Also note how the "block" attributes (e.g., p, s, or r) can be input 
   # on the same line or on different lines. Note "wrap around" line for 
   # last "block" statement.

     # sets entire grid
     block   p=2.00 s=1.55 r=3.0           

     # sets top layer as a gradient (goes from 0. to 5. km)
     block   p=1.00 gradient=0.20  z2=5.  
     block   s=0.78 gradient=0.154 z2=5. 

     # puts in "water box"
     block   p=1.50 s=0.00 r=1.0 \
             x1=20. x2=23. m1=40 m2=60     

   # use a moment source with the off-diagonal element set
     source  type=4 amp=1.e27 freq=0.5 m=50 x=15. mxz=1.0

   # Output (run time visualization), trace files (sac format), and image 
   # files. The trace files and image files are commented out.
     visual  movie=10
   # sac     x=10. z=2. file="sac1"
   # sac     x=20. z=2. file="sac1"
   # image   t=100 mode=1



                          PARAMETER USAGE
                          ---------------

The following describes some of the features and the specific methods 
to handle various E3D options. 

Grid Dimensions
---------------
The grid dimensions are specified using the attributes l, m, n, x, y, z, and 
dh for the "grid" input line. The parameter dh specifies the grid spacing
(in km), the parameters l, m, and n specify the number of grid points in
the grid, and the parameters x, y, and z specify the physical dimensions
of the grid (in km). Generally, dh is input along with l or y, m or z, 
and n or x. You do not need to specify n and x, for example (it's possible
to specify the number of grid points and the physical dimension and not
specify dh, but this isn't recommended).

The coordinate system is left-handed. The x and n coordinates are in the
horizontal direction, positive to the right. The y and l coordinates are
positive into the plane of the board. The z and m coordinates are positive
down (depth).

Here are two examples of grid line-types, each for the same grid: 

  grid x=20. y=10. z=10  dh=0.10
  grid n=201 l=101 m=101 dh=0.10

In this case the grid is 20 km long (a total of 201 grid points), 10 km wide 
(101 grid points), and 10 km deep (101 grid points). The grid spacing is
0.10 km.

For 2-D problems, y/l is not specified.

The rule of thumb is that 10 grid points are needed for the shortest usable
wavelength in the model. For example, to simulate 1 hz seismic energy in
a medium with a low (non-zero) shear wave velocity of 1.0 km/s, the grid
spacing should be 0.1 km or less.

The default finite-difference grid is uniformly spaced. Static grid 
refinement is possible, where grid spacing is mapped onto the velocity
(non-uniformly spaced grid). However, static grid refinement is not
discussed here.

Time-Stepping
-------------
The "time" line-type defines the finite-difference timing parameters. There
are 2 attributes - "t" and "dt". The t attribute defines the number of
time-steps and the dt attribute defines the time-step increment. The
total simulated time will be "dt*t". 

The time-step increment must satisfy the Courant condition. This means that
  
  dt < factor*dh/Vmax

where dh is the grid spacing and Vmax is the maximum velocity in the grid 
(the P velocity). The "factor" constant is 0.606 for 2-D problems and 0.494 
for 3-D problems. If dt is greater than this relationship, then the solution 
quickly becomes unstable. For most applications, it's best to set dt as 
close to the Courant limit as possible. This improves the run-time, and it 
improves the accuracy. A 1% buffer is reasonable.

The 0.606 (2-D) and 0.494 (3-D) factors are calculated from the relationship
factor = 1./(sqrt(d)*(c1 + c2)) where d is the dimension (2 or 3), c1 = 9/8,
and c2 = 1/24 (c1 and c2 are the factors used in the 4th order differencing
stencils.

Model
-----
There are two types of physical models: elastic (and viscoelastic) and
acoustic. The model is specified using the "model" attribute on the "grid"
input line. The elastic model (model=1) is the default. The acoustic model
(model=2) is a formulation that includes the P velocity and the density. 
Note that acoustic simulations are still possible with the elastic model 
by setting the appropriate shear velocities to zero. The difference is that 
the acoustic model is computationally pure acoustic, which is much more 
efficient in both speed and memory. Typically, the acoustic model cuts 
memory usage in half and is 2.5 to 3.0 times faster. If b=3 is specified
as the boundary conditions under the acoustic model, the boundary conditions
will operate like they are b=4.

The acoustic model is more limiting, however. Currently, it is available only 
for 3-D simulations and attenuation cannot be used. Certain output options 
are not available such as mode=2 with sac output and mode=6,7,8,17,18,19 with
image output. The only available sources are types=1,5. The input model must
be specified with "block", "vfile", or "poly3d" input (shear wave velocity
and ps ratios specified with these inputs will be ignored). The "join" and
"restart" options are not available. 

Again, note that the acoustic model includes both P velocity and density
(it is not P velocity only), and acoustic simulations are still possible
with the elastic model by setting the S velocity to zero (although this
is computationally more expensive).

Propagating Grid
----------------
Propagating grid mechanisms are available for improving the run-time 
performance of E3D. These are initiated by the "active" attribute for
the "grid" line-type. If "active=0" then none of the propagating grid 
mechanisms are used.  This is useful for benchmarking. The effectiveness
of the propagating grid options is limited for problems run on massively
parallel processors.

If "active=1" (the default), then those regions of the finite-difference 
grid void of non-numeric seismic energy are not involved in the calculation. 
This occurs at the early times of a calculation. For typical problems, the 
run-time can be improved by a factor of 2. 

If "active=2", regions of the grid void of seismic energy are not involved 
in the calculation (like "active=1"). In addition, regions of the grid 
that do not participate in surface output are not computed. This occurs 
during the final stages of a calculation. This option modestly improves
performance over "active=1". 

If "active=3", a specified region of the grid is active in the calculation. 
This region moves with selected seismic phases of interest. The speed and
size of this moving region is defined by the "vred" (reducing velocity) 
and "tlen" (time length) attributes of the "grid" line-type. The front
of the region will travel at a speed corresponding to the fastest velocity
in the grid. The back of the region is defined by vred and tlen.  The 
reducing velocity must be greater or equal (default) to the maximum P wave 
velocity in the grid.  At present, this option can be effectively used only 
when the source is at the left side of the grid (low x value) and the 
receivers are on the right side of the grid (high x value). It's most useful 
in 2-D applications, such as those typically done for regional seismic 
studies. For some problems, "active=3" can improve performance by a factor 
of 5 (or more). 

If "active=4", a user time-variable specified region of the grid is active
in the calculation. This region (or box) can be set by the user to move
with desired phases. The moving box is defined and read from an input file
given by the "afile" attribute (specified on the "grid" line). Each line 
of the input file indicates a timestep, followed by minimum and maximum 
l,m,n (y,z,x) grid point coordinates that define the box at that timestep 
(and perhaps future timesteps). The format of each line of the input 
file is:

timestep lmin lmax mmin mmax nmin nmax

where timestep is the timestep number for that line (the first timestep
is 1, not 0). The minimum and maximum values define the activated 
grid for that timestep. They are defined as global grid point values
(l,m,n corresponds to the y,z,x dimensions, respectively). Note that
the (l,m,n) minimum and maximum values are inclusive (starting at 0
and ending at the number of grid points minus 1). If the minimum 
value falls outside the minimum value of the grid (i.e., < 0) then
it will be set to 0. If the maximum value falls outside the maximum 
value of the grid then it will be set to the maximum grid value. 
Lines with timestep values less than 1 or greater than the number of
timesteps in the simulation will be ignored. Note that the input file
cannot contain blank lines or comment statements. Each line must be a 
valid representation (unless the timestep numbers are out of bounds). 
Note that all timesteps do not need to be specified in the input file. 
If a timestep is not specified then the minimum and maximum coordinates 
at the previous timestep will be used (in other words, the most recent
timestep specified in the input file). Timesteps in the input file must 
be in increasing order (otherwise the behavior is not specified). If 
coordinate values for the first timestep are not specified, then 
coordinate values for the first specified timestep value will be used 
for the initial timesteps. The "active=4" option can be used with 2D or 
3D simulations. Note that with 2D simulations the minimum and maximum l (y) 
value is still required in the input file, although the values can be set 
to any value (to be safe, both lmin and lmax should be set to 0).

As an example of the "active=4" option, the e3d input file might 
contain the lines:

grid l=100 m=100 n=100 dh=1.0 active=4 afile="active.box"
time dt=0.100 t=1000

and the file "active.box" might contain the lines (remember no blank or 
comment lines):

  1   50  55    20  25     0   99    
  2   50  60    20  25     0   99    
 10   50  60    20  30     0   99    
100   50  70     0  50     0   99    
500    0  99     0  99     0   99    

Note that this is only a hypothetical example. The initial activated 
region of the grid goes from l=50 to 55 (inclusive), m=20 to 25, and
n=0 to 99 (the whole grid in the x dimension). For the 2nd timestep 
the grid increases in size in the l (y) dimension extending from grid
points 50 to 60. Since timesteps 3-9 aren't specified, the activated 
region will be the same as timestep 2. The activated grid region
changes in size (or moves) again at timesteps 10, 100, and 500. The
activated region stays the same from timesteps 500 to 1000. Note that
the input file is free formatted.

Note that reflections can occur at box boundaries.

Velocity Model Input
--------------------
There are several ways to input a velocity model into the grid. This is
done by one or more combinations of the "block", "vfile", and "poly3d"
line-types. These line-types are sequential, hence those that occur later
in the input file will over write previous information. This is useful
for building velocity models (e.g., one can lay a more complicated model
covering part of the grid over a base model).

In general, the velocity information for each line can span all or a 
portion of the grid. When it covers the entire grid, then coordinate 
information is not used. That is, the attributes "l,l1,l2 ... z,z1,z2" are 
not included on the appropriate input line. If these coordinates are 
specified, then only a fraction of the grid will be affected.

For example, rectangular regions corresponding to all or a specified portion
of the grid can be input with the "block" line-type. Typical block line-types 
look like:

  block p=6.0 s=4.0 r=3.0
  block p=5.0              n1=10 n2=20 z1=5.

The first line will fill the entire grid with a P velocity of 6.0 km/s,
an S velocity of 4.0 km/s, and a density of 3.0 gm/cm**3. The second 
line will over-write that portion of the grid between grid points 10
and 20 (inclusive) in the horizontal direction, from 5 km depth to the
bottom of the grid (since z2 is not specified) in the vertical direction, 
and the entire y direction (since l,l1,l2,y,y1,y2 are not specified) with a 
compressional velocity of 5.0 km/sec. The S wave velocity and the shear wave 
velocity will be unaffected in this region since they were not specified
on the second input line.

Other medium parameters that could be used are ps (for the p to s ratio),
gradient (for the gradient of any specified parameter in units per km),
or attenuation  coefficients (Q, Qs, Qp, Qf). The gradient only works in
the vertical direction. For example, the following
block line-type:

  block p=6.0 gradient=0.5 ps=1.7 z1=20. z2=24.

will set the p velocity to 6.0 km/s at a depth of 20 km, and it will be
linearly increasing to 8.0 km/s at a depth of 24 km. The S wave velocity
in this region will be the P wave velocity divided by 1.7. The density
and attenuation coefficients will not be set or changed in this example.

A simple geologic model, such as layered models, can be built using these 
"block" line-types.

A gridded velocity model can be read from a separate file (binary) using the 
"vfile" line-type. A grid representing one medium parameter (such as p, s, r, 
or attenuation coefficients) is included in a single vfile. However, other
parameters can be set or manipulated with the "vfile" option. For example,
a vfile input line will look something like:

  vfile type=p file="vel.p" ps=1.73 r=3.0 n1=0 n2=99 m1=0 m2=99 l1=0 l2=99

This line will read a gridded file (named "vel.p"), which contains P seismic 
velocities (defined by the "type=p" attribute), from a region of the grid
from grid points 0 to 99 in the x direction (inclusive), grid points 0 to 99 
in the y direction, and grid points 0 to 99 in the z direction. The S wave 
velocity will be set to "P/1.73" in this region, and the density will be set 
to 3.0 gm/cm**3. The "ps" and "r" attributes are optional, but in this case
they're used to avoid reading another gridded velocity file representing 
their values. 

Possible types are p or P for P wave velocity, s or S for S wave velocity, 
r or R for density, Q for attenuation (both Qp and Qs), Qp for P wave
attenuation, Qs for S wave attenuation, and Qf for the attenuation central
frequency range. If the "type" attribute is not include, then the type
can be determine from the file extension (assuming the extension is one
of the types).

There is and must be a one-to-one mapping of the points in the velocity
file and the finite-difference grid nodes within the specified region.
That is, the number of velocity values in the gridded velocity file must
equal "(l2-l1+1)*(m2-m1+1)*(n2-n1+1)".

However, the reverse is not true. That is, the velocity file does not have
to span all (or even any) of the finite-difference grid. The coordinate 
variables n1, n2 (and m1, m2, l1, l2) allow the specification of where the 
velocity file maps to the finite-difference grid. The velocity file can
be a subset or superset of the grid.

For example, consider the x dimension with the n index. Then "n2-n1+1" will
be the dimensions of the model specified in the velocity file. The index 
"n1" defines the starting grid point in the finite-difference grid and the 
index "n2" defines the ending grid point. Suppose the finite-difference
grid is 100 grid points long. If the velocity file exactly maps to the 
finite-difference grid then "n1=0" and "n2=99". This will be the normal case.
If the velocity file represents a fraction of the finite-difference grid then
the values of "n" might be something like "n1=20" and "n2=50". If the velocity
file represents a region larger than the finite-difference grid then
the values of n might be "n1=-50" and "n2=150".

Hence, a large velocity file can be used in a smaller finite-difference grid,
or a small velocity file can be used in a section of a larger finite-difference
grid. Therefore, a complicated region of the grid, such as a basin, might be 
specified by a gridded velocity file, with other methods (such as "block" 
input) being used to define the simpler regions of the finite-difference grid.

It is possible to use a "rle" compressed version of the binary gridded 
velocity file. This can save considerable disk space since many times there 
are regions in the model with the same value. The "c" attribute of the "vfile" 
line-type is used to specify whether the gridded velocity file is compressed 
or not. If the file is rle compressed, set "c=1"; if it's not, set "c=0" 
(which is the default).  A gridded velocity file can be compressed with the 
auxiliary routine "getrle". A compressed gridded velocity file can be 
uncompressed with the auxiliary routine "getunrle".

The format of the (uncompressed) binary gridded velocity files is (4 bytes per 
value):

  x index varies fastest
  z index varies next fastest
  y index varies last

Thus, if there is an array of values V[i][j][k], where i corresponds to the
y index, z corresponds to the j index, and k corresponds to the x index, then
the following loop structure will be used to write the binary gridded
velocity file (C notation):

  for(i = 0; i < ymax; i++)
    for(j = 0; j < zmax; j++)
      for(k = 0; k < xmax; k++)
        fwrite(&(V[i][j][k]), 4, 1, stdout);

Again, if there are large regions in the grid that are homogeneous then the
gridded file can be rle compressed, and this rle compressed file can be
utilized by E3D. 

The vfile option has an undocumented special attribute for SEG filling. 
If seg=1, then the p velocities already in the gridded region specified
on the vfile line (which presumably were just read into the grid with 
the file specified on the vfile line) will be used to create s and density 
values as specified by the relationships suggested by Leon Thompson and 
Walter Kessinger (see House et al. abstract for 2000 SEG meeting).

The "poly3d" line-type allows one to input a velocity model via 3-D polygon
solids. The input file is complex and non-documented (the method and algorithm
comes from the University of California at Berkeley). Hence, documentation on
using this option currently is not provided.

The "medium" line-type is an old and unsupported method, which has been
superseded by the "block", "vfile", and "poly3d" line-types. The only 
exception is the polygon input method for 2-D simulations. In the future,
this method will be redesignated to a "poly2d" line-type, but for now
it is briefly described here. For example, to input a simple 2-D polygon
shape, sample input lines are:

  medium type=3 p=5 alpha=6.0 beta=4.0 rho=3.0 gradient=1.0
    2.0 2.0
    3.0 2.0
    3.0 4.5
    3.3 3.3
    2.0 2.0

This would input a 5-corner polygon with vertices in kilometers at the
points indicated, and with a P wave velocity of 6.0 km/s, an S wave velocity
of 4.0 km/s, and a density of 3.0 gm/cm**3 at the top of the solid. A 
gradient of 1.0 km/s will be applied to the P wave velocity in this solid
(the S wave velocity will be adjusted accordingly). These points could
be specified in a separate file using the "file" line-type (following
the "medium" line-type declaration). The coordinates of the polygon corners 
are in km.

Keep in mind that p, s, and r must be specified someplace in the grid
(with any input line type). There can be no null values. Also note that
all of the input types ("block", "vfile", "poly3d", medium) can be used
together to easily build complicated models. For example, "block" input
might be used to describe a simple geology, and "vfile" input used to 
describe only those regions where the geology is complicated (such as
near surface basins).

Note that with the acoustic model (model=2 on the grid line), the only
input types are "block", "vfile", and "poly3d". Shear wave and ps ratio's 
specified with these options will be ignored.

Source
------
Several source types can be used. These include point sources, and sources 
representing finite-faults. The parameters that describe the source are
it's type, location, orientation, amplitude, frequency, and start-time.
Multiple sources can be included into the same run.

The user has the option of reading the source time-function from a predefined
file (e.g., SAC file) using the "file" attribute. The number of points in
the file must equal the number of time-steps in the finite-difference 
simulation (no longer true). If the "file" attribute is not used, then a Ricker
wavelet is used as the source time-function. The frequency of the Ricker 
wavelet (cycles/sec or Hz) will be defined with the "freq" attribute. The 
default frequency is 1 Hz. The Ricker peak (trough) will be offset from 0 by
"1./freq + t0" (in seconds), where t0 is a user defined attribute with a 
default of 0 seconds. The Ricker forcing function, which drives the stresses
or velocities is given by:

  factor = (PI*freq*(t - 1./freq - t0))**2
  f = amp*(2*factor - 1)*exp(-factor)

where "amp" is the amplitude attribute required on the source line. If the 
file attribute is used, and the source time history is defined in a separate
file (SAC file), then the "t0" and "freq" attributes will be ignored. However,
the input source time-history will be multiplied by the "amp" attribute.

The source amplitude is defined by the "amp" attribute of the source 
line-type. It is required, although it can equal 0. The amplitude is 
usually thought of in terms of moment (dyne-cm). Since a point source
in 2-D is really a line-source in 3-D, the absolute source amplitude value
doesn't have true physical significance for 2-D modeling problems.
Note that the source drives the stresses (or velocities) in the solution.
Hence, it is a velocity source.

The source location is defined by the "x, y, z" (km) or "l, m, n" (grid points)
attributes (or a mixture of the two). In either case, the location is not 
interpolated but is centered around the Tnn (normal stress variable) of
a given grid point. Hence, if "l=10 m=10 n=10" the source will be centered
around the Tnn variable at grid point "l=10 m=10 n=10". Since the Tnn variable
is shifted 1/2 grid point in the x direction from the reference point (Vx
variable), the source will be shifted by 1/2 grid point in the x direction.
When "x, y, z" are used to define the source location, the source is centered
around the nearest Tnn grid point. In the future it will be possible to 
specify the source at an exact position in space (through interpolation).
If a fault coordinate is not given, a default coordinate will be used
(center of the grid).

The type of source is described with the "type" attribute. The simplest type
("type=1") is a compressional or P source (explosive). A pure shear source 
can be described by "type=2", but for reasons that won't be explained here,
the derivative of the time-history will be used. 

When "type=4", the source is described in terms of the moment tensor 
components using the "Mxx, Myy, Mzz, Mxy, Mxz, Myz" attributes (all lower 
case works too).  Each component will be multiplied by "amp" value; any 
components not listed will have 0. value. For example "mxx=1. myy=1. mzz=1." 
will give a compressional source, just like type=1; "mxy=1." will give a 
shear source. 

When "type=5" the source is input as forcing functions that drive the
velocities in any given direction. The attributes "Fx", "Fy", and "Fz" 
define the forcing functions (which are multiplied by the "amp" value.
In this case, the amplitudes can be thought of in terms of force. The
units are in 0.01 dynes. That is, an amplitude of 1.0 equals 0.01 dynes.
Note that the source time function for a simple explosion source should
be described as a ramp function, not a Gaussian pulse. This takes some
thinking with a rubber band to understand why this is the case.

When "type=6" the source is described in terms of fault parameters.
The fault parameters are "strike", "dip", and "rake", with aki and richards 
convention being used (p. 117) (except that the X and Y axes are interchanged.
In this convention, x is east, y is north, and z is down. Strike is clockwise
from north (0 <= strike < 360); dip is measured down from the right of the 
strike (0 <= dip <= 90.); rake is the angle between the slip and strike
directions, where the slip is the direction of the hanging wall relative to
the foot wall (a positive rake is a reverse fault, a negative rake is a 
normal fault). A strike of 0., a dip of 90., and a rake of 0. corresponds
to a north-south left-lateral fault or an east-west right-lateral fault.
The fault will be a point source. The strike, dip, and rake (all of which
must be specified) will be used to calculate moment tensor components.

When "type=7" the source is described in terms of a finite-length fault
with uniform moment. The length of the fault is defined by the "length" 
attribute, and the "width" of the fault is defined by the "width" attribute 
(both in km). The depth to the top of the fault plane is defined by the
"depth" attribute (in km). The fault coordinates ("x, y, l, or n") define
the horizontal position of the center-top of the fault plane. Note" the
"z" (or "m") attribute will define the depth to the top of the fault plane
if the "depth" attribute is not given. The orientation of the plane is 
defined by the "strike" and "dip" attributes. The slip direction is defined
by the "rake" attribute. The hypocenter is defined by "s0" and "d0" (in km); 
s0 is the distance from the fault center in the positive strike direction,
and d0 is the distance from the fault top in the down dip direction. If s0
extends beyond the fault plane (i.e., -length/2 > s0 > length/2) then it
will be set to the nearest fault end. If d0 is negative or if it is greater
than the fault width, it will be set to the nearest fault end. The rupture
velocity is defined by the "v" attribute (km/s). The rupture will propagate
in a circular pattern outward from the hypocenter. The fault will pass
through (or near) grid nodes, and these nodes will act as independent 
point sources. The sum of these point sources will equal displacement
on a finite-length fault. Note that if the rupture velocity v is not
defined (or is less than or equal to 0.), then the rupture velocity will
be set to 0.8 times the minimum S velocity in the model (or 0.5 times
the minimum P velocity if the S velocity is 0.).

Ideally, an earthquake fault will have a physically reasonable rupture 
time-history. The default Ricker wavelet is not physically reasonable, 
because it will require the slip to change from positive to negative and
back to positive. Therefore, it is prudent to use a predefined time-history
as given by the "file" attribute (it will be a SAC file). A Gaussian 
time-history is reasonable. Keep in mind that the time-history drives the
stresses (or velocities), so a Gaussian source function will be equivalent 
to a ramp-like function in displacement.

The all important question is how the input amplitude ("amp") is related to 
the seismic moment of the fault. This is a little tricky. The input 
time-history function drives the stresses (or velocities), and hence
represents a displacement rate. The total displacement (which is proportional
as the seismic moment), will be be proportional to the integral of the
time-history function. Hence, the amplitude attribute in the input file
will be set equal to: 

  amp = desired_seismic_moment/points_on_fault/time_history_integral

where the units are in "dyne-cm".

For a vertical fault that is parallel to a grid face, the number of points 
will equal "(length/dh + 1)*(width/dh + 1)", where dh is the spacing of the 
grid. For a vertical fault not aligned with the grid, the number of points
will equal "(length'/dh + 1)*(width/dh + 1)", where length' is equal to
"MAX(abs(x2-x1), abs(y2-y1))" where "x1,x2,y1,y2" are the coordinates at
the ends of the faults. For a non-vertical fault, width' is used in the
above equations where width' is equal to "ztop - zbot" (for faults dipping 
greater than 45 degrees), or width' equals the horizontal distance 
representing the projection of the fault width to the surface (for faults 
dipping less than 45 degrees).

For "mode=7", it's important to realize that the moment (or moment rate) is
constant over the fault plane. If the rigidity is non-uniform across the fault 
plane then the slip will be non-uniform. This can create problems for 
low-rigidity regions near the surface, since the slips can be non-physically
large (because moment is constant). Hence, it's often necessary to use
distributed slip (really distributed moment) models defined by modes 8 and 9.

For "mode=8", a file is read which defines the location, initiation
time, and moment for each point acting as a source along the fault.
The locations are given at nodal values, the initiation time is given
in time-steps, and the amplitude is unitless, since the units come from
the "amp" attribute. This file is defined by the "dfile" attribute, and the 
number of lines in the file is defined by the "p" attribute (it must exactly 
match the number of lines in the file, and there can be no blank lines or 
comments). The first 5 lines of a distributed slip file might look like this:

 192   48  379     0 2.327e+02
 192   47  378     8 2.302e+02
 192   48  378    93 2.251e+02
 192   49  378     8 2.208e+02
 192   49  379     8 2.208e+02

The first column is the l grid number (y direction) of the point source 
(or sub-fault) to be activated, the second column is the m grid number 
(depth direction), the third column is the n grid number (x direction), 
the 4th column is the time-step of source initiation, and the last column 
is the amplitude of the point source. The initiation time is an offset time
(in time-steps) from the time-history file (defined with the "file" 
attribute). An offset of 0 means that the time-history file will be used
as is. An offset of 8 means that the first value in the time-history file 
won't be used until 8 time-steps into the simulation. The amplitude of
each point source (or sub-fault) is given. The simulated earthquake moment
will equal the amplitude of the source ("amp" attribute), multiplied by the 
sum of the moments in the fault distribution file, divided by the integral of 
the source time function.

For "mode=9", everything is similar except that the rise-time is also
included in the fault distribution file. In this case, the time-history
file ("file" attribute) is not used. Instead, the source time-history is 
modeled as a Gaussian pulse. The Gaussian pulse is given by:

  g = amp*exp(-(PI/rise_time*t)**2)

where the time t is determined from 

  t = (timestep - source_initiation)*dt - 0.5*rise_time - t0

where source_initiation is the 4th column and t0 is the "t0" input attribute.

An example fault distribution file with rise time is given by:

 189    0  173   254 1.250e+19   1.000
 189    0  174   257 6.690e+19   1.000
 189    0  175   259 1.217e+20   1.000
 189    0  176   262 1.766e+20   1.000
 189    0  177   264 2.316e+20   1.000

In this case, the rise-time is 1.0 seconds. This is approximately the
width of the Gaussian pulse. For "mode=9", the seismic moment is equal to
"amp*SUM(input_moments)". Because the rise-time can be variable, it's integral
needs to be calculated within the code and hence it has already been factored
in the amplitudes (there is no need to divide by the integral of the 
time-history file). It's typical to have "amp=1.0" on the source line-type
in this case. (Special note: the integral of the Gaussian function defined
above equals "rise_time/sqrt(PI)".)

It's recognized that the input parameters for a finite-length fault are
not straightforward to implement. In the future, the input format will be 
modified so the parameters (slip, fault location, etc) can be input more
easily (for example, eliminate the need to worry about variations in 
rigidity).

Sometimes high frequency noise is entered into the source (and hence the 
solution). This can happen if the source initiates too close to the beginning 
of the simulation. The jump from 0 amplitude at timestep 0 to 
a source amplitude at timestep 1 acts as a step function, and introduces
high frequency noise. This can be avoided by using (or increasing) the
"t0" attribute. When time-history files are used, this can be avoided by
modifying the time-history so that the source initiates at a later time.

Note that the only sources available with the acoustic model (model=2 on
the grid line), are type=1,5. 

Run-Time Visualization
----------------------
A run-time visualization option is available, which is called with the 
"visual" line-type. This allows the wavefield to be viewed (via X Windows)
as the simulation progresses. The P and S potentials are displayed, so
one can inspect mode conversions at the surface and medium boundaries.
The P potential is red (red-blue for negative values) and the S potential
is green (green-blue for negative values). The "movie" attribute sets
the frequency (in timesteps) at which the visual display is updated. 
Currently, the amplitude of the display is not automatically adjusted. 
Hence, the display can be void of any meaningful information. The "scale"
attribute can increase/degrease the display amplitude so it can be 
observed. The scale value must often be determined by trial and error.
The "sample" attribute disseminates the visual display by the sample
factor. This is useful for large models, when the display would otherwise
cover the entire screen. An "x or n", "y or l", or "z or m" attribute
determines the dimension to display, and the plane position. This is
useful only for 3-D problems. For example, "x=10." will plot the YZ 
plane at a distance of 10 km from the x origin; "m=0" will plot the plane
at the surface. The "mode" attribute is not used. The runtime visualization 
is useful for debugging (both the code and the simulation model). It can be 
computationally expensive, so it should not be used for large problems.

Sometimes there is an "arithmetic exception" with the run-time visualization.
This is due to numeric overflow when the wavefield is converted into a raster
image. Reducing the "scale" on the "visual" line-type can eliminate this,
as can turning off the run-time visualization.

Seismogram (SAC) Output
-----------------------
Seismograms at individual points (traces) can be output. The files are output
in SAC format. For each point, several independent files will be written. For
2-D simulations, SAC files representing the Vx velocity, Vz velocity, 
P potential, and S potential will be written. For 3-D simulations, the Vy
velocity also will be written. These files have extensions x, y, z, p, 
and s. The "sac" line-type is used to write seismograms. The attributes are
the x, y, z, l, m, or n coordinates, "scale" (not used), "sample", "mode",
and "file". The velocity files will have units of m/s (the potential file
are unitless).

If the coordinates are given as l, m, or n, the Vx, Vy, and Vz velocities will 
be defined at the Vx variable for that grid node. The P and the S potentials
will be defined at the Tnn (normal stress) variables for that grid node. If 
coordinates are given as x, y, or z (in km), the Vx, Vy, and Vz velocities will
be defined at the exact point in the grid. Interpolation will be used if the
point does not fall on a grid point (which will likely be the case). The P and
S potentials will be defined at the Tnn (normal stress) variables for that
grid node.

The "sample" attribute will sample the trace data at that rate (default is 1,
which is every time-step). The "mode" attributed flag defines how the SAC files
will be output. If "mode=1", then velocities and potentials will be written
(the default). If "mode=2", then moment tensor components will be written
(this is good for utilizing seismic reciprocity). Currently, the "mode=2" 
option is not supported and shouldn't be used. The "file" attributed, 
defines the header for the file name. If this attributed is not specified,
then the header name "sac" will be used.

File names have the format headername.filenumber.filetype. Header name is
defined with "file" (or is sac). The filenumber is a unique number defining
the file name, starting at 0. The filetype is x, y, z, p, or s. Typical
file names might be "sac.000.x", "hywd.031.p", etc.

For example, the following line will write a sac file for a 2-D simulation:

  sac x=10. m=0 sample=2 file="hywd"

The selected point will be at the surface and at a distance 10 km from the
grid origin. Every other time step will be recorded. The files will be named
"hywd.000.x", "hywd.000.z", "hywd.000.p", and "hywd.000.s". Output units are
in m/s.

Care must be used when defining the trace location when topography is used.
The point should be at least one grid point below the topographic surface.

The auxiliary routine "readsac" will read a SAC file and convert it to ASCII.

In 2-D, the potentials at grid point (m, n) are defined as: 

  P = dVx/dx + dVz/dz    
  S = 0.5*(dVx/dz - dVz/dx)

where dx = dz = dh, 
and in 3-D:

  P = dVx/dx + dVy/dy + dVz/dz
  S = 0.5*((dVz/dy - dVy/dz) + (dVx/dz - dVz/dx) + (dVydx - dVxdy))

Note that Vz is positive down.

Note that mode=2 cannot be used with the acoustic model (model=2 on the
grid line).

Traces Output (Multiple Seismograms)
------------------------------------
In addition to seismogram output at individual points, seismograms at 
multiple points (e.g., hundreds) can be specified using the "traces"
input line. This is useful for the type of problems involved in oil
exploration, where one is interested in the output along a densely
sampled array or two dimensional grid. The Vx, Vy (3-D), Vz velocity
components can be output, as well as the P and S potentials. The
characteristics of these seismogram traces are similar to the SAC
output option.

Currently, there are two ways to specify traces: 1) listing the x,y,z
coordinates of the trace locations in an ASCII input file; 2) specifying
the starting and ending x,y,z coordinates of a "trace line", and the number 
of traces that fall on that line. Both methods can be used together.

The user can specify an input file that contains the x, y, and
z coordinates of the traces to be output (in kilometers from the grid
origin). Each trace coordinate location is listed on a separate line,
and any line not beginning with a data value is regarded as a comment
and is ignored. For example, a trace coordinate file might be 
given by ...

   45.0 85.4 25.4
   46.0 86.4 25.4
   47.0 87.4 25.4
   48.0 88.4 25.4
   49.0 89.4 25.4
   50.0 80.4 25.4
   51.0 91.4 25.4
  # comment
   another comment: following line is blank

  10.001 15.0 0.05
  10.001 15.0 0.10
  10.001 15.0 0.15
  10.001 15.0 0.20
 
This will define a horizontal line of 7 seismogram traces starting at
(x, y, z) = (45.0, 85.4, 25.4) (in kilometers from grid origin) and 
ending at (51.0, 91.4, 25.4). There is an additional vertical array
of 4 seismogram traces starting at a depth of 0.05 km and extending to
a depth of 0.20 km. Normally, the trace coordinate file will contain 
hundreds or thousands of trace locations.

In 2-D problems, the y coordinate is still specified but should be set
equal to 0.0. In 2-D or 3-D problems, if a value falls outside the grid,
it will not be directly output, although it will be output as a zero
trace in the demultiplexed file after post-processing (see below).

This trace coordinate file is specified on the "traces" line of the
parameter input file with the "tfile" attribute (tfile="trace_locations").

In the other method the user specifies the x, y, and z coordinates
of the endpoints of a line (km). The user also specifies the number of
traces that will fall on that line (inclusive of the endpoints).
For example, if the user specifies ...

  traces x1=1.0 x2=10. y1=3.0 y2=12.0 z1=0.1 z2=1.0 p=10

In this case, 10 traces will be constructed starting at (x, y, z) = 
(1.0, 3.0, 0.1) and ending at (x, y, z) = (10.0, 12.0, 1.0). In this
example the traces will be spaced 1.0 km in the x and y directions and
0.1 km in the z direction. Note that p is inclusive of the endpoints
and therefore p must be greater than 1.

The traces file method and the coordinates method can be used together.
In this case, the traces specified by the coordinate method will be 
written first in the output file.

Other attributes on the "traces" line of the input parameter file are
"file", "sample", "scale", and "mode". The "file" attribute specifies
the header of the single binary output file that is created (default
is "traces"), "sample" specifies the sampling increment in timesteps
to write each trace (default is 1), "scale" is not used, and "mode"
specifies what components to output in a bitwise logical output format:
mode=1 is Vx component, mode=2 is Vy component, mode=4 is Vz component,
mode=8 is P potential component, and mode=16 is S potential component.
For example, mode=7 will output the Vx, Vy, and Vz components, but will
not output the P and S potential components. The default is to output
all components which is mode=31.

Each component will be output to a separate file, which will end 
with ".Tcomponent" (e.g., .TVx, .TVy, .TVz, .TP, .TS).

Multiple traces files can be specified on the parameter input line. 
For example, one might specify a large 2-D horizontal grid of traces
on one traces line in the parameter input file, and a vertical array
on another traces line of the parameter input file. An example is ...

  traces file="outH" tfile"horizontal_grid" sample=1 mode=8
  traces file="outV" tfile"vertical_array"  sample=4 mode=7

In this case, 4 files will be created. The P potentials at the horizontal
grid of trace coordinates will be output to a file called "outH.0.TP",
and the Vx, Vy, and Vz components in the vertical array will be output
(sampled every 4th timestep) in the files "outV.1.TVx", "outV.1.TVy",
and "outV.1.TVz". Note that each grouping has an increment number, 
which avoids the potential for name clashes.

The output file(s) created with be multiplexed. That is, it will be
in timestep order. The files can be demultiplexed into trace ordering 
using the auxiliary routine "getdtraces". This routine can also extract 
a single trace and output it as a SAC file. Note that trace locations
falling output the grid will not be included in the multiplexed file,
but will be included as a set of zero seismograms in the demultiplexed 
file.
 
Image Output
------------
Various types of "image files" can be written using the "image" line-type.
Attributes to this option are "t", "movie", "scale", "sample", "file", "mode",
and "x, y, z, l, m, n" coordinates.

In general, the "t" attribute defines a specific time at which an image 
(or series of images) will be written. The "movie" attribute defines the 
increment in timesteps at which a series of images will be written. The 
"scale" attribute should not be used. The "sample" attribute defines the 
spatial sampling factor for the image (useful for big images to save disk 
space). The "file" attribute defines the header name for the image (the 
default name is "image"). For 3-D simulations, an "x, y, z, l, m, n" attribute 
defines the dimension and plane position of the image to be written. For 
example, "x=5." would output an image in the YZ plane shifted 5 km from the 
x origin.  The "mode" attribute defines the type of image to be written.

With one exception (mode 3), the images are output as floating point sets
of data. Conversion or auxiliary routines take this data and convert it
into true image formats, such as Sun Raster Files.

Types of "images" include: P and S potential wavefields (mode 1), Vx, Vy, and
Vz wavefields (mode 2), 8-bit raster files (mode 3: old option, don't use),
maximum horizontal velocity map (mode 4), cumulative energy map (mode 5), 
p, s and density images of the velocity model (modes 6-8), and combined 
"framefiles", which contain information at several time steps regarding the 
state variables (Vx, Vy, Vz, Txx, Tyy, Tzz, Txy, Txz, and Tyz) (modes 11-19),
and the potentials (P and S) (modes 20 and 21).

Typical image output lines are:

  image t=100   x=10.  mode=1  file="cross"
  image movie=500 m=0  mode=2  file="surface"

The first will produce a cross sectional image at timestep 100 of the p and s 
wavefields in the YZ plane, 10 km from the x origin. The files will be named
"cross.100.p" and "cross.100.s". The second line will produce a series of 
mapview images, every 500 timesteps, through the surface grid nodes ("m=0"). 
The images will be the x, y, and z wavefields. For 2-D problems, the plane 
coordinates (e.g., "x=10." and "m=0" do not need to be specified).

The P and S potential wavefield images (mode 1) are useful for viewing mode
converted energy. This is the default mode. They contain the P and S 
potentials at selected timesteps.  From a visualization standpoint, this is 
the primary form of output. These floating point image representations can be 
turned into raster images through the auxiliary routines "getimage", or more 
importantly "getps" (which combines the P and S wavefields into a single 
image). These files will be named "imageheader.timestep.p" and 
"imageheader.timestep.s". If the imageheader header isn't specified (using
the "file" attribute), then the default imageheader will be "image". Typical 
file names might be "hywd.0100.p" and "hywd.0200.p", which provide the p 
wavefield at timesteps 100 and 200.  The "t" attribute will write an image 
at a given timestep; the "movie" attribute will write an image every "movie" 
timestep.

In 2-D, the potentials are given by: 

  P = dVx/dx + dVz/dz    
  S = 0.5*(dVx/dz - dVz/dx)

where dx = dz = dh, 
and in 3-D:

  P = dVx/dx + dVy/dy + dVz/dz
  S = 0.5*((dVz/dy - dVy/dz) + (dVx/dz - dVz/dx) + (dVydx - dVxdy))

The X, Y, and Z wavefield images (mode 2) contain the x, y, and z velocities
(in m/s) at selected timesteps. Typical file names might be "hywd.0100.x",
"hywd.0100.y", and "hywd.0100.z". These files can be turned into raster 
images with the auxiliary routine "getimage". The "t" attribute will write 
an image at a given timestep; the "movie" attribute will write an image 
every "movie" timestep.

The maximum horizontal velocity (mode 4) and the cumulative energy array
(mode 5) computed during a simulation can be written as image files. These
modes can only write images in the XY planes (map view). The maximum 
horizontal velocity is determined by sqrt(MAX(Vx*Vx + Vy*Vy)), and is
in m/s. The cumulative energy is determined by the algorithm
0.5*density*SUM(Vx*Vx + Vy*Vy + Vz*Vz). The movie attribute determines how
often (frequency in timesteps) these algorithms are updated. For "mode=4",
a movie value of 5-10 is reasonable. For "mode=5", a movie value of 1 is
reasonable. The "t" attribute is not used. A single file is written at the 
end of the simulation, and is named "imagefile.totaltimesteps.V" for mode 4 
(e.g., "hywd.6000.V") and "imagefile.totaltimesteps.E" for mode 5 (e.g., 
"hywd.6000.V"). The auxiliary routine "getsubfield" is used to convert
this file into a standard image format, which can then be used with the
auxiliary routine "getimage" to convert it into a raster file (Note: for
2-D problems, a SAC file will be created representing the maximum horizontal
velocity or energy at some depth. This is an old option, and probably 
shouldn't be used.) The sample factor may need to be 1 (or unset) for
image mode 4 and 5 to work (possible software issue in the code).

Image files representing the medium parameters (P and S velocity, and the
density) can be written (modes 6, 7, and 8). The medium seismic velocities
will be in km/s, and the density will be in gm/cm**3. In a typical case, the 
code will be run for only 1 timestep, and the "t" attribute will be set to 1.
For example, the following line-types might be used:
 
  time  t=1 dt=0.015
  image t=1 mode=6 m=0
  image t=1 mode=8 x=10.

This will create a file "image.1.A", which represents a map view of the P 
wave velocity at the surface, and a file "image.1.R", which represents a 
cross section of the medium density in the YZ plane 10 km from the x origin.
These files will be named "imagefile.timestep.extension", where the 
extension is A (P wave velocity - alpha), B (S wave velocity - beta), or
R (density - rho). These files can be turned into raster images with the
auxiliary routine "getimage" (on massively parallel machines, these files
must first be converted with the auxiliary routine "getsubimage").

Note that on output, the velocities are converted from km/s into m/s, and
the density is converted from gm/cm**3 to kg/m**3. In addition, note that
the medium values may not be numerically equal to the input values. That is, 
there may be small roundoff issues. This is because the medium parameters are 
first converted into elastic parameters (mu and lamda) and the inverse of 
density, and on output are converted back into the medium parameters. This 
is normally not a problem, but can be a puzzling issue in certain 
circumstances.

Cumulative image data representing the state variables (Vx, Vy, Vz, Txx, 
Tyy, Tzz, Txy, Txz, Tyz) and potentials (P and S) can be written using 
modes 11-21. Each file contains time-slices of the selected variable 
written every "movie" time-step. These files are called "framefiles". These 
files will be named "imagefile.extension", where extension is one of the 
state variables. In essence, framefiles represent a catenation of several 
images files (or frames).  They are, however, useful for manipulating 
information. For example, suppose it's desire to compute a surface map of 
the static X displacement field at the end of a simulation. This option 
can be used to write the Vx velocity field along the surface every 10 time 
steps. An auxiliary routine can be used to read the framefile, and integrate 
the velocities to determine the displacements. Another example is if it is 
desired to have a map of the response spectra. This can be done at a single 
point using SAC files and seismogram output, but the framefile option allows 
an auxiliary routine to compute the needed map. The auxiliary routine 
"getframe" can extract a single image at a given time-step from the 
framefile, and the auxiliary routine "getframes" can read the framefiles 
to find an average value, the maximum value, or the integration. 

Note shear stresses (mode=17,18,19) cannot be output with the acoustic model 
(model=2 on the grid line), nor can the medium parameters (mode=6,7,8). Also,
note that since all three normal stresses with the acoustic model are 
identical, normal stress output (mode=14,15,16) will produce the same result.
The files still will have the Txx, Tyy, Tzz extension, however.

Some image files that are written on massively parallel machines must be
converted into a serial format using the auxiliary routine "getsubimage". 
This is for image modes 1, 2, 6, 7, and 8). They can then be converted
into regular images with the auxiliary routine "getimage". This procedure 
will not be needed in the future.

Volume Output
-------------
Volume output is like image output except that 3-D volumes of data are
output. Currently, only P and S potential files can be generated. These
have file modes 1 and 2, respectively (22 and 23 also work). The files
can be very large so the sample attribute should be considered. Volume
output can be generated only with 3-D simulations. One can output 
individual data volumes at selected timesteps (t attribute) or data
volumes at incremental timesteps (movie attribute).

Typical file naming convention would be something like: 

  volume.00900.3DP
  volume.00900.3DS

The first part is the user definable header (default is volume), the
second part is the timestep of output, and the third part defines it
as a 3-D data volume (3D) and the output type (P potential or S potential).

The auxiliary routine "getvslice" can be used to generate an image slice
of the data volume at a given plane and orientation dimension. These data 
slices can be viewed with "getimage" or "getps". 3-D volume files generated
on MPP machines need to be converted to the serial version with the
auxiliary routine "getsubvolume".

Note: on mpp machines the data volumes for each node can be written
to the same file on a global file system or to independent files on
a local file system. If they are written to independent files on a 
local file system then files for all nodes must be combined for later 
processing. They can be combined in any order.

Surface and Absorbing Boundary Conditions
-----------------------------------------
Clayton and Engquist boundary conditions are the primary means used to absorb
wave energy that hits the grid boundaries (this is the default). In addition,
unless specified otherwise, a free-surface condition is applied at the top 
grid boundary. The boundary condition mode is controlled with the "b"
attribute of the "grid" line-type. The default is "b=3" (Clayton and Engquist
on the sides, free surface on the top). If "b=1" then all side boundaries
will have reflecting conditions (not Clayton and Engquist boundary conditions); 
if "b=2" then all side boundaries will have Clayton and Engquist boundary
conditions (but no surface boundary condition). Clayton and Engquist boundary
conditions absorb much, but not all, of the energy. The surface and Clayton
and Engquist boundary conditions are applied at grid points external to the
specified grid. Depending on the boundary, 2 to 3 addition grid points are 
used (e.g., the actual grid is slightly larger than the virtual grid). In
addition, a b=4 boundary condition is available. It operates like b=3 
except that the boundary conditions have been modified to prevent some
numerical instabilities that occur when material with zero shear velocity
intersects the boundary. Reflections using this option will be slightly
larger than reflections with b=3. Also, with the acoustic modeling option
(model=2 on the grid line), b=3 will act as b=4. Note that the b=4 option
is a temporary fix used to prevent instabilities, and is also known as 
the "water" boundary conditions. In general, the boundary conditions should 
be improved.

"Sponge" boundary conditions (Kosloff) can be implemented using the "damp" 
and "adamp" attributes of the "grid" line-type. The damp attribute specifies 
the depth in grid points for the sponge layer around the boundaries. Unlike 
Clayton and Engquist boundary conditions, these points are internal to the 
grid (e.g., they take up the calculational domain). The adamp attribute 
specifies the damping factor. This factor should fall between 0 and 1
(normally, it should be close to 1.). If "adamp=1.", then no damping will 
be applied.  If "adamp=0.", then damping will be 100%. If "adamp < 0", then
a default value will be used which is equal to "1 - 1/(2*damp)". 
Normally, damp might be on the order of 20 and adamp on the order of 0.95. 
However, experimentation is necessary to determine what values of damp and 
adamp to use for the sponge parameters.

The sponge algorithm is to reduce value of the velocity variables at each
time-step by a certain factor. This factor is equal to 

  factor = 1. - (damp - distance)/damp*(1 - adamp)

where distance is the distance in grid points from the edge of the grid.
The velocity variables are multiplied by this factor.

The default is no sponge boundary conditions. In general, sponge boundary 
conditions are not used because they are non-physical, and hence tend to 
reflect energy back into the grid.  In addition, they are somewhat 
computationally intensive. Sponge boundary conditions can be used 
simultaneously with Clayton and Engquist boundary conditions. 

Another method to apply boundary conditions is to use seismic attenuation.
This is similar to using sponge boundary conditions except that the
attenuation algorithm is more physical and works considerably better.
To use this method, one runs the condition with the attenuation option
on, and specifies a region of low Q around the grid boundaries. Typically, 
a Q value on the order of 5 works well. The depth of the attenuating zone
might be 20 - 50 grid points. 

This type of boundary condition can work very well, and will be the
choice of the future. 

Attenuation
-----------
Attenuation is initiated by the "q" attribute of the "grid" line-type. 
If "q=1", then values for Q will be expected (required) as input in the 
velocity model. Q can be entered as either Q, Qp, or Qs (if Q is used then 
Qp = Qs). In addition, Qf should be entered. Qf is the central frequency of 
the source. The attenuation model approximates Q at a range of frequencies. 
Hence, it's necessary to specify the frequency range (or central frequency) 
of interest.  Q can be used as absorbing boundary conditions along the side 
boundaries of the grid (for 2-D and 3-D models). Note that attenuation 
currently can not be used with the pure acoustic model (model=2 on the grid 
line), although it can be used if the shear velocities are set to zero in a 
regular elastic model.

Topography
----------
Currently, it's possible to incorporate topography in 2-D simulations.
A 3-D algorithm will be implemented. Topography is initiated with the
"topo" line-type. The available attributes are "file", where the topographic 
information is read from an input file, and "depth", where the topographic
boundary is defined as a flat layer of a certain depth (in km). The depth
attribute is used for testing. If a file is specified, the depth attribute
will be ignored. Topographic modeling uses a density extinguishing approach.

The "file" attribute specifies a SAC file that contains the topographic
information. The number of data values in this file MUST equal the number
of grid points in the horizontal x direction of the finite-difference 
grid. The values in this file is the depth of the topography (in km).
That is, valley's will have a greater depth than mountains. The depths
can be negative, but in this case the topography will intersect the top
boundary of the grid.

To use topography, the E3D code must be compiled using the routines
paraxial2d.cslow, stencil2d.cslow, and update.cslow. This limitation
will be changed.

Restart
-------
Simulations that end abruptly (by choice or computer failure) can be 
restarted using the restart facility. The restart option has several
attributes. The "dump" attribute defines how often (in timesteps)
the simulation parameters are dumped. The "restart" attributed is
a flag indicating dump files from a previous simulation (or partial
simulation) have been created and that the simulation can be restart
at that point. This attribute is set equal to the timestep number
that created the last simulation dump. The "file" attribute is the
header name given to the dump and/or restart files. The "file" 
attribute is required. The file naming convention is: 

  filename.mppnode.timestep.E3D_RESTART

An example is:

  restart.240.4500.E3D_RESTART

After a dump file for one timestep has been written the dump file for
the previous timestep (if greater than 0) will be deleted. If it's 
desired that no dump/restart files exist at the end of the simulation
then the "nolast" flag attribute can be set to 1. 

Note: there could be problems if the dump interval is not an integral
fraction of the number of timesteps.

Note: currently, the restart option may only work for 3D simulations.

Note: currently, the restart option will not work for attenuation.

Note: currently, the restart option is not available for the acoustic
model (model=2 on the grid line).

On MPP machines it is suggested that the dumps be made to a directory
local to a given node (e.g., /var/tmp).

Because of the large volume of data that can be generated, the restart
facility should be used with caution. In addition, no dump information
regarding seismograms or certain image outputs (mode 4 and 5) are
maintained. In addition, if a crash occurs during the output of certain
data (i.e., traces data, image field data - image modes 11 to 21 - then
these files can be corrupted and must be returned to a noncorrupted
state. Currently, there is no automatic facility to do this.

The restart facility is currently under development and should be used
with caution. Boundary conditions may not be correctly restarted to
exact precision.

Parallelization
---------------
The code has been implemented on massively parallel processors and
work station networks using domain decomposition and a message passing
paradigm using mpsc, mpi, or pvm. Note, on shared memory parallel 
machines, user parallelization is not used. 

The "parallel" line-type is used to describe the parallel grid decomposition.  
The attributes "nx", "ny", and "nz" define how the parallel nodes will be 
distributed along the x, y, and z dimensions. The product "nx*ny*nz" must equal
the number of nodes being used on the parallel machine. For 2-D calculations,
ny must equal 1 (which is the default). On single processor machines (and
shared memory parallel machines), nx, ny, and nz are implicitly equal to 1
(i.e., they don't need to be specified).

Trial and error experimentation is needed to decide how to best decompose a
given problem on a given machine. However, because of the code architecture 
(message passing algorithm and array ordering), it's usually best to 
constrain nx to a small value (such as 1). In addition, other factors play 
a role. For example, the effectiveness of the propagating grid algorithm
is limited for problems run on massively parallel processors. For this
reason, it's often useful to leave 1 dimension undecomposed (e.g., nx=1).

The nodal decomposition directions are x, followed by z, followed by y.
Hence, parallel node 0 will be at the grid origin, and parallel node
(nx*ny*nz - 1) will be the node farthest from the origin.

When nx, ny, and nz are used, the grid is decomposed such that an equal 
number of grid points falls within the domain of a single parallel node.
This type of load balancing works well. For more precise load balancing,
the "file" attribute can be used. The file attribute specifies an input 
file where each line contains the number of grid points in a particular
dimension for a given node. For example, the input lines

  grid l=100 m=100 n=100 
  parallel nx=2 ny=3 nz=2 file="mpfile" 

would have an grid decomposition input file ("mpfile") something like

  50   # x dimension
  50   # x dimension
  35   # y dimension
  30   # y dimension
  35   # y dimension
  55   # z dimension
  45   # z dimension

This will divide the grid into two 50 grid point segments in the x direction,
three segments in the y direction (35, 30, 35), and two segments in the
z direction (55, 45). 

Again, automatic load balancing works well, so the "file" option normally
should not be used.

Some image files that are written on massively parallel machines must be
converted into a serial format using the auxiliary routine "getsubimage". 
This is for image modes 1, 2, 6, 7, and 8). They can then be converted
into regular images with the auxiliary routine "getimage". This procedure 
will not be needed in the future.

Hybridization
-------------
The hybridization option ("join", "xjoin", "tjoin" attributes of the "source" 
line-type) is not recommended at this time. Also, it cannot be used with
the acoustic model (model=2 on the grid line).

Static Grid Refinement
----------------------
The static grid refinement (variable density grid spacing) option is not 
recommended at this time. However, it's specified by using appropriate 
values for the "m" attribute for the "grid" line-type. That is, a value
such as m=10,20,50 has a 10 node top layer of grid spacing dh, a 20 node
middle layer of grid spacing 2*dh, and a 50 node bottom layer of grid
spacing 4*dh. This option is not fully parallelized, and input/output
and other parameter settings are poorly defined.



                          COMPUTATIONAL ISSUES
                          --------------------

Two-dimensional problems require 8 floating point values per grid point,
and three-dimensional problems require 12 floating point values per grid 
point (actually, it's currently 9 and 13 but this limitation will be 
changed). Two-dimensional problems require 64 floating point operations
(no divides) per grid point at every time step; three-dimensional problems
require 141 floating point operations. A speed of approximately 50 MFlops 
is observed on a Sparc Ultra 140 processor. Needless to say, a large 
problem tends to run at a higher MFLOPS rate than a smaller problem. That
is, problem size is an important issue for benchmarking.

Some features are number limited (as defined by the "elas.h" file). There
can be a maximum of 100000 timesteps, 100 seismograms, 100 sources, 100
separate image planes, 100 velocity model input types, 100 points in a 2D
polygon for the velocity model, 256 parallel decomposition nodes in any
one direction, a 3 variable density grids.

Binary input/output formats for the E3D code assume 4-byte integer words,
and 4-byte floating point words. I/O issues may need to be addressed when 
the code is implemented on alternate machine architectures (e.g., longer 
word sizes found on many Cray's).

Sometimes there is an "arithmetic exception" with the run-time visualization.
This is due to numeric overflow when the wavefield is converted into a raster
image. Reducing the "scale" on the "visual" line-type can eliminate this,
as can turning off the run-time visualization.

Note that there are two standard formats for binary files. One is "big endian"
and the other is "little endian". Both use IEEE arithmetic but have their
bytes in different orders. "big endian" is used on most Unix machines (e.g.,
Sun and SGI's, most mpp's, etc.) and "little endian" is used with DEC Alpha
chips and on PC's (e.g., Linux). E3D binary files used and generated on one 
type of machine ("big endian") cannot be used on another type of machine
("little endian") unless their bytes are swapped. Byte swapping can be done
with the auxiliary routine "swapbytes". This involves E3D files such as
sac files, image files, and gridded input files used with the vfile option.



                        INPUT AND OUTPUT FORMATS
                        ------------------------

Input file data formats are described in the various sections above. The only
exception are SAC files, which can be used to define source time-histories for
the "source" line-type, and which are used to define surface topography for
the "topo" line-type (2-D). (SAC files are also used to output maximum 
velocities and energy data for 2-D simulations.) SAC files are used to output 
seismograms, and their format is included below.

SAC Files
---------
SAC files are binary files, comprised of a 158 word header (4-bytes per word),
followed by a specified number of 4-byte floating point values representing
time-series data. The number of time series data is specified in the SAC 
header. The header includes other information, although there are only 
two values that are of significance. They are NPTS (the number of time-series
data that follow the header), and DELTA (the time sampling increment of the
data in seconds). DELTA is the first word in the header, and is represented 
by as a floating point value; NPTS is 80th word (offset 79) in the header,
and is represented as an integer value.

The SAC file units for velocity output (for E3D) is m/s.

SAC stands for Seismic Analysis Code, and is the most common seismic utility
tool used by research seismologists in the world.

Traces Files
------------
Files output using the "traces" seismogram output option are binary files.
They consist of a 3-word (4-byte/word) header, ordering information 
terminated with a negative value (ordering information is necessary for
parallelization issues), and the seismogram (trace) values in a multiplexed
or timestep ordered sequence. The format is given by ...

  global_number_of_traces  (int)
  timesteps (T)            (int)
  time_increment (seconds) (float)
  order0                   (int)
  order1                   (int)
    .
    .
    .
  orderN-1                 (int)
  -1                       (int)
  trace0___timestep0       (float)
  trace1___timestep0       (float)
         .
         .
         .
  traceN-1_timestep0       (float)
  trace0___timestep1       (float)
  trace1___timestep1       (float)
         .
         .
         .
  traceN-1_timestep1       (float)
         .
         .
         .
  trace0___timestepT-1     (float)
  trace1___timestepT-1     (float)
         .
         .
         .
  traceN-1_timestepT-1     (float)

Note that the number of global traces may be greater than the number of
traces given in this multiplexed output file. This is because some trace
locations could be outside the domain of the finite-difference grid. The
N integer values of ordering information contain the trace ordering as
listed in the original trace coordinate file specified on the traces line
of the parameter input file. This information is necessary because of the 
decomposition used on parallel machines. The output order may not be the 
same as the input order. The number of ordering values will indicate how
many trace locations are inside the grid, which may or may not be equal
to the global number of traces. This ordering sequence ends with a -1
(integer). This header information is followed by the seismic amplitudes
in a multiplexed or timestep order.

Note that if both the coordinate output method and the traces file method
are used to define trace coordinates for a single traces file, the traces
specified by the coordinates will be listed first.

The auxiliary routine "getdtraces" can be used to demultiplex this output
file (trace ordering) and reorder the traces so that they are in the orginal 
order as specified by the trace coordinate input file. In this reordered
file, traces that are outside the grid will be inserted with null values. 
Hence, this demultiplexed file can be larger than the original output file.
This demultiplexed file will not have any header information, and hence will
have a size equal to the number of traces times the number of timesteps 
(assuming the sampling rate is 1). The "getdtraces" auxiliary routine also
can be used to output a single trace as a SAC file. This is useful for
testing.

Image Files (non MPP machines)
------------------------------
Image files are created using the "image" line-type (and modes 1, 2, 6 - 8).
Also, they are created by auxiliary routines that are used on file formats
created by other image output modes. For example, "getframe" and "getframes"
operate on files created using image modes 11-21. The output from these 
routines are image files; "getsubfield" operates on files created using image
modes 4 and 5. The output from this routine is an image file. 

Image files are a floating point representation of an image (2-D visualization
plane representing a physical quantity). Image files are not raster images
themselves (such as Sun Raster files), although they can be converted into
raster images using routines such as "getimage" and "getps".

Image files have a simple binary format. There is a 2 integer header (4 bytes
per integer), followed by an xy grid of floating point values (4 bytes per
value). The size of the xy grid (number of floating point values) is given
by the header. The first integer value in the header gives the number of
points in the x dimension, the second integer value gives the number of 
points in the y dimension. In this case, xy refer to the convention of 
raster images (x is positive to the right across a computer image; y is 
positive down). Each "row" of x values (constant y value) is sequential
in the file. That is, the first row (top of the image) will be followed by 
the second row, which will be followed by the third row, and so on.

For example, if the image has a x dimension of 100 and a y dimension of 50,
then the first 4-bytes in the header will be an integer with a value of
100, and the second 4-bytes will be an integer with a value of 50. These
header values will be followed by 5000 4-byte floating point values 
representing the magnitude of the variable of interest. The first 100 floating
point values will be the top row of the "image". The final 100 floating point
values will be the bottom row of the "image". 

Again, it is necessary to run an auxiliary routine such as "getimage" or
"getps" to convert these floating point image files to true raster images
(Sun Raster images). These raster images can then be further manipulated
with public domain software such as "xv" or "imconv" (from the San Diego
Supercomputer Center). For example, although Sun Raster files are a standard
raster file format (they are recognized on many systems, not just Sun's),
"xv" or "imconv" can convert them into other formats such as "gif" or 
"SGI rgb".

When images are output in the xz or yz planes, they are correctly displayed.
When they are output in the xy plane (map view), their orientation is 
reversed. This is due to the difference between the finite-difference grid 
coordinate system and the coordinate system used to define images (y is down).
Hence, these images need to be reversed along the x axis. This is a minor
inconvenience. This reversal can be done with utilities like "xv" or "imconv". 
Nonetheless, this will not be needed in the future. 

Image Files (MPP machines)
--------------------------
On massively parallel processors, and on problems run with static grid 
refinement (variable density grids), the initial image file output is in 
the MPP image file format. This is done to facilitate efficient image output 
during the simulation. This MPP image file format is converted into regular 
image file format (see above) using the auxiliary routine "getsubimage".

The MPP image file format is as follows:

      x   : x dimension of unsampled global grid    (4   bytes int)
      y   : y dimension of unsampled global grid    (4   bytes int)

   followed by n "subgrids" of format

      x   : x dimension of unsampled subgrid        (4   bytes int)
      y   : y dimension of unsampled subgrid        (4   bytes int)
      x0  : x global starting coordinate of subgrid (4   bytes int)
      y0  : y global starting coordinate of subgrid (4   bytes int)
      f   : sampling factor                         (4   bytes int)
      data: field data                              (4*n bytes float)
            where n = ((x-1)/f + 1)*((y-1)/f + 1)

Each MPP node (or each static grid) writes it's own "subgrid". The
term "unsampled" refers to the grid size if a sampling factor of "1"
("sample=1") is used for the "image" line-type.

Again, these files are converted into standard image file format using
the auxiliary routine "getsubimage". These converted image files then can be
turned into raster images using auxiliary routines like "getimage" and 
"getps". Hence, the user need not be concerned with this more complex 
representation.

Field Files
-----------
Field files are written by modes 4 and 5 of the "image" line type. They
have the following format, irregardless of whether they are written on
MPP or non-MPP machines:

      x   : x dimension of unsampled global grid    (4   bytes int)
      y   : y dimension of unsampled global grid    (4   bytes int)

   followed by n (n = x*y) tuples of the format:

      x   : x global starting coordinate            (4   bytes int)
      y   : y global starting coordinate            (4   bytes int)
      v   : field value of interest                 (4   byte float)

These field files are converted into regular image files using the 
auxiliary routine "getsubfield". These converted image files then can be
turned into raster images using auxiliary routines like "getimage".
Hence, the user need not be concerned with this more complex representation.

Frame Files
-----------
It's currently beyond the scope of this documentation to describe the format
of "framefiles" (output created by modes 11-21 of the "image" line-type).
However, they are somewhat similar to the MPP image file format.

Volume Files
------------
Volume files on serial machines have the following format: 9 integer
header followed by floating point values representing the data volume.
The header values are gx, gy, gz, sx, sy, sz, x0, y0, z0. The values
gx, gy, and gz are the dimensions of the data volume that was created.
On serial machines sx=gx, sy=gy, and sz=gz, and x0=0, y0=0, and z0=0.
The floating point values are arranged such that the x dimensions varies
the fastest, then the z dimension, then the y dimension.

On MPP machines then the data are composed of subvolumes, where each
subvolume represents a node on the machine (i.e., process). These blocks
of data are catenated together one after another. Each subvolume has a 
9 integer header just like with serial output. Each subvolume has a
set of floating point values following the header just like with serial
output. In this case gx, gy, and gz are the dimensions of the global data 
volume (which will be the same for all subvolume headers), sx, sy, and sz 
equal the dimensions of that particular subvolume, and x0, y0, and z0 equal 
the offsets of that particular subvolume within the global 3-D data set.
That is, there will be sx*sy*sz floating point values following the 
header for each subvolume.

The auxiliary routines "getvslice" and "getsubvolume" can be used to
process the volume data.



                            AUXILIARY ROUTINES
                            ------------------

A number of auxiliary routines are used with the E3D code. These auxiliary
routines are used to display and manipulate image output, convert SAC files
to ASCII files (and the reverse), and to compress/uncompress gridded 
velocity files. There are several other routines, but they are not 
described here.

Typing the auxiliary routine with no arguments will usually print a brief
usage message about the routine.

getimage
--------
"getimage" converts floating point image files into Sun Raster 24/32-bit
images. These floating point image files are created with modes 1, 2, and 
4 - 6 using the "image" line-type option in the E3D simulation, and they
are output from auxiliary routines "getsubimage", "getsubfield", "getframe",
and "getframes" (these auxiliary routines deal with other file formats
and convert them into the standard image file format). 

Sun Raster files represent a standard image file found on many computer
systems (not just Sun's). Nonetheless, they can be converted into other
file formats (such as "gif" or "sgi rbg") using public domain software
such as "xv" or "imconv" (from the San Diego Supercomputer Center).

The usage of "getimage" is:

  getimage [optional parameters] < imageFile > SunRasterFile

  parameters:
  scale  -     scale factor (default 0: autoscale)
  ascale - autoscale factor (default 100.)
  sample - sampling factor (default 1)
  tfile  - topography file (default none)
  mode   - 32 bit mode (if mode = 32)
  file   - color table file (GMT format)

  example: getimage sample=2 scale=2.0 < image.1000.x > image.1000.x.ras

The "scale" option scales the image. If it is selected then the image will
be scaled by this amount. If it is not selected, then the scale factor
will be determined based on the maximum value in the image file, multiplied
by the "ascale" factor (whose default is 100). The scale factor used will
be output to the screen.

The "sample" parameter determines the spatial sampling factor of the image.
The "tfile" parameter references the topography file used when E3D is run
with the topography option (hence the image above the topography will be
whited out). The "mode" writes the image in either 24-bit or 32-bit format
(default 24-bit).

The example takes the image output from the E3D code that represents the 
velocities in x-direction, reduces the image size by a factor of 2 (in
each direction), scales the image by 2.0, and writes out a Sun Raster 
file. The Sun Raster file can be viewed with public domain utilities such
as "xv", or it can be converted into other raster image file formats (such
as "gif" or "sgi rgb") using "xv" or "imtools" (from the San Diego 
Supercomputer center).

On a first pass, it's usually best not to scale the image. Instead, use
the autoscale option (e.g., don't use the parameters "scale" or "ascale").

The default is to have the created images be shades of red. 

If the "file" attribute is specified, then a colortable file will be
input. In this case, the autoscale option will not be active and the default
scale factor will be 1. It's probably not a good idea to use the scale factor.
The format of the colortable file is GMT format. In this format, an ASCII 
file is read that defines the mapping between image value and color. An 
example colortable file is:

 # Example color table file
 #   V1    R1  G1  B1    V2    R2  G2  B2
   000.00 255 255 255  200.00 255 128 128
   200.00 255 128 128  400.00 255   0   0
   400.00 255   0   0  600.00   0 255   0
   600.00   0 255   0  800.00   0   0 255
   800.00   0   0 255 1000.00   0   0 255
   F          0   0 255
   B        255 255 255
   N        255 255 255

In this case, image values falling between the listed floating point values
(V1 and V2) will have a RGB (red/green/blue) color assignment linearly 
interpolated between the given color indexes. For example, an image value
of 200.0 will have a (red,green,blue) color mapping of (255,128,128). An
image value of 300. will have a color mapping of (255,64,64). An image
value less than the minimum value set in the color table file will have
a color mapping defined by F (default is (0,0,0)). An image value greater 
than the maximum value will have a color mapping defined by B (default is
(255,255,255)). An image value between the minimum and maximum values but
not defined by the colortable will have a color mapping defined by N
(default is (0,0,0)). The values specified in the colortable need not
be continuous, and there can be gaps. If there is an overlapping conflict, 
the latter value will be used. For each line, the 2nd floating point value
must be greater or equal to the first. The color indexes must be between
0 and 255. The colortable file can contain comments, which are any lines
not beginning with a number (0-9), ".", "-", "+", or "F", "B", or "N".
The maximum number of entries specified in the colortable file is currently
set at 1000.

getps
-----
"getps" takes image files representing the p and s potentials and converts
them into a single Sun Raster file that shows both potentials simultaneously.
These image files are created using "mode=1" from the "image" line-type
during the E3D simulation.

This is the primary way to view images, because the interaction and mode
conversions between P and S energy can be observed. The P potential will
be red (red-blue for negative values) and the S potential will be green
(green-blue for negative values).

Sun Raster files represent a standard image file found on many computer
systems (not just Sun's). Nonetheless, they can be converted into other
file formats (such as "gif" or "sgi rbg") using public domain software
such as "xv" or "imconv" (from the San Diego Supercomputer Center).

The usage is:

  getps [optional parameters] PImagefile SImagefile > SunRasterFile

  parameters:
  scale  -    scale factor (default 0: autoscale)
  pratio -  p scale factor (default 1: autoscale)
  sample - sampling factor (default 1)
  color  - color scheme (1 = black bg; 2 = white bg) (default 1)
  clip   - clipping color value (default 255 color=1; 0 color=2)
  mode   - 32 bit mode (if mode = 32)
  tfile  - topography file (default none)

  example: getps sample=2 image.1000.p image.1000.s > imageps.1000.ras

The "scale" option scales the image. If it is selected then the image will
be scaled by this amount. If it is not selected, then the scale factor
will be determined based on the maximum value in the image file, multiplied
by an "ascale" factor (which is 100). The scale factor used will be output 
to the screen.

The "pratio" parameter is a factor to multiply the P potential relative to 
the S potential (the default is 1). This would be done to highlight P energy,
for example (or highlight S energy if "pratio" is less than 1.). 

The "color" parameter determines the background color scheme (white background
or black background). The "clip" parameter represents the maximum raster value
(0 to 255) for the red, green, and blue colors in the raster image.  The 
default clip value is 255 for a black background ("color=1") and 255 - clip 
for a white background ("color=2").

The "sample" parameter determines the spatial sampling factor of the image.
The "tfile" parameter references the topography file used when E3D is run
with the topography option (hence the image above the topography will be
whited out). The "mode" writes the image in either 24-bit or 32-bit format
(default 24-bit).

The example takes the image file output from the E3D code that represents the 
P and S potentials, reduces the image size by a factor of 2 (in each 
direction), and writes a Sun Raster file which shows the P and S potentials
simultaneously. The Sun Raster file can be viewed with public domain utilities 
such as "xv", or it can be converted into other raster image file formats 
(such as "gif" or "sgi rgb") using "xv" or "imtools" (from the San Diego 
Supercomputer center).

On a first pass, it's usually best not to scale the image. Instead, use
the autoscale option (e.g., don't use the parameters "scale"). 

getsubimage
-----------
"getsubimage" converts image files created on an MPP machine into the standard
image file format. These files would be created by modes 1, 2, and 6 - 8 using
the "image" option during the E3D simulation. These converted files can be
turned into raster images using auxiliary routines such as "getimage" and
"getps". The usage is:

  getsubimage [optional parameters] < e3dSubImageFile > imageFile

  parameter:
  sample - sampling factor (default autosample)

  example:
  getsubimage sample=1 < subimage.0500.p > image.0500.p

This example converts the image (representing the P potential) into an
image file called "image.0500.p". A sampling factor of 1 means that
the image file will be sized as if "sample=1" was used on the "image" line
during the E3D simulation. The file "image.0500.p" then can be converted
into a raster image with an auxiliary routine such as "getimage", and can
be used with auxiliary routine "getps".

getsubfield
-----------
"getsubfield" converts files created by modes 4 and 5 of the "image" line-type
into standard image files. These standard image files then can be converted
into raster images using auxiliary routines such as "getimage". The usage is:

  getsubfield [optional parameters] < e3dSubFieldFile > imageFile

  parameter:
  sample - sampling factor (default autosample)

  example:
  getsubfield sample=1 < subimage.1000.V > image.1000.V

This example converts the image (representing maximum horizontal velocity)
into an image file called "image.1000.V". A sampling factor of 1 means that
the image file will be sized as if "sample=1" was used on the "image" line
during the E3D simulation. The file "image.1000.V" then can be converted
into a raster image with an auxiliary routine such as "getimage".

getframe
--------
"getframe" is used to extract individual timestep images from a "framefile"
(image output using modes 11-21). The output will be an image file. Usage is:
  
  getframe [optional parameters] < e3dFrameFile > imageFile
   
  parameters:
  frame  - frame number to extract (start 1; 0 will output header)
  sample - sampling factor (default 0 == autosample)

  example:
  getframe frame=5 sample=1 < image.Vx > image.1000.Vx

In this example, getframe will output the 5th image frame in the file
(which happens to be timestep 1000. The sampling factor is 1, so that 
the original image size is maintained.

The output image ("image.1000.Vx) can be viewed or manipulated with "getimage".

getframes
---------
"getframes" manipulates an entire "framefile" (image output using modes 11-21).
Average values (as an image) for the simulation can be determined, the values
can be integrated (e.g., to produce an image of the displacements from the
velocities), and the maximum value can be determined. The output will be an 
image file. Usage is:

  getframes [optional parameters] < e3dFrameFile > imageFile

  parameters:
  mode   - 1 average   (average at end)
         - 2 integrate (average at end)
         - 3 maximum   (absolute) [default]
  sample - sampling factor (default 0 == autosample)
  time   - number of frames to average at end (default nframes/20)

  example:
  getframes mode=2 sample=1 time=10 < image.Vx > image.Dx

This example takes a framefile representing images of the X component of 
velocity and integrates them, to produce an image file of the X component
of displacement. "time=10" means that the last 10 frames in frame file,
or more precisely their integral, will be averaged to determine the final
image of displacement. This is done because the velocities near the end of
the simulation still tend to oscillate, so an average of the last several
frames more accurately gives the displacement field.

The output image ("image.Dx) can be viewed or manipulated with "getimage".

getvslice
---------
"getvslice" generates an image slice from a volume dump of data. The
image slice will have the same format as an image created using the
image option. The "getvslice" must be used on a volume file generated
on a single process machine or one that has been converted to this
format using "getsubvolume". The image slice that is generated from
"getvslice" can be converted to a raster image using "getimage" or
"getps". The usage for "getvslice" is:

  getvslice filename dimension plane > imageSlice

In this case, filename is the file create using the volume output 
option (converted using "getsubvolume" if necessary), dimension is
the orientation of the image slice (1 = xy plane, 2 = xz plane,
3 = yz plane), and plane is the plane number in the volume that 
will be extracted 

An example usage is:

  getvslice volume.04000.3DP 2 100 > imagexz.04000.p

This will extract an image slice from the file "volume.04000.3DP". It
will be in the xz plane, and 100 grid nodes from the surface. The
file created, "image.04000.p" can be converted to a raster image with
the auxiliary program "getimage" (or "getps").

Typing "getvslice" with no arguments will print a brief usage message.

getsubvolume
------------
"getsubvolume" converts an e3d volume dump (volume line) generated on an
mpp machine (multiple processes) to one that would have been generated on
a serial machine. It is similar to "getsubimage" except that currently
"getsubvolume" accepts no arguments or parameters. An example usage is:

  getsubvolume < volumempp.04000.3DP > volume.04000.3DP

This converts (reorganizes) the volumempp.04000.3DP file created on a 
multiple process machine (e.g., mpp) into one that would have been created
on a single process machine (e.g., workstation). 

Note: getsubvolume handles everything in core so the volume dump size should
not be greater than the amount of memory on the machine. 

Note: there is no usage message with "getsubvolume".

readsac
-------
"readsac" reads a SAC file and outputs an ASCII file, which has a 1 line 
header and is followed by each time-series data in the SAC file (1 data value
per line). The usage for "readsac" is:

  readsac < SACFile > ASCIITraceFile

writesac
--------
"writesac" takes 2 arguments and reads an ASCII trace file (no header), 
and writes a SAC file. The ASCII trace file contains one value per line,
which represents time series data. The arguments are the number of points
in the ASCII file (NPTS), and the time-increment separating the data values
(DELTA). The usage for writesac is:

  writesac NPTS DELTA < ASCIITraceFile_no_header > SACFile

  example: 
  writesac 1000 0.015 < file.asii > file.sac

getdtraces
----------
"getdtraces" reads a multiplexed (timestep ordered) traces output file (with 
header) and converts it into a demultiplexed (trace ordered) file (without 
header).  It also can be used to to extract a single trace from the multiplexed
traces output file and output as a SAC file (the trace number starts at 0). 
The usages are:

  getdtraces       < E3D_traces_output_file_(multiplex) > demultiplexed_traces
  getdtraces trace < E3D_traces_output_file > SACFile

  examples:
  getdtraces   < outH.0.TVx > outH
  getdtraces 5 < outH.0.TVx > sac5.x

Typing "getdtraces" with no arguments or input/output files will print
the usage message.

getrle
------
"getrle" is used to compress a gridded velocity file. The gridded velocity
file will have the format described above, in the section on "Velocity
Model Input". The usage for "getrle" is:

  getrle < unCompressedFile > RLECompressedFile

getunrle
--------
"getunrle" is used to uncompress a rle compressed gridded velocity file. It 
will convert the file back into it's original binary format, which is described
above in the section on "Velocity Model Input". The usage for "getunrle" is: 

  getunrle < RLECompressedFile > unCompressedFile

swapbytes
---------
"swapbytes" is used to convert binary data files (sac files, image files, 
volumes, gridded input files for the vfile option, etc.) generated on a
"little endian" architecture (e.g., PC's under Linux or DEC Alpha chips)
to ones appropriate for "big endian" architectures (e.g., Sun and SGI
workstations, most other Unix machines). It can also convert files from
"big endian" to "little endian". This would be done, for example, when
using a gridded velocity file created on a Sun Workstation to one that
will be used for a simulation on a DEC Alpha. "swapbytes" does simple
byte swapping and assumes words (integers and floats) are all 4-bytes.
The usage for "swapbytes" is:

  swapbytes [optional parameter] < littleEndianFile >    bigEndianFile
  swapbytes [optional parameter] <    bigEndianFile > littleEndianFile

The optional parameter is used to swap only the first N bytes in the file.
This would be done, for example, to concert a Sun Raster file created on
a "little endian" machine to a real Sun Raster files which is assumed to be
created on a "big endian" machine. Since a typical Sun Raster format
(from e3d) has a 32 byte header followed by the bytes represented in the
image (which should not be swapped), one would only want to convert the
first 32 bytes.

Example usages for swapbytes are:

  swapbytes    < DEC_Alpha_sacfile.000.x    > Sun_sacfile.000.x
  swapbytes 32 < DEC_Alpha_rasterfile.100.p > Sun_rasterfile.100.p

Note: swapbytes currently does things in core. Currently, the files
can only be 100 MBytes in size. This is an easily changeable parameter
in the code.