Sven Ivansson Feb. 9, 1998
National Defence Research Establishment (FOA)
S-17290 Stockholm, Sweden
email: sveni@sto.foa.se
General introduction to the program rmodpfw
-------------------------------------------
The program rmodpfw belongs to the rpress package. It is intended for
computation of normal modes as functions of depth for fixed frequencies.
The indicator-function method of section 4.1 in Ref. 7 is used. Excitation
strengths for a given vertical source array can also be included, cf. Ref. 5
and Ref. 7. To find the modal horizontal slownesses, one should first run
the program rmodfaut. 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 rmodpfw
----------------------------
FILM
MHQ0 FREQMHQ0 IAZIMI
if ( iazimi not 0 ) then
ALFA1AZIM
end if
ICUTTYP
if (icuttyp=3) then
AIMCUT
end if
ISHEETTYP
EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX
EPSFSUB HSTEPFSUBPRE
FREQMIN FREQINT NFR
if (nfr>=0) then
if ( freqmin<0.0 and freqint=0.0 ) then
OFREQ(1:nfr)
end if
end if
ALADAUER PREALA
IUROT
MAINMODU
IUTTYP
if (iuttyp=0) then
KFWTYP
end if
MAINMODP
INOLX
IYSTAB
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.
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.)
ISHEETTYP 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.
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.
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.
However, there is a "special option", specified by a negative NFR, in which
the number of frequencies and the frequencies themselves will be read from
an infile (the infile MAINMODU, see below).
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 rmodpfw. 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.
MAINMODU is the prefix for infiles with 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.) These infiles are obtained as "mainmodu " outfiles
from the program rmodfaut.
In the "special option" characterized by a negative NFR, however, the
given MAINMODU string itself will be the name of an infile containing
frequency as well as 'horizontal slowness' information. There can be
several data sets in this infile, the "-NFR:th" will be taken. A good
way to study the variation of mode shape with frequency is to take NFR=-1
and MAINMODU as the "filfu" outfile from the program rmodfgr.
IUTTYP is used to control the output to files with name "MAINMODP plus a
suffix and an extension": 1 will provide displacement output in the files
with suffix _1 and _2, 2 will provide stress output in the files with
suffix _3 and _4, 3 will provide energy- and power-density output in the
files with suffix _7 and _8 and _9. One may combine 1 and 2 by writing 12,
and 1 2 and 3 by writing 123. The extension will specify the frequency,
in the same way as for the MAINMODU infiles.
In the "special option" characterized by a negative NFR, however, there
will not be any extension specifying frequency. Instead, the frequency
may vary among the modes in the same file.
In each of the MAINMODP outfiles, results for the pertinent number of
modes will be provided in sequence. For each mode, one line is given for
each interface depth with the depth value appearing first on the line. Note
that results are provided at each interface depth, in the solid as well as
fluid regions, regardless of whether or not a receiver has been placed
at that depth. In some cases two results are provided for the same depth,
for example for the horizontal displacement at a fluid-solid interface
where it may be discontinuous. For each mode, multiplication with the
appropriate complex factor is made to make the absolutely greatest element
of the computed v(z) vector (see Ref. 7) equal to 1.
It is also possible to have IUTTYP equal to zero (unless the "special
option" has been chosen).
When IUTTYP is zero, KFWTYP specifies the field component that is to be
provided in the outfiles "MAINMODP plus an extension specifying frequency
in the usual way". There are four options: 1 for horizontal displacement,
2 for vertical displacement, 3 for tangential stress in horizontal cuts,
4 for normal stress in horizontal cuts. The outfiles in question will give
the excited mode component scaled with the appropriate excitation factor
for subsequent field synthesis using the programs rpress and/or rpressfw
(with the "npuarr()=0,xyonenorm()=0.0" option there). The underlying
theoretical formula for the excitation factor is provided by Theorem 2 of
Ref. 7 (see equation (32) of Paper VI in Ref. 0). When IUTTYP=0, results
are only provided at the receiver depths, and the horizontal displacement
is left uncomputed within the fluid and set to zero there. The source
position(s) and strength(s) governing the excitation are given in the medium
file FILM as described below.
A nonzero INOLX will indicate that propagation through layers within a
lower homogeneous half-space will be performed in a special way,
analytically removing the undesired waves there.
IYSTAB concerns stabilization during the backward propagation. Stabilization
will only be performed if IYSTAB is nonzero. A negative value will imply
ordinary stabilization (section 5 of Ref. 6), whereas a positive value will
imply modified stabilization (Appendix A of Ref. 6, Appendix B of Ref. 7).
In the latter case, which is strongly recommended, the modified technique
is also used for stable computation of the "energy integral".
Format for the medium infile FILM
---------------------------------
ZSOURC0 ZSOURC1 NS NSHC
NSOURCFW
if (nsourcfw>0) then
NSOURCFW
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 vanishes, there is a single symmetric point source (with
strength unity as specified in the information file rpress.txt in this
directory rpress) located at the specified matching depth in the (interior)
fluid region.
When NSOURCFW is positive, on the other hand, there is a vertical source
array with NSOURCFW individual symmetric point sources. Point source number
n is located at the bottom of layer NSOFW0(n), referring to the layer
numbering prior to possible splitting as governed by the parameters
MHQ0,FREQMHQ0. The NSOFW0() must be positive integers and the medium
parameters must be continuous at the corresponding depths. The relative
strength (and phase) of point source number n is given by the complex
number AMSOURCFW(n), see the information file rpressfw.txt in this
directory rpress. However, the source information will be relevant only
when "IUTTYP=0".
In other respects, things are essentially as for a medium infile FILM to
either rpress or rpressfw, see the information files rpress.txt and
rpressfw.txt in this directory rpress. We only need to make a few remarks:
* The receiver depths are relevant only when "IUTTYP=0".
* There is an additional option 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 NSOURCFW must be positive
if "IUTTYP=0".
* The specified matching depth is actually only relevant for setting the
source depth when NSOURCFW vanishes and "IUTTYP=0". The "actual" matching
depth is chosen adaptively for each mode according to the
indicator-function method of section 4.1 in Ref. 7.
Format for the outfiles
-----------------------
Each outfile will consist of a number of data sets, each data set
corresponding to a particular mode. There is an initial line for each
data set, specifying the number of lines (depths) with corresponding data
that will follow. Each data line gives the depth and a complex number for
the mode quantity in question.
Before the 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 number of data sets (modes) is
written.
Subroutine directories
----------------------
The source code for the program and its subroutines can be found in the
following directories:
rpress/cpm/fw (program and main routines)
rpress/cpm/fw2 (other crucial routines)
rpress/cpm (other crucial routines)
rpress/cpm/single (routines that have different "multi" 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 rmodpfw was used for computing the mode functions in
Fig. 1 of Ref. 7:
ref7fig1.mx ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
2 ! icuttyp
1 ! isheettyp
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsubpre
100.0 0.0 1 ! freqmin freqint nfr
1.0 1.0 ! aladauer,preala
0 ! iurot
ref7fig1u ! mainmodu
2 ! iuttyp
ref7fig1p ! mainmodp
1 ! inolx
1 ! iystab
The medium infile ref7fig1.mx had the following appearance:
0.030 0.040 10 10 ! matching depth 0.040 km
0
0.000 0.000 0.0 0.0 0.0
0.000 0 0
3.500 2.000 0.9 0.0 0.0 0 0
0.010 0 0
3.500 2.000 0.9 0.0 0.0 10 0
0.010 0 0
4.000 2.200 0.95 0.3 0.2 0 0
0.025 0 0
4.000 2.200 0.95 0.3 0.2 15 0
0.025 0 0
3.500 1.600 0.98 0.1 0.1 0 0
0.030 0 0
3.500 1.600 0.98 0.1 0.1 5 0
0.030 0 0
1.500 0.000 1.0 0.0 0.0 0 0
0.040 0 0
1.400 0.000 1.0 0.1 0.0 -10 10
0.060 0 0
1.400 0.000 1.0 0.1 0.0 20 20
0.100 0 0
1.500 0.000 1.0 0.2 0.0 -40 40
0.100 0 0
1.400 0.500 1.2 0.3 0.1 0 0
0.110 0 0
1.400 0.500 1.2 0.3 0.1 10 0
0.110 0 0
3.500 2.000 2.0 0.2 0.3 0 0
0.150 0 0
3.500 2.000 2.0 0.2 0.3 40 0
0.150 0 0
5.000 2.67581763 3.2 0.0 0.0 0 0
0.400 0 0
5.000 2.67581763 3.2 0.0 0.0 250 0
0.400 0 0
5.000 2.67581763 8.2 0.0 0.0 0 0
0.600 0 0
5.000 2.67581763 8.2 0.0 0.0 200 0
0.600 0 0
8.000 4.000 2.0 0.1 0.1 0 0
0.750 0 0
8.000 4.000 2.0 0.1 0.1 150 0
0.750 0 0
8.000 4.000 2.0 0.1 0.1 0 0
The 'horizontal slowness' infile ref7fig1u.100 was computed using the
program rmodfaut. To enhance speed, it was here advantageous to avoid layer
splitting in the medium infile (used in ref7fig1.mx above to get output at
many depths). After some editing, for example choice of the desired modes
and removal of the three last data sets, the file ref7fig1u.100 got the
following appearance:
0 100.0 ! "lines to skip", freq
1 2 1 .0 0
4 ! nzero; now up; nzleft= 0 ! "isheet" rel "icuttyp=2"
0.5292851337591499 5.369426644550051E-04 1
0.8320830142430394 1.348341219028073E-03 1
2.39439017180364 4.917060644640056E-03 1
0.3775987732731447 2.951967798507614E-10 1
The indicator functions, as given in the right panels in Fig. 1 of Ref. 7,
are computed by subroutines in the file rpress/cpm/fw/fwsubytz.f and
rpress/cpm/fw2/dhivanfw2.f. These subroutines are always called by
rmodpfw, but the indicator functions are not normally written to disk.
Normally, only an "optimal matching depth" is needed (see section 4.1 of
Ref. 7). It as stored by the variable NSOOPT and the minimum value of the
indicator function is stored by the variable QIOPT.
We may note that in order to compute the power-density results of panels
c-d in Fig. 1 of Ref. 4, rmodpfw was used with IUTTYP=3.
Alternative methods
-------------------
Three methods for stable computation of normal modes as functions of depth
are actually proposed in Ref. 7:
1) The indicator-function method, on which the present program rmodpfw is
based.
2) The "inverse iteration" method, which is the alternative general method
of section 4.1 in Ref. 7. It is described in some more detail in section
5.3 of Paper V in Ref. 0.
3) The special method explained in connection with equation (49) in Ref. 7.
It is based on the proportionality theorem in Appendix C there.
All three methods were implemented and tested on a previous VAX/VMS
implementation of the rpress-package programs. The corresponding routines
can be found in the merged file rpress/cpm/fw/oldcode.
References
----------
See the information file read.me in this directory rpress.