grn_sph.f for computing greens function for synthetic seismograms using summed free oscillations for spherical earth models

click title to see Google Drive directory JParkCodes

grn_sph.f

The task of computing synthetic seismograms from summed free oscillations for an aspherical earth model is typically accomplished in two or three steps. The first step computes the hybrid free-oscillations for a chosen earth model, comprised of the "simple" spherical-earth modes whose particle motions follow single vector spherical harmonics. The second step takes those hybrid modes and computes greens functions for particle motion at a given source location (allowing six different time series, one for each component of the seismic moment tensor) and station location (for a given sensor and three time series for three components: vertical/horizontal). The last step is the computation of a synthetic seismogram by summing the six time series for each motion component by weighting with the components of the moment tensor.

Historically, I have used the greens function routine as the location for adding the instrument response, but I have reconsidered this choice. Computers are fast enough for the FFTs an instrument response to be computed on the fly in the last step. This also avoids the problem of updating the instrument response codes (notoriously unstable, owing to the various formats for instrument responses) in multiple greensfunction codes. My final code (sac_syndat.f) is used for all, and should be the place where instrument response either is (or is not) applied. The output of this routine has units of displacement in meters.

So, for making these codes more public, lets start with the simple code for synthetic seismograms on a spherical earth model.
The code is grn_sph.f, which has existed in one form or another for over 20 years (this explains the dead code inside it).

print *,'typical input for this program:'
print *,'0 20 --> freq interval to sum (mHz)'
print *,'/home/jjpark/Mineos/Models/premiss --> 1-D earth model'
print *,'/home/jjpark/Mineos/premissnl --> spheroidal info'
print *,'/home/jjpark/Mineos/premiss_eig --> spheroidal modes'
print *,'/home/jjpark/Mineos/premistnl --> toroidal info'
print *,'/home/jjpark/Mineos/premist_eig --> toroidal modes'
print *,'1000 --> # of overtone branches (here: all of them)'
print *,'0 --> # of data points to synthesize (0-> #=data: NPTS)'
print *,'/park/nobackup/Sumatra/VH/PET.VHZ --> data file'
print *,'.....'
print *,'. . . more data files'
print *,'.....'
print *,'stop'
print *,' '
print *,'The format of the greens fct files is impure SAC-format'
print *,'readable only by sac_syndat and related JPark codes'
print *,'It has a SAC header with one ASCII variable set to GRN'
print *,'and six*NPTS data values that follow the header'
print *,'NPTS for each component of the source moment tensor'
print *,'The greens-function file names are created from your input file'
print *,'by prepending gsph_ to the input filename'
print *,'and replacing the last character with z,r,t'

This code reads files full of free-oscillation displacement functions for a standard earth model, computes the strain of the mode at the desired source depth, and sums 6 traces for surface displacement at a specified distance from the source. Because the earth model is a spherical-earth model, only the source depth and source-receiver distance affect the waveforms that represent oscillatory displacement. The code, however, takes a data file to specify the source and receiver locations and creates synthetics for the geometry. A subsequent processing step (syndat.f) will take the greens function file, sum with the components of a moment tensor, and apply the instrument responses appropriate to the station's sensors and channel.

The start time for the earthquake determines where the decaying sinusoids begin. Older versions of green-fct codes used an explicit time that the user had to supply. This version assumes that the information is in the omarker time in the SAC header. There is also a bmarker time that is used in case the start time of the record does not correspond to the reference time in the data header.
c the reference time of the file could be before or after the first data point
c the ymdhmsec parameters define the reference time
c the origin quake of the quake and start of record are wrt reference time
origin=ahead(8) ! event start is the omarker in sac header
recstart=ahead(6) ! record start is the bmarker in sac header