Seems that there are 3 ways to do this:

1. Place an array of gauges across the area of interest (how the 2d code works, but 3d obviously requires significantly more gauges)
2. Store the displacement values in the q vector, using b4step to update the values
3. Do post-processing on the results with a python script

I'm leaning towards (2) instead of (1), but what do you think (in particular, @rjleveque and @smoe1)?

• move process_gauges.py to /src/python/seismic if it will be used for multiple 2d seismic examples
• delete plot_okada.py if the same functionality now exists in Fault.plot_okada(...) in /src/python/seismic/dtopotools_horiz_okada_and_1d.py
• remove "reload(dtopotools)" from python files if not needed

It looks like fdisc.f is a carry over from another code and should be removed. A lot of the code in process_gauges.py also appears in setplot_surface.py... so it may also be a candidate for removal.