Sven Ivansson Feb. 23, 1998
National Defence Research Establishment (FOA)
S-17290 Stockholm, Sweden
email: sveni@sto.foa.se
General introduction to the program rpress (rpress = Reflectivity PRESSure)
------------------------------------------
The program rpress belongs to the rpress package. It can be used for cw
("continuous wave" with a monofrequency time dependence) computations as
well as computation of transients. The underlying theory is described in
the references. A range-independent (laterally homogeneous) fluid-solid
medium is assumed, with a fluid region interspaced between two (possibly
vanishing) solid regions. The lower and upper boundaries may independently
be either free, rigid, or of "homogeneous half-space" type (either fluid or
solid). The wave-field is created by a symmetric point source in the
(interior) fluid region, and the receivers are also located there. Hence,
backward propagation through solid regions will not be needed. (To get the
full wave-field, the program rpressfw should be used.)
Units are s for time, km for length and "10**12 kg" for mass. However, dB
is "re 1 m" in the way explained below. The output concerns the pressure
field in the fluid region given in "(10**15 N)/(km**2)".
For the cw computations, an underlying harmonic time dependence s(t)
according to "s(t)=exp(-i*omeg*t)" is understood, omeg>0 being the angular
frequency and t being time. For computation of transients, the source
wavelet s(t) can be arbitrarily specified. In any case, the point source
has a diagonal moment tensor with all three diagonal elements equal to
" s(t) " "(10**15 N)*(km)". The output for the cw computations concerns
the complex pressure field, and phases are given in addition to the
amplitudes.
For typical applications, the source strength as specified would be
unrealistically large, and one would have to scale down the results. However,
the output in dB "re 1 m" would of course not be affected.
As specified in the information file read.me in this directory rpress, there
are other programs in the rpress package. Some important general remarks on
the implementation are also given there.
Input to the program rpress
---------------------------
LX LI(TOP) ICWORKDIM
FILM
MHQ0 FREQMHQ0 IAZIMI
if ( iazimi not 0 ) then
ALFA1AZIM
end if
ICUTTYP
if (icuttyp=3) then
AIMCUT
end if
EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX
EPSFSUB HSTEPFSUBPRE
PREEPS0 IXTYP0 NRUS0PRE IRUS0 JRUS0
if (lx<=1) then
FREQMIN FREQINT NFR
if ( freqmin<0.0 and freqint=0.0 ) then
OFREQ(1:nfr)
end if
else
TMIN DAUER URED
PREALA
IUROT
FREQDOM
if (freqdom<=0.0) then
CR9FIL
end if
end if
RMIN RINT
if ( rmin<0.0 and rint=0.0 ) then
OR(1:li)
end if
FILO
if (lx>1) then
MAINTR
end if
FILG ! filg=' ' allowed only when lx>1
if ( filg not ' ' ) then
IDB
MAINAM
MAINPH
end if
INTREST
IFILONPRE
BESSFIL
IADINTXFI
if ( iadintxfi not 0 ) then
SMAXPRE NELMAXPRE NBACKMAX
IXTYP NRUS IRUS JRUS
ARGDIFLIM1 ARGDIFLIM0 ARGDIFLIM2
end if
NPARTS
do jp = 1, nparts
IBESSARR(jp) NPUARR(jp)
if ( npuarr(jp) > 0 ) then
UPHOERN(1:npuarr(jp),jp)
else if ( npuarr(jp) = 0 ) then
MAINMOD(jp)
XYONENORM(jp)
else if ( npuarr(jp) < 0 ) then
AIMEND(jp)
end if
end do
if (nparts>0) then
HMAXDH HMAXRELWP
if ( iadintxfi not 0 ) then
EPSH01Q PREEPS
HKVARPART LUNKVAR FAKTKVAR
if ( hkvarpart>0.0 and faktkvar>0.0 ) then
DDSUPMIN DDSUPMAX DDSUPFAKTMAX
DERIVDRELMAX DRELMAX DRELGOALFAKT
ARGDIFSTORE
end if
end if
end if
ICONTTYP
if ( iconttyp not 0 ) then
ITAIL(1) ITAIL(2)
if ( itail(1) not 0 ) then
IBESSTAIL(1) ARGRKTAIL(1) UP0TAIL(1)
end if
if ( itail(2) not 0 ) then
IBESSTAIL(2) ARGRKTAIL(2) UP0TAIL(2)
end if
ANPHPMINTERM
if ( iadintxfi not 0 ) then
PREEPSTERM
end if
if (iconttyp=1) then
IRCONT PREEPSCONT
else if (iconttyp=2) then
KMAXEU PREEPSEU
end if
end if
Explanation of the parameters
-----------------------------
LX is less than or equal to one for cw ("continuous wave" with a mono-
frequency time dependence) computations whereas it is greater than one for
computation of transients. In the latter case, LX specifies the number of
time samples to be computed and it should be an integer power of two.
Transients are computed by frequency synthesis. Hence, cw output can be
obtained here as well.
LI is the number of vertical receiver arrays, i.e. receiver ranges.
ICWORKDIM is a parameter to give the dimension of certain work-space
arrays.
FILM, character*32, is the name of the infile that specifies the
range-independent medium. The medium is built up by a number of layers,
each of which is divided into a number of sublayers. If so is indicated in
the infile FILM, the number of sublayers for a given layer will be
multiplied by MHQ0*max(1,f/FREQMHQ0) where f is the current frequency.
IAZIMI is nonzero to specify frequency dependence of the characteristic
medium slownesses (velocities) according to Azimi's dispersion law (e.g.
Aki and Richards 1980 p. 176). (In the presence of anelastic attenuation,
causality for transients requires some kind of frequency dependence.) The
parameter ALFA1AZIM is a time constant given in s, see (5.76) in the book
by Aki and Richards.
ICUTTYP specifies the type of branch-cuts to be used if the lower or upper
boundary of the medium is a homogeneous half-space. The values 1 and 2 will
imply hyperbolic and vertical branch-cuts, respectively. (The value 3
will imply a kind of mixture, with a change of type where the imaginary
part of horizontal slowness is +-AIMCUT s/km.)
The parameters
EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX
EPSFSUB HSTEPFSUBPRE
are relevant only when numerical integration, as specified in the infile
FILM, is to take place through certain layers. This option has not yet
been fully tested, however, and we omit the detailed presentation of these
parameters.
The field results (for each frequency) can be computed with error control,
implying that the computations can be refined adaptively by extrapolation
until a specified accuracy is achieved. There are two levels of
extrapolation: one "inner" for the 'horizontal slowness' integrals and one
"outer" concerning refinement of the medium model, see section 5.1 in
Ref. 0. For some details on the (inner) extrapolation process, see
section II in Ref. 1.
PREEPS0, IXTYP0, NRUS0PRE, IRUS0 and JRUS0 are control parameters for the
outer extrapolation. (For many medium models, however, outer extrapolation
will not be needed.)
PREEPS0 is the specified relative accuracy. For a vertical receiver array
at the range r km, it is transferred to an absolute accuracy requirement
given by eps0=PREEPS0*pm0/r. (The reference pressure pm0 will be introduced
below in connection with the parameter IDB. For r equal to zero, eps0 is
set to PREEPS0*pm0.)
IXTYP0 specifies the type of extrapolation to be used. The values
"1,-1,0 / 2,-2 / 3,-3" will imply "polynomial / rational /
inverse-polynomial" extrapolation, respectively. Within each of these
groups, the different numbers are used to specify different ways of
estimating errors in the extrapolation tables (see p. 1571 in Ref. 1). The
choice in this respect is not very important in general and we omit a
detailed description.
NRUS0PRE gives the depth of the extrapolation tables (typically about 8).
A minus sign will imply that log output for the outer extrapolation is
provided in outfiles ftn100, ftn101 and ftn102.
IRUS0 is the minimum depth after which extrapolation is allowed to
terminate (typically 2).
JRUS0 is the maximum number of columns in the extrapolation tables
(typically about 7).
For cw computations, NFR frequencies OFREQ(1:nfr) Hz are considered in
sequence. They can be given explicitly, or they can be obtained as
OFREQ(jf)=FREQMIN+jf*(FREQINT/NFR) for jf=1,..,nfr.
For computation of transients, the nominal start time is TMIN s and the
width of the time window is DAUER s. In particular, the time between
successive samples will be DAUER/LX s. URED is the reduction slowness
in s/km. Hence, for a receiver at the range r km from the source, the
actual start time for the output trace will be TMIN+r*URED s.
PREALA is a "temporary reduction factor during DAUER". For PREALA
less than one, complex frequencies will be used during the computations
in order to mitigate the effects of fft wrap-around. Recommended values
are PREALA=0.01 or 0.001.
IUROT is used to indicate whether or not given values for horizontal slowness
are to be rotated in the complex plane because of nonreal frequencies. IUROT
is irrelevant when PREALA=1.0. Rotation is recommended, it is indicated by
IUROT=1 as opposed to IUROT=0. After rotation, a given real horizontal
slowness will correspond to a real horizontal wavenumber.
When FREQDOM is positive, it is related to the dominant frequency of the
source in Hz. The spectrum for the moment-tensor wavelet s(t) of the source
is 0.5*(1+cos(pi*0.5*f/freqdom)) for abs(f) < 2*freqdom and zero elsewhere.
Hence, a kind of weighted sinc-function sum is obtained in the time domain.
The time function is not strictly causal, but it is translated in time so
that it appears to start at time zero. (It can be plotted using a program
called sigrpress.)
If FREQDOM is given as a nonpositive number, however, the source spectrum
is taken from the file CR9FIL. (A program called cr9rpress can be used to
compute the required spectrum of a desired wavelet given in the time
domain.)
Equivalent "freqmin freqint nfr" are written on the terminal. This can be
useful for initial runs with the programs rmodfaut and rmodpfw if mode
contributions are to be included.
The LI vertical receiver arrays are located at ranges OR(1:li) km from the
source. These ranges can be given explicitly, or they can be obtained as
OR(ir)=RMIN+ir*(RINT/LI) for ir=1,..,li. In any case, they must form an
increasing sequence of nonnegative values.
FILO, character*32, is the name of an outfile which essentially contains
a listing of the input parameters.
For computation of transients, the traces will appear in outfiles with
prefix MAINTR (character*32). There will be one outfile for each receiver
depth as specified in the infile FILM, the extensions will simply be
1,2,3... Each such outfile will contain LI traces, one for each range.
FILG, character*32, is the name of an outfile which contains the receiver
geometry (depths and ranges in km).
IDB specifies the type of output that is to be given for the cw point-source
results. Consider our point source located in a homogeneous acoustic (no
shear) whole-space with the same medium parameters as at the source depth
in our original medium. Here our point source will require a certain power
(which is in general different from the power required in our original
medium). Next, consider a point source of the same frequency and power but
in a homogeneous acoustic (no shear) whole-space with velocity 1.5 km/s, no
anelastic attenuation and density "10**12 kg"/(km**3) (i.e.
1.0 kg/(dm**3)). Let "pm0" be the corresponding peak pressure, in
(10**15 N)/(km**2), at a distance of 1 km. For IDB>=0, the field results will
be given relative "pm01m=1000*pm0" and for IDB>0 in dB. For IDB =-1,-2 they
will not be scaled in this way, IDB=-1 will imply output in single precision
whereas IDB=-2 will imply double-precision output. For IDB=-3, the complex
field results will be multiplied with the value of the Lame parameter lambda
at the matching depth (as given in the infile FILM, see below). This is
useful for tests of reciprocity.
MAINAM is the prefix for the field-amplitude outfiles. There will be one
file for each frequency, the frequency will be given as extension. (For the
case when the frequency is not an integer in Hz, the extension will be the
integer part followed by X and the sequence order of the frequency. For
example, if MAINAM is the string outampl and the third frequency to be
considered is 25.3 Hz, the outfile for this frequency will be outampl.25X3.)
In the same way, MAINPH is the outfile prefix for the phase values of the
complex field. The phases will be given in degrees.
A nonzero INTREST will imply an analytical computation of part of the
field, as described at the end of section 6.3 in Paper IV of Ref. 0. This
option should only be used for computations of the whole field without
"mode parts" and "branch parts" (positive NPUARR() as described below),
the integrands will be reduced by the correct amount.
IFILONPRE indicates whether or not the Filon rule will be used for the
primary trapezoidal integral estimates. The value 1 will imply Filon-rule
estimates, whereas the value 0 will imply the ordinary trapezoidal rule.
(Some other options are available as well, but they are not normally
used.)
BESSFIL gives the name of an infile (created using a program called
besseltab) containing asymptotic approximations of a specified order to
certain Bessel and Hankel functions. Normally, however, the character % is
given, indicating an option to use standard first-order asymptotic
approximations.
IADINTXFI controls the type of numerical integration to be performed. A
nonzero value will imply adaptive integration with error control, whereas
the value zero will imply trapezoidal integration with a fixed step-size.
For positive IADINTXFI, log output from the adaptive integration process
will be provided on a ftn??? file with ??? replaced by the number IADINTXFI.
(The numbers 6 and 66 will indicate the terminal, however.)
Integration-step histograms as given in Fig. 5 (bottom frame) of Ref. 1 are
written on ftn"???+1" provided IADINTXFI is 6 or greater than or equal to 80.
Numerical integration takes place along a number of straight paths in
complex 'horizontal slowness' space. For adaptive integration with error
control (i.e. nonzero IADINTXFI), a number of control parameters are needed.
The parameters we have in mind are
SMAXPRE NELMAXPRE NBACKMAX
IXTYP NRUS IRUS JRUS
ARGDIFLIM1 ARGDIFLIM0 ARGDIFLIM2.
A negative value of SMAXPRE will imply the corresponding positive value as a
pre-set upper limit for the maximum allowed length of the "planned path
segments" into which each straight path is divided. A nonnegative value of
SMAXPRE will give a computed upper limit, this is the normal case.
A positive value of NELMAXPRE gives the maximum allowed number of
simultaneously active subintervals obtained by interval splitting for
a planned path segment. Zero or a negative value will result in a
computed maximum number, this is the normal case.
NBACKMAX is the maximum number of back steps with wasted evaluations
of the integrands for each planned path segment. (Normally, each
evaluation of the integrands is stored in a stack. If the allocated
space for the stack becomes occupied, the computations can nevertheless
be continued provided a "back step" is taken and a new effort is made
with a smaller planned path segment. The computations will be stopped
when the number of back steps exceeds NBACKMAX.)
IXTYP, NRUS, IRUS, and JRUS control the inner extrapolation process
for computation of the integral along each planned path segment. They
correspond in the obvious way to the parameters IXTYP0, NRUS0PRE,
IRUS0 and JRUS0 for the outer extrapolation. (No minus sign for NRUS
is allowed, however.)
ARGDIFLIM1, ARGDIFLIM0 and ARGDIFLIM2 are certain limits for the maximal
argument differences in degrees of the dispersion function between
consecutive points used at a particular extrapolation depth for a particular
subinterval used in the process of adaptive integration. One should
have 0.0 <= argdiflim1 < argdiflim0 < argdiflim2. These parameters can be
used to bring about an early subinterval splitting (which can be
economical). The logic is adapted to cases of simple poles of the
integrands, i.e. zeros of the dispersion function, close to the
integration path. Early subinterval splitting will be enforced if the
maximal argument differences in both halves of the subinterval exceed
ARGDIFLIM0 and also if the maximal argument difference in one half exceeds
ARGDIFLIM2 whereas no argument difference in the other half exceeds
ARGDIFLIM1.
In the computations, the field is built up from different contributions or
parts. There will be "mode-path-branch parts" and "tail integrals".
NPARTS is the number of path-mode-branch parts. For each such part, values
must be assigned to some parameters.
IBESSARR() specifies the type of Bessel/Hankel function to be used.
The value 1 or 2 will prescribe J0 whereas -1 or -2 will prescribe
0.5*H0(1). The tables from BESSFIL will be used for the values +-1,
whereas +-2 will imply very accurate computations using subroutines
provided by Anders Bostrom. The value -3 for IBESSARR() is a special
option implying quick computations of 0.5*H0(1) using first-order
asymptotics and simple multiplications to move out in range. This option
can only be used for vertical receiver arrays that are uniformly
distributed in range.
NPUARR() is positive to indicate a "path part", zero to indicate a "mode
part", and negative to indicate a "branch part".
When NPUARR() is positive, there is a path part and NPUARR() will give the
number of corner points for the corresponding polygonal integration path.
The complex horizontal slowness coordinates for the corner points, in their
correct order, will be contained in UPHOERN(,). The P and S branch-points of
a lower homogeneous half-space are specified by +-991 and +-992,
respectively. Here, the plus sign indicates the first quadrant and the
minus sign indicates the third quadrant.
When NPUARR() is zero, there is a mode part and the string MAINMOD(),
character*32, is the prefix for infiles containing modal horizontal
slownesses and related information. The extensions are the frequencies
in the same way as for the outfile prefix MAINAM. When XYONENORM() is
nonzero, these infiles are obtained as the "mainmodu " outfiles from the
program rmodfaut. When XYONENORM() vanishes, they are obtained as the
"mainmodp " outfiles from the program rmodpfw run with the "iuttyp=0" option.
The mode parts (residue contributions) can be computed in three different
ways, depending on the value of XYONENORM(). In each case, a positive
orientation of a surrounding integration path is assumed.
-> When XYONENORM() is positive, the mode parts are computed by integration
along a circular path around each modal horizontal slowness as specified in
the MAINMOD() infile, see Appendix B of Paper IV in Ref. 0. The planned
radius in s/km of the integration circle in the complex horizontal slowness
plane is given by XYONENORM(). However, for each particular mode, the
radius may be modified in order to ascertain that one and only one modal
horizontal slowness is encircled.
-> When XYONENORM() is zero, the appropriately excited mode parts as
functions of depth are directly read from the MAINMOD() infile. Very
little remains to be done, essentially only scaling with the appropriate
range-dependent Hankel function.
-> When XYONENORM() is negative, the mode parts are computed according to
Theorem 1 of Ref. 7. This option should be used with caution, however, since
the particular choice of matching depth can be critical. (In fact, different
matching depths may be needed for different modes.)
When NPUARR() is negative, finally, there is a branch part and branch-cut
integrals will be computed. The value -1 indicates the P-cut of a lower
homogeneous half-space, -2 indicates the S-cut and -12 indicates both cuts.
(It is required that the P- and S-cuts are separated. Branch-cut integrals
for an upper homogeneous half-space have not been implemented.) The
integration will be performed along the cut(s), away from the real
horizontal slowness axis, until the imaginary horizontal slowness part
becomes AIMEND() s/km (AIMEND() is assumed to be positive). A change of
integration variables is made, and this is essential close to a
branch-point (see Appendix A of Paper IV in Ref. 0).
HMAXDH and HMAXRELWP determine the maximum allowed (final) integration
step-size in s/km along the integration path. This step-size is obtained
as the minimum of HMAXDH and HMAXRELWP/(f*OR(LI)), where f is the current
frequency in Hz. Hence, HMAXRELWP provides a dimension-less specification
in terms of the minimum asymptotic Hankel-function period. (When OR(LI) is
zero, HMAXRELWP will be irrelevant and HMAXDH will be taken.) To be
cautious, cf. the Appendix to Ref. 1, HMAXRELWP should be taken somewhat
smaller than 1.0.
When adaptive integration with error control is performed, the available
error tolerance is uniformly distributed along the integration path.
However, even for a very small subinterval, over which integration is
to be performed, the share of the available error tolerance will be at least
EPSH01Q. The example in Remark 4 (ii) of section II in Ref. 1 corresponds to
setting EPSH01Q to 0.01.
PREEPS is a relative error tolerance for adaptive integration. It
corresponds in the obvious way to the parameter PREEPS0 for the outer
extrapolation process.
If HKVARPART is positive, then the option for adaptively changing the
integration path will be put into action; see the final paragraph of
section III in Ref. 1. When the current integration step-size for the
trapezoidal-sum estimates falls below HKVARPART times the "maximum
allowed (final) integration step-size" (as determined from HMAXDH and
HMAXRELWP), a path segment along the positive real axis will be moved into
the fourth quadrant. If positive, HKVARPART might typically be 0.01.
The absolute value of FAKTKVAR (typically 1.0) determines the relation
between the depth of the roundabout path and the width of the original path
segment. A minus sign for FAKTKVAR will indicate that the search for poles
of the integrand within the enclosed region be omitted. The path segments
to be moved will be temporarily stored on the file ftn???, where ??? is to
be replaced by the number LUNKVAR.
However, particular caution is needed when computing the field in a medium
without anelastic attenuation but with some solid region (Ivansson and
Karasalo, 1998, manuscript in preparation).
A number of additional parameters are required if the check for poles of the
integrand (appearing among the zeros of the dispersion function) is to be
performed. This check is made using the argument principle for analytic
functions (see Ref. 2) to determine the number of enclosed zeros. The
parameters in question are DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,
DRELMAX,DRELGOALFAKT and ARGDIFSTORE. Since these parameters are taken from
the program rmodfaut, we refer to the information file rmodfaut.txt (in
this directory rpress) for the explanation.
As mentioned, the computed field is built up from "mode-path-branch parts"
and "tail integrals". Tail integrals will be present if ICONTTYP is nonzero.
They can be computed in two different ways, depending on whether ICONTTYP is
1 or 2.
There can be two tails, the first one (present if ITAIL(1) is nonzero) in
the first quadrant and the second one (present if ITAIL(2) is nonzero) in
the second quadrant; see equation (16) and Fig. 3 in Ref. 1. A nonzero
ITAIL() can be equal to 1 or 2. For ITAIL() equal to 1, an asymptotic form
of the Green's function will be used. When ITAIL() equals 2, the Green's
function will be accurately computed. Each tail is specified by
IBESSTAIL(), ARGRKTAIL() and UP0TAIL().
IBESSTAIL() specifies the type of Bessel/Hankel function to be used, in the
same way as IBESSARR() for the mode-path-branch parts.
ARGRKTAIL() gives the direction of the tail path, in degrees relative to
the direction of the positive real horizontal slowness axis. Typical values
are ARGRKTAIL(1)=+30.0 degrees and ARGRKTAIL(2)=-30.0 degrees. (Note:
Compared to equation (16) and Fig. 3 in Ref. 1, there is a sign change and a
change of direction for the second tail.)
UP0TAIL(1) gives the start point of the first tail, whereas UP0TAIL(2) gives
the end point of the second tail. Both these parameters are complex
horizontal slowness values in s/km.
The tail integrals are computed successively along asymptotic
quarter-period segments for the Hankel functions. The maximum allowed
(final) integration step-size in s/km along the integration path is the
length of the current quarter-period segment divided by ANPHPMINTERM
(typically 7.0). For adaptive integration, PREEPSTERM determines the
accuracy requirements for each quarter-period computation in the same
way as does PREEPS for the mode-path-branch parts.
For ICONTTYP equal to 1, the chosen quarter-periods are those for the
IRCONT:th vertical receiver array. (Typically, IRCONT is set to 1.) In this
case, the tail computations will continue with quarter-period segments
farther and farther away from UP0TAIL() until the contribution from the last
quarter-period falls below a limit determined from PREEPSCONT in the same
way as accuracy requirements are determined from PREEPSTERM.
For ICONTTYP equal to 2, separate tail computations are performed for each
vertical receiver array. (This will slow down the computations for large
values of LI.) In this case, Euler's transformation is used for convergence
acceleration (e.g. Dahlquist, Bjorck, Anderson 1974, Problem 3 in section
7.6). KMAXEU will denote the maximum number of columns in the "extrapolation
table" whereas PREEPSEU will determine the error tolerance in the same way
as accuracy requirements are determined from PREEPS.
Format for the medium infile FILM
---------------------------------
ZSOURC0 ZSOURC1 NS NSHC
0
AA BB RRHO DBLA DBLB
a number of data blocks of the following type
ZZ IIFSMTYP IIMHQANI
AA BB RRHO DBLA DBLB NSH NH
Explanation of the format for the medium infile FILM
----------------------------------------------------
The infile FILM is read by the subroutine rmodel, and the parameters of
the fluid-solid medium are set up along with the matching and receiver
depths. Concerning the matching depth, see section 5 of Ref. 6. For reasons
of simplicity, it is restricted in the program to be in the (interior)
fluid region (see below).
For rpress, the source depth will be at the matching depth. Hence, it will
be in the (interior) fluid region, and the outgoing pressure field in the
vicinity of the source appears as
-exp(i*omeg*R/alpha) * (omeg**2)/(4*pi*(alpha**2)*R)
where omeg is the angular frequency and R is the spherical distance to the
source, while alpha is the medium velocity at the source. Concerning
computation of transients, the factor omeg**2 will of course imply a
sign-changed second derivative of the primarily given moment-tensor
wavelet s(t).
The initial "AA BB RRHO DBLA DBLB" data concern the upper boundary. A
free or rigid boundary is specified by a vanishing or negative RRHO,
respectively. An upper homogeneous half-space is specified by a positive
RRHO, and the medium parameters are obtained from the data in the same
way as described described below.
ZZ is the depth to which the medium parameters in a data block refer.
The ZZ value of a subsequent data block must be greater than (or equal to)
that of the preceding data block.
IIFSMTYP should be zero. (A nonzero value is possible to indicate a layer
with continuous parameter variation, but we do not describe that option
here.)
IIMHQANI should be zero or one, the value one indicating
frequency-dependent division into sublayers using the parameters MHQ0
and FREQMHQ0. (Preparation for transversely isotropic layers has been
made, using other values of IIMHQANI as an indicator.)
AA is the P-velocity and BB is the S-velocity. A vanishing value of BB
indicates fluid.
RRHO is the density.
DBLA and DBLB are values for absorption (anelastic attenuation) of P- and
S-waves, respectively, given in "dB per wavelength".
NSH specifies the primary number of layers into which the depth interval
starting from the ZZ value of the previous data block shall be divided.
(As mentioned, further division may take place if IIMHQANI is one.) A
vanishing NSH is used to indicate a discontinuity with a "new start". In
the fluid case (i.e. when BB vanishes), NSH can be negative which indicates
a nonhomogeneous region starting from the ZZ value of the previous data
block. In this nonhomogeneous region, the square of the inverse P-velocity
will be a linear function of depth. When abs(NSH) is greater than one,
division into homogeneous (or nonhomogeneous with the linear inverse-square
P-velocity) layers is done by interpolation for the medium parameters.
NH is the number of receiver depths within the depth interval starting
from the ZZ value of the preceding data block. If NH does not vanish, NSH
must be divisible by NH. The receiver depths will be written in the outfile
FILG. For rpress, receivers are only allowed in the (interior) fluid region.
The last data block specifies the lower boundary. There is a lower
homogeneous half-space if RRHO is positive. If RRHO of the last data block
is negative or vanishes, on the other hand, the lower boundary will be
rigid or free, respectively.
The matching depth is specified by the parameters ZSOURC0, ZSOURC1, NS and
NSHC. ZSOURC0 and ZSOURC1 must coincide with the ZZ values of two
consecutive data blocks for the fluid region. NSHC must be equal to
abs(NSH) of the latter data block. NS, finally, can take the values
0,1,..,NSHC and indicates precisely at which interface the matching depth
is to be placed.
Format for the outfiles
-----------------------
Let NHYD be the number of receiver depths (i.e. the sum of the NH of the
infile FILM). The output in a file with prefix MAINAM or MAINPH is given as
a matrix of size NHYD x LI which is written row by row. Before the matrix,
however, there appear several lines (their number given as the first integer
in the initial line of the file) with repeated input data and an additional
line where the integers NHYD and LI are written.
Outfiles may also be created by the subroutine utspec. The idea is that the
user could edit this routine to get his own desired output format.
Concerning computation of transients, the outfiles with prefix MAINTR are
"unformatted" in the Fortran sense (in order to reduce their size). Their
structure can hopefully be understood from the Fortran code in the file
rpressmain.f in the directory rpress/cpm.
Subroutine directories
----------------------
The source code for the program and its subroutines can be found in the
following directories:
rpress/cpm (program and main routines)
rpress/cpm/single (routines that have different "multi" versions)
rpress/cpm/dhcomp (routines that have different "full-wave" versions)
rpress/int (general integration and extrapolation routines)
rpress/rut (general routines)
rpress/bess (Bessel- and Hankel-function routines)
rpress/slatec (slatec routines for complex Airy functions)
rpress/ik (routines to allocate space for arrays)
To conveniently produce the Green's function curves (essentially the
Hankel-transform integrand) in Fig. 23 of Paper IV in Ref. 0, trapezoidal
integration was used for simplicity (IADINTXFI set to zero) and special
output was written to disk. Normally, trapezoidal integration is performed
by the routine in the file rpress/int/trapezfi.f, but in this case the
modified routine in the file rpress/int/trapezfiSPEC.fORIG was used.
Example
-------
The following input to rpress was used for computing the pressure curves in
Fig. 1a of Ref. 4:
(Note: The dB output from rpress is "re 1 m", whereas the results in the
figure referred to were "re 1 km" and hence 60 dB higher.)
0 80 200000 ! lx,li(top),icworkdim
ref4fig1a.m ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
2 ! icuttyp
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsub
0.00001 2 0 2 5 ! preeps0, ixtyp0,nrus0pre,irus0,jrus0
1.0 80.0 80 ! freqmin,freqint,nfr
-1.0 0.0 ! rmin,rint
0.095 0.097 0.100 0.104 0.106
0.190 0.194 0.200 0.208 0.212
0.380 0.388 0.400 0.416 0.424
0.760 0.776 0.800 0.832 0.848
1.520 1.552 1.600 1.664 1.696
2.340 2.404 2.500 2.628 2.692
3.040 3.104 3.200 3.328 3.392
4.680 4.808 5.000 5.256 5.384
6.080 6.208 6.400 6.656 6.784
7.680 7.808 8.000 8.256 8.384
9.680 9.808 10.000 10.256 10.384
12.160 12.416 12.800 13.312 13.568
14.360 14.616 15.000 15.512 15.768
17.360 17.616 18.000 18.512 18.768
20.360 20.616 21.000 21.512 21.768
24.320 24.832 25.600 26.624 27.136
ref4fig1a.dat ! filo
ref4fig1a.g ! filg
1 ! idb
ref4fig1aa ! mainam
ref4fig1ap ! mainph
0 ! intrest
1 ! ifilonpre
% ! bessfil
-1 ! iadintxfi
1.0 0 8 ! smaxpre,nelmaxpre,nbackmax
2 9 2 7 ! ixtyp,nrus,irus,jrus
0.0 10.0 430.0 ! argdiflim1,argdiflim0,argdiflim2
1 ! nparts
2 2 ! ibessarr(),npuarr()
(0.0,0.0) (4.0,0.0)
0.01 0.95 ! hmaxdh,hmaxrelwp
0.01 0.000001 ! epsh01q,preeps
0.005 111 -1.0 ! hkvarpart,lunkvar,faktkvar
1 ! iconttyp
2 2 ! itail(1),itail(2)
-2 30.0 ( 4.0,0.0) ! ibesstail(1) argrktail(1) up0tail(1)
-2 -30.0 (-4.0,0.0) ! ibesstail(2) argrktail(2) up0tail(2)
7.0 ! anphpminterm
0.0000001 ! preepsterm
1 0.000001 ! ircont,preepscont
The medium infile ref4fig1a.m had the following appearance:
0.000 0.050 1 1 ! matching depth = source depth = 0.050 km
0
0.000 0.000 0.000 0.0 0.0
0.000 0 0
1.450 0.000 1.000 0.0 0.0 0 0
0.050 0 0
1.450 0.000 1.000 0.0 0.0 1 1
0.100 0 0
1.450 0.000 1.000 0.0 0.0 1 0
0.100 0 0
1.550 0.120 1.600 0.4 0.6 0 0
0.103 0 0
1.550 0.120 1.600 0.4 0.6 1 0
0.103 0 0
4.000 2.300 2.500 0.2 0.3 0 0
Note the specification of the integration path with consistent "uphoern(,)"
and "up0tail(1:2)"; cf. equation (16) and Fig. 3 in Ref. 1, although there
is a sign change and a change of direction for the second tail there. The
breakaway points +-4.0 s/km are not far enough out to include the Scholte
wave, which may have a modal slowness of more than 1/0.120 s/km in this
example. However, the Scholte wave should be negligible because the sediment
layer is thin and the common source-receiver depth is not close to the
water-sediment interface.
References
----------
See the information file read.me in this directory rpress.