Sven Ivansson Feb. 24, 1998
National Defence Research Establishment (FOA)
S-17290 Stockholm, Sweden
email: sveni@sto.foa.se
General introduction to the program rmodfaut
--------------------------------------------
The program rmodfaut belongs to the rpress package. It is intended for
computation of modal horizontal slownesses for fixed frequencies. The modal
horizontal slownesses are found as zeros of the analytic dispersion
function, using the argument principle of complex analysis (see Ref. 2).
Group-slownesses (inverse group-velocities) are computed as well. 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). In
addition, an entirely solid medium (without fluid) can be handled.
Units are s for time, km for length and "10**12 kg" for mass. In particular,
the unit for slowness is s/km.
An underlying harmonic time dependence "exp(-i*omeg*t)" is understood,
omeg>0 being the angular frequency and t being time.
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 rmodfaut
-----------------------------
FILM
MHQ0 FREQMHQ0 IAZIMI
if ( iazimi not 0 ) then
ALFA1AZIM
end if
ICUTTYP0
if (icuttyp0=3) then
AIMCUT
end if
ISHEETTYP0
EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX
EPSFSUB HSTEPFSUBPRE
XYTOL0 IXTYP0 NRUS0 IRUS0 JRUS0
if (nrus0>0) then
XYONENORM
end if
FREQMIN FREQINT NFR
if ( freqmin<0.0 and freqint=0.0 ) then
OFREQ(1:nfr)
end if
ALADAUER PREALA
IUROT
NRECT
do irect = 1, nrect
UPREAL1(irect) UPREAL2(irect) UPIMAG1(irect) UPIMAG2(irect)
end do
IOPTION
ICHECK
RIKTSORTANG
MAINMODU
MAINBOX
DDSUPMIN DDSUPMAX DDSUPFAKTMAX
DERIVDRELMAX DRELMAX DRELGOALFAKT
ARGDIFSTORE
SLIMLEFT
REDUCFAKT DDELT XYTOL
Explanation of the parameters
-----------------------------
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.
ICUTTYP0 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.) In the presence of an upper
homogeneous half-space, only the case of vertical branch-cuts has been
implemented.
ISHEETTYP0 specifies the type of Riemann sheet to be used if the lower or
upper boundary of the medium is a homogeneous half-space. The value 1,2,3,4
will imply the ++,-+,+-,-- sheet, respectively. (Here, the P-wave sheet is
given before the S-wave sheet.) In the presence of two homogeneous
half-spaces, the possibility of different sheets for the lower and upper
half-spaces has not been implemented. (This would be an almost trivial
extension, but at present the complete list of leaky modes will not be
obtained.)
One should note that internally, the computations are always done with
"icuttyp=2" (vertical cuts). For the case ICUTTYP0=1, repeated
computations may be needed (done automatically). An idea by Brian Maranda
(pers. comm.) to consider a dispersion-function product to get all modes
on all sheets simultaneously has not been implemented.
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.
Extrapolation is possible for refinement of the medium model, see section
5.1 of Ref. 0. XYTOL0, IXTYP0, NRUS0, IRUS0 and JRUS0 are control
parameters for the extrapolation. (For many medium models, however, this
extrapolation will not be needed.)
XYTOL0 is the specified final accuracy for the modal horizontal slownesses.
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.
NRUS0 gives the depth of the extrapolation tables (typically about 8).
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).
XYONENORM is given when NRUS0 is positive. It is, essentially, a proposal
for the half-side of the search rectangles for computation of each modal
horizontal slowness as the medium is refined.
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.
PREALA is a "temporary reduction factor during the time ALADAUER". For
PREALA less than one, complex frequencies will be used during the
computations.
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, which is the normal case
when running rmodfaut. Rotation is indicated by IUROT=1 as opposed to
IUROT=0. After rotation, a given real horizontal slowness will correspond
to a real horizontal wavenumber.
NRECT is the number of rectangles in the complex horizontal slowness plane,
within which modal horizontal slownesses are to be sought. The coordinates
of the vertices of the rectangles are obtained from UPREAL1(),UPREAL2(),
UPIMAG1(),UPIMAG2().
IOPTION can be 0, 1 or 2. The value 0 means that only the number of modes
(zeros of the dispersion function) are computed, whereas 1 and 2 means
that the modal horizontal slownesses will actually be computed. In the case
2, some log information will be provided in the file ftn099.
ICHECK can be 0 or 1. The value 1 implies that an additional argument
variation check will be made, around each found zero of the dispersion
function.
RIKTSORTANG gives the order in which the modal horizontal slownesses will be
written in the outfile(s). It is an "angle of direction", in degrees, in the
complex horizontal slowness plane. For example, -180.0 will imply that the
modal horizontal slownesses are sorted according to decreasing real part,
and -90.0 will imply that they are sorted according to decreasing imaginary
part.
MAINMODU is the prefix for the outfiles with the computed modal horizontal
slownesses. 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 MAINMODU is the string outmodu and
the third frequency to be considered is 25.3 Hz, the outfile for this
frequency will be outmodu.25X3.)
In the same way, MAINBOX is the outfile prefix for other outfiles. It may
happen that a certain modal horizontal slowness is not computed, for example
if there are two very close modal horizontal slownesses, closer than
SLIMLEFT (see below). Limits for the "nzleft" modal horizontal slownesses
which are left are given in these outfiles.
DDSUPMIN, DDSUPMAX and DDSUPFAKTMAX are control parameters for the
computation of argument variation of the dispersion function along a line
segment. Appropriate steps for the horizontal slowness are adjusted
adaptively (see Ref. 2).
DDSUPMIN (typically 0.0 s/km) is the minimal step that is handled when the
path is stepped through to determine the argument variation. If steps
shorter than DDSUPMIN are suggested, the computations are stopped, this
may indicate a zero right on the path.
DDSUPMAX (typically 0.05 s/km) is the maximal step. If a larger step is
suggested, the next step-size will nevertheless be DDSUPMAX.
DDSUPFAKTMAX (typically 5.0) is the maximal step-change factor for two
subsequent steps. If a larger step is suggested, the next step-size will
nevertheless be DDSUPFAKTMAX times the previous step-size.
DERIVDRELMAX, DRELMAX and DRELGOALFAKT are further control parameters for
the computation of the argument variation.
When the argument variation computations are started, an appropriate initial
step for horizontal slowness is found by making a difference approximation
to the derivative of the dispersion function with respect to horizontal
slowness.
DERIVDRELMAX (typically 0.12, always less than 1.0) is the maximum tolerable
relative error in this estimate of the derivative of the dispersion function
at the starting point. (The requirement for a small relative error is relaxed
for small derivatives. One would otherwise get problems for starting points
at the origin, where the derivative of the dispersion function vanishes.)
DRELMAX (typically 0.4, always less than 1.0) is the maximum allowed
relative change of the dispersion function in each step. (If it is exceeded,
a change to a smaller step will be made.)
DRELGOALFAKT (typically 0.75, always less than 1.0) is a factor for choosing
the next step-size. The largest step-size that is expected to fulfil the
step-size requirements is reduced by multiplication with DRELGOALFAKT.
ARGDIFSTORE (typically 150.0 degrees, always less than 180.0 degrees) is a
lower limit for dispersion-function argument difference in degrees between
points along a line segment at which the argument is stored. Storage is
needed to get the correct multiple of 360 degrees when a rectangle is
divided, cf. Ref. 2, but the memory requirements can be reduced if
ARGDIFSTORE is increased.
As explained in Ref. 2, each search rectangle is in general split until
there is at most one modal horizontal slowness in each part. However, this
splitting process is stopped prematurely for a subrectangle if the length of
its longest side falls below SLIMLEFT.
REDUCFAKT, DDELT and XYTOL are tolerance parameters for the final
determination of each modal horizontal slowness using the secant rule.
REDUCFAKT (less than one) is the required reduction factor in each step for
the dispersion function when the secant rule is applied.
DDELT is the 'horizontal slowness' distance-difference limit for choosing a
second point to initiate the secant rule. (The first point is chosen as the
mid-point of the particular subrectangle.) Since the required derivatives are
actually available at little extra computational work, Newton's method could
be used as an alternative.
XYTOL is the desired accuracy for the modal horizontal slownesses.
Format for the medium infile FILM
---------------------------------
ZSOURC0 ZSOURC1 NS NSHC
NSOURCFW
if (nsourcfw>0) then
do i = 1, nsourcfw
NSOFW0(i) AMSOURCFW(i)
end do
end if
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
----------------------------------------------------
When NSOURCFW is positive, the file contains data for a vertical source
array (see the information file rpressfw.txt in this directory rpress).
NSOURCFW and NSOFW0() are integers, while AMSOURCFW() are complex numbers.
However, the source information is not relevant for running rmodfaut. The
receiver depths are not relevant either.
In other respects, things are as for a medium infile FILM to rpress, see
the information file rpress.txt in this directory rpress. There is an
additional option, however, to handle an entirely solid medium (without
interior fluid). This option is specified by NS=NSHC=0. In this case, one
must take ZSOURC0=ZSOURC1="the initial ZZ", and the matching depth will be
at this depth.
For reasons of computational efficiency, the number of layers should be
kept as small as possible. Splitting of homogeneous layers may be needed
in infiles FILM to the programs rpress and rpressfw to locate receivers
at desired depths. It can then be useful to have a version of FILM without
layer splitting when running rmodfaut for the same medium model.
For the same reason, a homogeneous half-space should "start" as early as
possible. (This may even be essential in order to make the problem
well-conditioned, since "leaky" boundary conditions far out will not be
stable after transportation to the depths of interest.) A word of caution
is needed here, however: If an rmodfaut MAINMODU outfile is to be used as
infile to rpress or rpressfw with the "npuarr()=0, xyonenorm()<0" option,
the boundary-condition depths (half-space starting depths) in the infiles
FILM to rmodfaut as well as rpress/rpressfw must coincide! The reason is
that derivatives of the dispersion function will be read from the MAINMODU
file when running rpress/rpressfw, hence the computed dispersion functions
must coincide.
Format for the outfiles
-----------------------
There will be four data sets in the outfiles with prefix MAINMODU (in
the case of extrapolation for medium refinement there will be a fifth
data set). Each data set starts with a line specifying the number of
lines (modes) with corresponding data that will follow.
The first data set contains the real and imaginary parts of the modal
horizontal slownesses (possibly rotated as specified by the parameter
IUROT). In addition, the type of Riemann sheet "relative vertical
branch-cuts" is provided for each mode. This information is needed when
modal field synthesis is subsequently performed by running rpress (or
rpressfw) with some "npuarr()=0 and corresponding xyonenorm() nonzero".
The second data set gives the estimated accuracy and the distance to the
border of the final subrectangle for each modal horizontal slowness.
The third data set gives the derivative of the dispersion function with
respect to (unrotated) horizontal slowness at each modal horizontal
slowness. In order to avoid overflow (and underflow), a real number "dlndu"
is given together with a complex number "derivdu"; the derivative will be
dexp(dlndu)*derivdu.
The fourth data set shows the (unrotated) complex group-slownesses, i.e.
derivatives of the modal wavenumbers with respect to angular frequency.
Possible frequency dependence of the characteristic medium slownesses, as
implied by layer splitting according to MHQ0,FREQMHQ0 or a nonzero IAZIMI,
is not taken into account here.
(The fifth data set gives the modal horizontal slownesses, possibly rotated,
as obtained from the extrapolation for medium refinement.)
Before the four (or five) data sets, 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 integer four (or
five) is written.
If "nzleft" is nonzero, one should in addition take a look at the
corresponding outfile with prefix MAINBOX.
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/slatec (slatec routines for complex Airy functions)
Example
-------
The following input to rmodfaut was used for computing the modal horizontal
slowness results in Fig. 3c of Ref. 2:
(Note: Horizontal wavenumbers in m**(-1) were shown in the figure, whereas
the rmodfaut output as it stands gives horizontal slowness in s/km.)
ref2fig3c.m ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
2 ! icuttyp0
1 ! isheettyp0
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsubpre
0.0 -2 0 2 7 ! xytol0, ixtyp0,nrus0,irus0,jrus0
25.0 0.0 1 ! freqmin,freqint,nfr
1.0 1.0 ! aladauer,preala
0 ! iurot
1 ! nrect
-2.0 4.0 -0.001 5.0 ! upreal1(),upreal2(),upimag1(),upimag2()
1 ! ioption
1 ! icheck
-180.0 ! riktsortang
ref2fig3c ! mainmodu
ref2fig3cx ! mainbox
0.0 0.05 5.0 ! ddsupmin,ddsupmax,ddsupfaktmax
0.12 0.4 0.75 ! derivdrelmax,drelmax,drelgoalfakt
150.0 ! argdifstore
0.0 ! slimleft
0.7 0.00001 0.000000001 ! reducfakt,ddelt,xytol
The medium infile ref2fig3c.m had the following appearance:
0.000 0.040 1 1 ! matching depth 0.040 km
0
0.000 0.000 0.0 0.0 0.0
0.000 0 0
1.4165 0.000 1.0 0.0 0.0 0 0
0.040 0 0
1.4165 0.000 1.0 0.0 0.0 1 0
0.040 0 0
1.650 0.390 1.3 0.165 1.85 0 0
0.063 0 0
1.650 0.390 1.3 0.165 1.85 1 0
0.063 0 0
2.000 0.500 2.6 20.0 25.0 0 0
0.263 0 0
2.000 0.500 2.6 20.0 25.0 1 0
0.263 0 0
2.000 0.500 -2.6 20.0 25.0 0 0
References
----------
See the information file read.me in this directory rpress.