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.