Sven Ivansson Feb. 24, 1998 National Defence Research Establishment (FOA) S-17290 Stockholm, Sweden email: sveni@sto.foa.se General introduction to the program rmodfgr ------------------------------------------- The program rmodfgr belongs to the rpress package. It is intended for tracking the horizontal slowness of a particular mode as a function of frequency. To find the modal horizontal slowness at the initial frequency, one should first run the program rmodfaut. With the modal slowness obtained at one frequency, the derivative of the slowness with respect to frequency (which is closely related to the group-slowness) can be used as an initial prediction for the modal slowness at the next frequency. The required derivative is computed from derivatives of the dispersion function, see sections 5.2 and 3.1 of Ref. 7. An accurate slowness value can subsequently be computed using the secant rule (which is also used for final estimates in 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. Note, for example, that "solid-fluid-solid" media with surrounding half-spaces of either fluid or solid type can be handled. The program rmodfgr was a useful tool for verifying the theoretical results in Ref. 9 numerically. One should be aware that the mode tracking can fail. Automatic shifts of Riemann sheets as a branch-cut is approached have not been implemented. One should look manually for branch-cut crossings and restart the tracking as desired. (An undesired jump across a branch-cut cannot be excluded in the automatic procedure.) Furthermore, it is well known that dispersion curves of different modes may pass very close to each other and that they may even cross. Even if the frequency steps are chosen adaptively, there is a risk of jumping to the dispersion curve of another mode during the tracking process. Hence, it may be useful to check the 'modal slowness' pattern at isolated frequencies using the program rmodfaut. (No argument-variation checks are made within rmodfgr.) 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 rmodfgr ---------------------------- FILM MHQEFF ICUTTYP if (icuttyp=3) then AIMCUT end if ISHEETTYP EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX EPSFSUB HSTEPFSUBPRE FREQMIN FREQINT NFR if ( freqmin<0.0 and freqint=0.0 ) then OFREQ(1:nfr) end if ALADAUER PREALA IUROT0 UPSTART FILFU DSFMIN DSFMAX DSFFAKTMAX DSUMAX DRELGOALFAKT RSUG2 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 MHQEFF. Note that frequency-dependent layering and/or characteristic medium slownesses, as made possible in the other programs by the parameters "mhq0,freqmhq0,iazimi", is not allowed here. 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. PREALA is a "temporary reduction factor during the time ALADAUER". For PREALA less than one, complex frequencies will be used during the computations. IUROT0 is used to indicate whether or not given values for horizontal slowness are to be rotated in the complex plane because of nonreal frequencies. IUROT0 is irrelevant when PREALA=1.0, which is the normal case when running rmodfgr. Rotation is indicated by IUROT0=1 as opposed to IUROT0=0. After rotation, a given real horizontal slowness will correspond to a real horizontal wavenumber. One should note that internally, the computations are always done with "iurot=0" (no rotation). When IUROT0=1, rotation is performed to get the initial horizontal slowness "ustart" from UPSTART (see below), and rotation is also performed in connection with the output. UPSTART is the horizontal slowness (possibly rotated as specified by the parameter IUROT0) for the desired mode at the initial frequency as given by FREQMIN or OFREQ(1). FILFU is the name of the outfile in which three data sets will be provided: modal horizontal slownesses (possibly rotated as specified by the parameter IUROT0), derivatives of (unrotated) horizontal slowness with respect to frequency, and (unrotated) complex group-slownesses (i.e. derivatives of the modal wavenumbers with respect to angular frequency). For each data set, results for the different frequences will be provided in sequence. A good way to study the variation of mode shape with frequency is to run the program rmodpfw with "nfr=-1" and "mainmodu" as the FILFU outfile from rmodfgr. DSFMIN, DSFMAX and DSFFAKTMAX are control parameters for the frequency steps that are taken during the computations. The frequency step is adjusted adaptively to achieve a reasonable trade-off between efficiency and safety in the mode-tracking process. DSFMIN (typically 0.0 Hz) is the minimal frequency step that is handled. DSFMAX (typically 1.0 Hz) is the maximal frequency step. If a larger step is suggested, the next step-size will nevertheless be DSFMAX. DSFFAKTMAX (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 DSFFAKTMAX times the previous step-size. DSUMAX and DRELGOALFAKT are further control parameters for the frequency steps. The derivative of the modal slowness with respect to frequency is used as an initial prediction for the (unrotated) slowness at the next frequency. If the implied change of slowness is absolutely greater than DSUMAX (typically 0.01 s/km), a smaller frequency step is chosen. 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. RSUG2, REDUCFAKT, DDELT and XYTOL are tolerance parameters for the final determination of the modal horizontal slowness at the next frequency using the secant rule. RSUG2 (typically 0.01 s/km) is the half-width for a tolerance square that is put out centered at the prediction of (unrotated) horizontal slowness obtained from the derivative of the modal slowness with respect to frequency. The size of the square is reduced if it is crossed by a branch-cut, but it is not checked automatically if an undesired jump across a branch-cut has appeared. If the required zero (modal slowness) is not found within the square by the secant rule, one goes back and tries with a smaller frequency step. 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 of course chosen as the prediction from the derivative of the modal slowness with respect to frequency.) 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 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 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 rmodfgr. 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". The specified matching depth is not relevant for the rmodfgr computations, and the "actual" matching depth will be at "the initial ZZ". As mentioned in the information file rmodfaut.txt in this directory rpress, the number of layers should be kept as small as possible. This will enhance computational speed. Format for the outfiles ----------------------- The outfile FILFU will consist of three data sets. There is an initial line for each data set, specifying the number of lines with corresponding data that will follow. Each data line is for a particular frequency that is given together with a complex number for the desired quantity at that frequency. Before the three 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 (three) is written. 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/rut (general routines) rpress/slatec (slatec routines for complex Airy functions) Example ------- The following input to rmodfgr was used for tracking the Scholte mode in Fig. 1c of Ref. 7 from the original frequency 100 Hz down to 1 Hz: ref7fig1.mx ! film 1 ! mhqeff 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 -100.0 100 ! freqmin freqint nfr 1.0 1.0 ! aladauer,preala 0 ! iurot0 (2.39439017180364,4.917060644640056E-03) ! upstart slowfigfu ! filfu 0.0 1.0 5.0 ! dsfmin,dsfmax,dsffaktmax 0.01 0.75 ! dsumax,drelgoalfakt 0.01 ! rsug2 0.7 0.00001 0.000000001 ! reducfakt,ddelt,xytol The medium infile ref7fig1.mx has already been presented in the information file rmodpfw.txt (in this directory rpress), and we do not repeat it here. To enhance speed, however, it would be advantageous to remove the layer splitting when rmodfgr is run. To study the variation of mode shape with frequency, we may then run the program rmodpfw with the following input, note that "nfr=-1": 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 0.0 0.0 -1 ! freqmin freqint nfr 1.0 1.0 ! aladauer,preala 0 ! iurot slowfigfu ! mainmodu 1 ! iuttyp slowfigp ! mainmodp 1 ! inolx 1 ! iystab The results are actually very interesting, and they are included in the two figures that can be printed from the two postscript files slowfig.ps and slowfigp.ps, respectively, in this directory rpress. In both figures, the horizontal scale is for frequency and it covers -10 Hz to 110 Hz, with a tick mark every 10 Hz. The first figure, the one in slowfig.ps, shows the dispersion curve for our mode. The vertical slowness scale covers 0 to 4 s/km. Actually, three curves are shown. The curve for modal horizontal slowness is the one with a minimum of about 1 s/km at about 11 Hz. The corresponding group-slowness curve is also shown, with a minimum at about 10 Hz and a maximum at about 18 Hz. These two curves were obtained by running the program rmodfgr in the way shown above, tracking our mode from 100 Hz down to 1 Hz. For a solid-fluid-solid medium of this type, case (iv) of Theorem 6.1 in Ref. 9 predicts the existence of a slow wave such that the modal slowness tends to infinity as the inverse 2/3:th power of frequency when the frequency tends to zero. The third curve in our slowfig.ps figure shows the corresponding theoretical low-frequency slowness asymptote, with the constant factor obtained from equation (6.7) of Ref. 9. (As stated in Ref. 9, this equation is for the case with a rigid lower boundary for our fluid, but it applies to our actual case with a lower homogeneous solid half-space as well.) At low frequencies, there is a close agreement to the "rmodfgr" modal slowness curve, suggesting that our 100 Hz Scholte mode is transformed to this kind of slow mode when the frequency is lowered. To study the transformation process, we consider the second figure, the one in slowfigp.ps. It gives a "waterfall plot" for the modal vertical displacement, as obtained by running the program rmodpfw in the way shown above. The vertical scale is for depth, from 0.000 km at the surface down to 0.250 km at the bottom. The Scholte character of the mode persists down to 11 Hz. At lower frequencies, however, the main characteristic of the mode is an almost constant nonzero vertical displacement within the upper solid (the upper 0.030 km). This is consistent with the bending-wave appearance of the slow mode, as predicted by Theorem 6.6 of Ref. 9. Hence, we conclude that the Scholte wave for the lower fluid-solid interface of our solid-fluid-solid medium is transformed to a slow bending wave involving the upper solid as the frequency tends to zero. Both our figures show the real parts of slowness and vertical displacement, respectively. However, the corresponding imaginary parts are negligible. For a physically realistic visco-elastic medium model, absorptions in dB/wavelength should actually tend to zero as the frequency tends to zero. In Ref. 9, this was called the "truly" visco-elastic case. With vanishing anelastic attenuation in the medium, our imaginary parts would vanish identically. Using our programs, it turns out to be difficult to follow the mode in this example to lower frequencies than those of the order of a Hz. Numerical cancellations emerge among the terms of the propagator matrices, and a different numerical implementation would be needed. (The Runge-Kutta procedure implied by nonvanishing "iifsmtyp" in the medium infile "film" might be useful in this context.) References ---------- See the information file read.me in this directory rpress.