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.