delim $$
NAME
hspec91 - evaluate F(f,r) for 16 Green's functions
VERSION
(2.0) - New Formulation
(3.0) - Added Kjartansson (1979) constant Q operator and
wavefield separation.
AUTHOR
R. B. Herrmann, Saint Louis University, 1991
LANGUAGE
FORTRAN 77
LIBRARY REQUIRED
Access to the following nonstandard routines is required:
getarg(), numarg().
Computations are performed in complex double precision. The
following non-standard double precision complex routines are
required: cdsqrt(), cdabs(), dcmplx(), and cdlog().
SYNOPSIS
hspec91
[ -H]
[-A arg ]
[-F fref ]
[-Z]
[-R]
[-K]
[-SD]
[-SU]
[-SPUP]
[-SPDN]
[-SSUP]
[-SSDN]
[-RD]
[-RU]
[-RPUP]
[-RPDN]
[-RSUP]
[-RSDN]
[-?]
DESCRIPTION
The purpose of this program is to evaluate the medium
response for sixteen basic Green's functions as a function
of distance and frequency. The Green's functions are for a
point vertical downward force, a point horizontal force, a
point explosion, and double couple earthquake sources.
The meaning the the command line arguments are as follow:
-H Use the Hankel function for argument kr > arg.
-A arg
Argument for controlling the switch from Bessel to
Hankel function. (arg = 6.0 is the default).
-F fref
Set the reference frequency for causal Q. (fref = 1.0
is the default).
-Z
-R These flags are used to control the time of the first
sample of the resultant time series. The data input on
UNIT 05 permits the use of an initial time shift as
well as a time shift proportional to distance. If nei-
ther -Z nor -R are given, then the distance is the epi-
central distance. if -Z is given, then the distance is
the vertical source-receiver distance. If -R is given,
the distance is the hypotenuse between the source and
receiver locations.
[-K ]
Use Kjartansson Constant Q operator instead of the
Futterman operator (default = Futterman). Causal Q
means that the layer velocity is complex.
For Kjartansson the complex velocity is given by
v over { left [
( omega ^/^ omega sub r ) sup { - gamma } ^ left ( 1 ^-^ i
tan ( pi gamma ^/^ 2 ) right ) right ] }
$
where $ gamma = { 1 over pi } tan sup -
1 left ( 1 over Q right )
and $ v $ is the elastic velocity, $ omega $ the angu-
lar frequency and $ omega sub r $ the reference fre-
quency, e.g., $ omega sub r ^=^ 2 pi ^ fref $.
For Futterman, the complex velocity is
v ^ left [ {
1 ^+^ { 1 over { pi Q } } ln left ( { omega ^/^ omega
sub r } right ) ^+^ i over { 2 Q } } right ]
[-SD]
Perform wavefield separation at source to consider only
downgoing waves from the source. Default no separation
done at the source.
[-SU]
Perform wavefield separation at source to consider only
upgoing waves from the source. Default no separation
done at the source.
[-SPUP]
[-SPDN]
[-SSUP]
[-SSDN]
These permit wavefield separation at the source for
upgoing and downgoing P and upgoing and downgoing S,
respectively. The values are logically additive. The
-SUP is equivalent to -SPUP -SSUP together. Default no
separation done at the source. Due to the nature of the
wavefield at low frequencies mathematically correct
linear trends may be seen if the P-wave contribution is
desired independently of the corresponding S-wave.
[-RD]
Perform wavefield separation at receiver to consider
only downgoing waves from the source. Default no
separation done at the receiver.
[-RU]
Perform wavefield separation at receiver to consider
only upgoing waves from the source. Default no separa-
tion done at the receiver.
[-RPUP]
[-RPDN]
[-RSUP]
[-RSDN]
These permit wavefield separation at the receiver for
upgoing and downgoing P and upgoing and downgoing S,
respectively. The values are logically additive. The
-SUP is equivalent to -RPUP -RSUP together. Default no
separation done at the source. Due to the nature of the
wavefield at low frequencies mathematically correct
linear trends may be seen if the P-wave contribution is
desired independently of the corresponding S-wave.
Note that this may not be too useful for receivers at
the free surface since these options do not permit
separation of the wavefield for free surface motion for
incident P or incident S because of the difficulty of
separating the composite reflected P and S wavefields.
The use of these flags is a filtering operation applied
to the true motion. These receiver flags may have value
for receivers at depth, as for VSP.
[-?]
Do not run the program. Provide online help for command
line arguments.
INPUT (UNIT 05)
The control input is on UNIT LIN:
LINE 1
read(LIN,18)name4
18 format(a)
name4
- name of output file containing the F(f,r) spectra
LINE 2
19 format(2i5)
c-----
c read in time domain damping, sampling interval
c-----
read(LIN,*)alpha,dt
alpha
- time domain damping. The purpose of this is to reduce
the effect of periodicity of the Fast Fourier Transform
by working with a time series multiplied by e-at, where
t is measured from the beginning of the time series.
After computation, the effect of this operator is taken
out by multiplying by a factor e+at in the program
rhfoc91(VI). There is a tradeoff between eliminating
the wrap around and enhancing noise in the final step.
If the length of the time series is TMAX seconds, a
suitable value is to keep a TMAX = 2.5, which
corresponds to a factor of 10 at the end of the time
series.
dt - sampling interval is seconds
LINE 3
read(LIN,*)n,n1,n2
n - total number of samples in time series. This must be
a power of 2. This requirement is not checked in the
program.
n1
n2 - control for frequency limits of computation. Given an
n point time series, the frequency sampling interval is
df = 1 / ( n dt ), and frequencies corresponding to n1
and n2 are f1 = (n1 -1 ) df and f2 = (n2 - 1) df. The
Nyquist frequency is given by a value n/2 + 1. To get
all frequencies, set n1 = 1 and n2 = n/2 + 1.
LINE 4
read(LIN,*)ieqex
ieqex
- control for Green's functions computed. The choices
are as follow:
= 0 : earthquake and explosion
= 1 : point force and explosion
= 2 : earthquake, explosion and point force
= 3 : explosion only
= 4 : earthquake only
= 5 : force only
LINE 5
read(LIN,*)jbdry
jbdry
- parameter to control conditions at top and bottom
surfaces of the model. jbdry is of the form
jbdry = 10*surface + halfspace
surface = 0 - elastic
1 - free
2 - rigid
halfspace = 0 - elastic
1 - free
2 - rigid
e.g., jbdry = 11 consists of a model with both surfaces
free, and jbdry = 00 consists of a model inserted into
a wholespace, with the property of the upper halfspace
equal to that of the first layer in the model, and that
of the lower halfspace equal to that of the last layer
of the model.
LINE 6
read(LIN,*)mmax
mmax - number of layers in the model
LINE 7 (mmax lines)
read(LIN,*) (d(i),a(i),b(i),rho(i),qai),qb(i),i=1,mmax)
d - layer thickness - if d(mmax) is zero, the layer
thickness will be internal increased to permit sources
at deeper depths, however, this may not be what is
desired if the base boundary condition is not that of
an elastic halfspace.
a - compressional wave velocity
b - shear wave velocity
rho - layer density
qa - compressional wave quality factor, or its inverse,
for the layer. The program ASSUMES Q $ >=$ 1. If a
smaller value is input, it is ASSUMED to be inverse Q.
qb - shear wave quality factor, or its inverse, for the
layer. The program ASSUMES Q $ >=$ 1. If a smaller
value is input, it is ASSUMED to be inverse Q.
LINE 8
read(LIN,24) xleng, xfac
xleng
- control of wavenumber sampling following Bouchon
(1981). To avoid spatial periodicity, xleng must be at
least two times larger than any distance at which a
synthetic is computed and the distance traveled by the
fastest compressional velocity that arrives in the sig-
nal time window. Mathematically, if $ r $ is the epi-
central distance, $ z $ is the vertical separation
between source and receiver, and $L$ = xleng, then the
conditions for a wholespace (and also for layered
media) are
L ^>^ 2 r ( L ^-^ r ) sup 2 ^+^ z sup 2 ^>^ v sup 2 t
sup 2
xfac - control on upper limit of integration. Basically, the
integration over wavenumber will go up to xfac * w
/vmin , where vmin is the smaller of the lowest shear
wave velocity in solid layers or the lowest compres-
sional wave velocity in fluid layers.
LINE 9
read(LIN,*)mdpths
mdpths
- number of source depths for which synthetics are
desired
LINE 10 (mdpths lines)
do 2008 i=1,mdpths
read(LIN,*)depths(i)
2008 continue
depths(i)
- depth of source
LINE 11
read(LIN,*)mdpthr
mdpthr
- number of receiver depths for which synthetics are
desired
LINE 12 (mdpthr lines)
do 2009 i=1,mdpthr
read(LIN,*)depthr(i)
2009 continue
depthr(i)
- depth of receiver
LINE 13
read(LIN,*)ndist
ndist
- number of distances at which to compute synthetic
LINE 14 (ndist lines)
do 2010 i=1,ndist
read(LIN,1)rr,tshf,vr
r(i)=rr
tshift(i) = tshf
vred(i) = vr
2010 continue rr - distance at which synthetic is desired
tshift
vred - controls to set time of first output sample. If the
reduction velocity vred is zero, then the time of the
initial output sample is tshift. Otherwise, a reduced
travel time is used.
The initial time sample for a given distance i is given
by
t0 = tshift(i) + dist/vred(i) if vred(i) > 0, or
t0 = tshift(i) if vred == 0.0
The parameter dist = the epicentral distance r(i) if
neither the -Z nor the -R flags are given, dist =
abs(depths - depthr) if the -Z flag is given, and dist
= sqrt[ r(i)**2 + (depths - depthr)**2 ] if the -R flag
is given. This option is introduced to assist in making
VSP time histories.
LINE 15
read(LIN,*,end=4010,err=4010)cmax,c1,c2,cmin
cmax,c1,c2,cmin
- phase velocity limits for wavenumber filter. If
these are all < 0 or if this LINE is not in the input
file, NO wavenumber filtering is done. Wavenumber
filtering is a non-causal process and will lead to
numerical noise in the resultant time histories. Even
though cleaner seismograms arise when this is not done,
there are times when this is appropriate for computa-
tional efficiency. One case is in which there is a very
thin, low velocity material at the surface of the model
with the source not located in this layer. Rather than
integrate to very small wavenumber, a cutoff is
appropriate.
The wavenumber filter is set up as a cosine-tapered
bandpass filter. The output of the filter is [ 0,
cosine-taper, 1.0, cosine-taper, 0] for wavenumber in
the range [ 0, w / cmax, w / c1, w / c2, w / cmin, oo].
If c2 or cmin <= 0, then the filter is a high pass in
wavenumber with the taper at low wavenumber defined by
cmax and c1. If c1 or cmax <= 0, then the filter is
lowpass with the high wavenumber taper defined by c2
and cmin.
If wavenumber filtering is done, the trick of using a
known analytic form to improve the integral at short
distance is not used.
OUTPUT
Diagnostic output is placed on UNIT 06. The desired spectra
REFERENCES
Kjartansson, E. (1979). Constant Q-wave propagation and attenua-
tion, J. Geophys. Res. 84, 4737-4748.
SAMPLE DATA SETS
Three runs are made to present the results and to indicate
the limits of the program. The first run considers an elas-
tic medium problem with the receiver radial distance non-
zero, the second considers the same model with the radial
distance zero, and the next considers the fluid problem, and
the last consider wavefield separation.
ELASTIC WHOLESPACE, NON-ZERO DISTANCE:
The following shell script, controls the computations
through the plotting. The plots are given in Figure 1.
#!/bin/sh # ELASTIC WHOLESPACE cat > hspec91.dat <<EOF
swe0010
1.5600000E-01 0.1250000e+00
0128 1 065
2
00
1
0.40000+02 0.6000e+01 3.4640e+00 0.2800e+01 0.1000e-03
0.1000e-03
2.0000000E+02 0.3000000E+01 1 20 17 0.0
2.5 5.0 7.5 10.0 12.5 15.0
17.5 20.0 22.5 25.0 27.5 30.0
32.5 35.0 37.5 40.0 1
10.0 0.00 0.0 EOF hspec91 < hspec91.dat >
hspec91.out rhfoc91 -f swe0010 -p -l 2 > file91 prof91 <
file91 > prof91out mv plot* PLOTHSPC1 rm file91 rm swe0010
hspec91.dat prof91out hspec91.out
ELASTIC WHOLESPACE, ZERO RADIAL DISTANCE:
The following shell script, controls the computations
through the plotting. The plots are given in Figure 2.
#!/bin/sh #ELASTIC WHOLESPACE ZERO DISTANCE cat >
hspec91.dat <<EOF swe0010
1.5600000E-01 0.1250000e+00
0128 1 065
2
00
1
0.40000+02 0.6000e+01 3.4640e+00 0.2800e+01 0.1000e-03
0.1000e-03
2.0000000E+02 0.3000000E+01 1 20 17 0.0
2.5 5.0 7.5 10.0 12.5 15.0
17.5 20.0 22.5 25.0 27.5 30.0
32.5 35.0 37.5 40.0 1
00.0 0.00 0.0 EOF hspec91 < hspec91.dat >
hspec91.out rhfoc91 -f swe0010 -p -l 2 > file91 prof91 <
file91 > prof91out mv plot* PLOTHSPC2 rm file91 rm swe0010
hspec91.dat hspec91.out prof91out
FLUID LAYER OVERLYING ELASTIC LAYER; TOP IS HALFSPACE, BOTTOM IS
HALFSPACE:
The following shell script, controls the computations
through the plotting. The plots are given in Figure 3.
#!/bin/sh #FLUID TOP LAYER, NOTE MODEL REQUIRES SOLID LAYER
AT BASE cat > hspec91.dat <<EOF swe0010
1.5600000E-01 0.1250000e+00
0128 1 065
2
00
2
0.40000+02 0.6000e+01 0.0000e+01 0.2800e+01 0.1000e-03
0.1000e-03
0.10000+02 0.6000e+01 0.3464e+01 0.2800e+01 0.1000e-03
0.1000e-03
2.0000000E+02 0.3000000E+01 1 20.0 21 0.0
2.5 5.0 7.5 10.0 12.5 15.0
17.5 20.0 22.5 25.0 27.5 30.0
32.5 35.0 37.5 40.0 42.5 45.0
47.5 50.0 1
10.0 0.00 0.0 EOF hspec91 < hspec91.dat >
hspec91.out rhfoc91 -f swe0010 -p -l 2 > file91 prof91 <
file91 > prof10out mv plot* PLOTHSPC3 rm file91 rm swe0010
hspec91.dat hspec91.out prof91out
The next two runs test the other options of the program,
e.g., have the free surface free, and the base rigid. The
direct arrivals should be identical to the wholespace solu-
tions, except when complicated by reflections.
ELASTIC MODEL:
Figure 4 shows the results of the following computations.
#!/bin/sh # TOP INTER FACE IS FREE, BASE IS RIGID cat >
hspec91.dat <<EOF swe0010
1.5600000E-01 0.1250000e+00
0128 1 065
2
12
1
0.40000+02 0.6000e+01 3.4640e+00 0.2800e+01 0.1000e-03
0.1000e-03
2.0000000E+02 0.3000000E+01 1 20 17 0.0
2.5 5.0 7.5 10.0 12.5 15.0
17.5 20.0 22.5 25.0 27.5 30.0
32.5 35.0 37.5 40.0 1
10.0 0.00 0.0 EOF hspec91 < hspec91.dat >
hspec91.out rhfoc91 -f swe0010 -p -l 2 > file91 prof91 <
file91 > prof91out mv plot* PLOTHSPC4 rm file91 rm swe0010
hspec91.dat hspec91.out prof91out
FLUID TOP LAYER.
Figure 5 shows the result of the following computations.
#!/bin/sh # TOP SURFACE IS FREE, BASE IS RIGID, FLUID TOP
LAYER, SOLID LOWER LAYER cat > hspec91.dat <<EOF swe0010
1.5600000E-01 0.1250000e+00
0128 1 065
2
12
2
0.40000+02 0.6000e+01 0.0000e+01 0.2800e+01 0.1000e-03
0.1000e-03
0.10000+02 0.6000e+01 0.3464e+01 0.2800e+01 0.1000e-03
0.1000e-03
2.0000000E+02 0.3000000E+01 1 20.0 21 0.0
2.5 5.0 7.5 10.0 12.5 15.0
17.5 20.0 22.5 25.0 27.5 30.0
32.5 35.0 37.5 40.0 42.5 45.0
47.5 50.0 1
10.0 0.00 0.0 EOF hspec91 < hspec91.dat >
hspec91.out rhfoc91 -f swe0010 -p -l 2 > file91 prof91 <
file91 > prof10out mv plot* PLOTHSPC5 rm file91 rm swe0010
hspec91.dat hspec91.out prof91out
ELASTIC MODEL WAVEFIELD SEPARATION AT RECEIVER.
Figure 6 uses a multilayered model to compute the complete
synthetic and the component of the wavefield separated at
the receiver. The source is at a shallow depth, above the
receivers.
cat > hspec91.dat <<EOF swe0010
0.7800000E-01 0.1250000e+00
0256 1 129
1
10
3
0.15000+02 0.6000e+01 3.4640e+00 0.2800e+01 0.1000e-03
0.1000e-03
0.15000+02 0.7000e+01 4.0000e+00 0.3000e+01 0.1000e-03
0.1000e-03
0.00000+02 0.8000e+01 4.7000e+00 0.3300e+01 0.1000e-03
0.1000e-03
3.0000000E+02 0.3000000E+01 1 1.0 16 2.5
5.0 7.5 10.0 12.5 15.0 17.5
20.0 22.5 25.0 27.5 30.0 32.5
35.0 37.5 40.0 1
10.0 0.00 0.0 EOF
hspec91 -RD < hspec91.dat > hspec91.outrd rhfoc91 -f swe0010
-p -l 2 > file91rd prof91 < file91rd > prof91out mv plot*
PLOTHSPC6rd
hspec91 -RU < hspec91.dat > hspec91.outru rhfoc91 -f swe0010
-p -l 2 > file91ru prof91 < file91ru > prof91out mv plot*
PLOTHSPC6ru
hspec91 < hspec91.dat > hspec91.outnr rhfoc91 -f swe0010
-p -l 2 > file91rr prof91 < file91rr > prof91out mv plot*
PLOTHSPC6
rm file91* rm swe0010 hspec91.dat prof91out
ELASTIC MODEL WAVEFIELD SEPARATION AT SOURCE.
Figure 7 uses a multilayered model to compute the complete
synthetic and the component of the wavefield separated at
the receiver. The receiver is at a shallow depth, above the
sources.
cat > hspec91.dat <<EOF swe0010
0.7800000E-01 0.1250000e+00
0256 1 129
1
10
3
0.15000+02 0.6000e+01 3.4640e+00 0.2800e+01 0.1000e-03
0.1000e-03
0.15000+02 0.7000e+01 4.0000e+00 0.3000e+01 0.1000e-03
0.1000e-03
0.00000+02 0.8000e+01 4.7000e+00 0.3300e+01 0.1000e-03
0.1000e-03
3.0000000E+02 0.3000000E+01 16 2.5 5.0 7.5
10.0 12.5 15.0 17.5 20.0 22.5
25.0 27.5 30.0 32.5 35.0 37.5
40.0 1 1.0 1
10.0 0.00 0.0 EOF hspec91 -SD < hspec91.dat >
hspec91.outsd rhfoc91 -f swe0010 -p -l 2 > file91sd prof91
-R < file91sd > prof91out mv plot* PLOTHSPC7sd
hspec91 -SU < hspec91.dat > hspec91.outsu rhfoc91 -f swe0010
-p -l 2 > file91su prof91 -R < file91su > prof91out mv plot*
PLOTHSPC7su
hspec91 < hspec91.dat > hspec91.out rhfoc91 -f swe0010
-p -l 2 > file91ss prof91 -R < file91ss > prof91out mv plot*
PLOTHSPC7
rm swe0010 hspec91.dat prof91out
Man(1) output converted with
man2html