Sven Ivansson Feb. 23, 1998 National Defence Research Establishment (FOA) S-17290 Stockholm, Sweden email: sveni@sto.foa.se General introduction to the program rpress (rpress = Reflectivity PRESSure) ------------------------------------------ The program rpress belongs to the rpress package. It can be used for cw ("continuous wave" with a monofrequency time dependence) computations as well as computation of transients. The underlying theory is described in the references. 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). The wave-field is created by a symmetric point source in the (interior) fluid region, and the receivers are also located there. Hence, backward propagation through solid regions will not be needed. (To get the full wave-field, the program rpressfw should be used.) Units are s for time, km for length and "10**12 kg" for mass. However, dB is "re 1 m" in the way explained below. The output concerns the pressure field in the fluid region given in "(10**15 N)/(km**2)". For the cw computations, an underlying harmonic time dependence s(t) according to "s(t)=exp(-i*omeg*t)" is understood, omeg>0 being the angular frequency and t being time. For computation of transients, the source wavelet s(t) can be arbitrarily specified. In any case, the point source has a diagonal moment tensor with all three diagonal elements equal to " s(t) " "(10**15 N)*(km)". The output for the cw computations concerns the complex pressure field, and phases are given in addition to the amplitudes. For typical applications, the source strength as specified would be unrealistically large, and one would have to scale down the results. However, the output in dB "re 1 m" would of course not be affected. 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 rpress --------------------------- LX LI(TOP) ICWORKDIM FILM MHQ0 FREQMHQ0 IAZIMI if ( iazimi not 0 ) then ALFA1AZIM end if ICUTTYP if (icuttyp=3) then AIMCUT end if EPSGLOBLN EPSGLOB EPSFAKT HFAKTMAX HOFF EPSMAGNIFMAX EPSFSUB HSTEPFSUBPRE PREEPS0 IXTYP0 NRUS0PRE IRUS0 JRUS0 if (lx<=1) then FREQMIN FREQINT NFR if ( freqmin<0.0 and freqint=0.0 ) then OFREQ(1:nfr) end if else TMIN DAUER URED PREALA IUROT FREQDOM if (freqdom<=0.0) then CR9FIL end if end if RMIN RINT if ( rmin<0.0 and rint=0.0 ) then OR(1:li) end if FILO if (lx>1) then MAINTR end if FILG ! filg=' ' allowed only when lx>1 if ( filg not ' ' ) then IDB MAINAM MAINPH end if INTREST IFILONPRE BESSFIL IADINTXFI if ( iadintxfi not 0 ) then SMAXPRE NELMAXPRE NBACKMAX IXTYP NRUS IRUS JRUS ARGDIFLIM1 ARGDIFLIM0 ARGDIFLIM2 end if NPARTS do jp = 1, nparts IBESSARR(jp) NPUARR(jp) if ( npuarr(jp) > 0 ) then UPHOERN(1:npuarr(jp),jp) else if ( npuarr(jp) = 0 ) then MAINMOD(jp) XYONENORM(jp) else if ( npuarr(jp) < 0 ) then AIMEND(jp) end if end do if (nparts>0) then HMAXDH HMAXRELWP if ( iadintxfi not 0 ) then EPSH01Q PREEPS HKVARPART LUNKVAR FAKTKVAR if ( hkvarpart>0.0 and faktkvar>0.0 ) then DDSUPMIN DDSUPMAX DDSUPFAKTMAX DERIVDRELMAX DRELMAX DRELGOALFAKT ARGDIFSTORE end if end if end if ICONTTYP if ( iconttyp not 0 ) then ITAIL(1) ITAIL(2) if ( itail(1) not 0 ) then IBESSTAIL(1) ARGRKTAIL(1) UP0TAIL(1) end if if ( itail(2) not 0 ) then IBESSTAIL(2) ARGRKTAIL(2) UP0TAIL(2) end if ANPHPMINTERM if ( iadintxfi not 0 ) then PREEPSTERM end if if (iconttyp=1) then IRCONT PREEPSCONT else if (iconttyp=2) then KMAXEU PREEPSEU end if end if Explanation of the parameters ----------------------------- LX is less than or equal to one for cw ("continuous wave" with a mono- frequency time dependence) computations whereas it is greater than one for computation of transients. In the latter case, LX specifies the number of time samples to be computed and it should be an integer power of two. Transients are computed by frequency synthesis. Hence, cw output can be obtained here as well. LI is the number of vertical receiver arrays, i.e. receiver ranges. ICWORKDIM is a parameter to give the dimension of certain work-space arrays. 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.) 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. The field results (for each frequency) can be computed with error control, implying that the computations can be refined adaptively by extrapolation until a specified accuracy is achieved. There are two levels of extrapolation: one "inner" for the 'horizontal slowness' integrals and one "outer" concerning refinement of the medium model, see section 5.1 in Ref. 0. For some details on the (inner) extrapolation process, see section II in Ref. 1. PREEPS0, IXTYP0, NRUS0PRE, IRUS0 and JRUS0 are control parameters for the outer extrapolation. (For many medium models, however, outer extrapolation will not be needed.) PREEPS0 is the specified relative accuracy. For a vertical receiver array at the range r km, it is transferred to an absolute accuracy requirement given by eps0=PREEPS0*pm0/r. (The reference pressure pm0 will be introduced below in connection with the parameter IDB. For r equal to zero, eps0 is set to PREEPS0*pm0.) 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. NRUS0PRE gives the depth of the extrapolation tables (typically about 8). A minus sign will imply that log output for the outer extrapolation is provided in outfiles ftn100, ftn101 and ftn102. 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). For cw computations, 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. For computation of transients, the nominal start time is TMIN s and the width of the time window is DAUER s. In particular, the time between successive samples will be DAUER/LX s. URED is the reduction slowness in s/km. Hence, for a receiver at the range r km from the source, the actual start time for the output trace will be TMIN+r*URED s. PREALA is a "temporary reduction factor during DAUER". For PREALA less than one, complex frequencies will be used during the computations in order to mitigate the effects of fft wrap-around. Recommended values are PREALA=0.01 or 0.001. 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. Rotation is recommended, it is indicated by IUROT=1 as opposed to IUROT=0. After rotation, a given real horizontal slowness will correspond to a real horizontal wavenumber. When FREQDOM is positive, it is related to the dominant frequency of the source in Hz. The spectrum for the moment-tensor wavelet s(t) of the source is 0.5*(1+cos(pi*0.5*f/freqdom)) for abs(f) < 2*freqdom and zero elsewhere. Hence, a kind of weighted sinc-function sum is obtained in the time domain. The time function is not strictly causal, but it is translated in time so that it appears to start at time zero. (It can be plotted using a program called sigrpress.) If FREQDOM is given as a nonpositive number, however, the source spectrum is taken from the file CR9FIL. (A program called cr9rpress can be used to compute the required spectrum of a desired wavelet given in the time domain.) Equivalent "freqmin freqint nfr" are written on the terminal. This can be useful for initial runs with the programs rmodfaut and rmodpfw if mode contributions are to be included. The LI vertical receiver arrays are located at ranges OR(1:li) km from the source. These ranges can be given explicitly, or they can be obtained as OR(ir)=RMIN+ir*(RINT/LI) for ir=1,..,li. In any case, they must form an increasing sequence of nonnegative values. FILO, character*32, is the name of an outfile which essentially contains a listing of the input parameters. For computation of transients, the traces will appear in outfiles with prefix MAINTR (character*32). There will be one outfile for each receiver depth as specified in the infile FILM, the extensions will simply be 1,2,3... Each such outfile will contain LI traces, one for each range. FILG, character*32, is the name of an outfile which contains the receiver geometry (depths and ranges in km). IDB specifies the type of output that is to be given for the cw point-source results. Consider our point source located in a homogeneous acoustic (no shear) whole-space with the same medium parameters as at the source depth in our original medium. Here our point source will require a certain power (which is in general different from the power required in our original medium). Next, consider a point source of the same frequency and power but in a homogeneous acoustic (no shear) whole-space with velocity 1.5 km/s, no anelastic attenuation and density "10**12 kg"/(km**3) (i.e. 1.0 kg/(dm**3)). Let "pm0" be the corresponding peak pressure, in (10**15 N)/(km**2), at a distance of 1 km. For IDB>=0, the field results will be given relative "pm01m=1000*pm0" and for IDB>0 in dB. For IDB =-1,-2 they will not be scaled in this way, IDB=-1 will imply output in single precision whereas IDB=-2 will imply double-precision output. For IDB=-3, the complex field results will be multiplied with the value of the Lame parameter lambda at the matching depth (as given in the infile FILM, see below). This is useful for tests of reciprocity. MAINAM is the prefix for the field-amplitude outfiles. 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 MAINAM is the string outampl and the third frequency to be considered is 25.3 Hz, the outfile for this frequency will be outampl.25X3.) In the same way, MAINPH is the outfile prefix for the phase values of the complex field. The phases will be given in degrees. A nonzero INTREST will imply an analytical computation of part of the field, as described at the end of section 6.3 in Paper IV of Ref. 0. This option should only be used for computations of the whole field without "mode parts" and "branch parts" (positive NPUARR() as described below), the integrands will be reduced by the correct amount. IFILONPRE indicates whether or not the Filon rule will be used for the primary trapezoidal integral estimates. The value 1 will imply Filon-rule estimates, whereas the value 0 will imply the ordinary trapezoidal rule. (Some other options are available as well, but they are not normally used.) BESSFIL gives the name of an infile (created using a program called besseltab) containing asymptotic approximations of a specified order to certain Bessel and Hankel functions. Normally, however, the character % is given, indicating an option to use standard first-order asymptotic approximations. IADINTXFI controls the type of numerical integration to be performed. A nonzero value will imply adaptive integration with error control, whereas the value zero will imply trapezoidal integration with a fixed step-size. For positive IADINTXFI, log output from the adaptive integration process will be provided on a ftn??? file with ??? replaced by the number IADINTXFI. (The numbers 6 and 66 will indicate the terminal, however.) Integration-step histograms as given in Fig. 5 (bottom frame) of Ref. 1 are written on ftn"???+1" provided IADINTXFI is 6 or greater than or equal to 80. Numerical integration takes place along a number of straight paths in complex 'horizontal slowness' space. For adaptive integration with error control (i.e. nonzero IADINTXFI), a number of control parameters are needed. The parameters we have in mind are SMAXPRE NELMAXPRE NBACKMAX IXTYP NRUS IRUS JRUS ARGDIFLIM1 ARGDIFLIM0 ARGDIFLIM2. A negative value of SMAXPRE will imply the corresponding positive value as a pre-set upper limit for the maximum allowed length of the "planned path segments" into which each straight path is divided. A nonnegative value of SMAXPRE will give a computed upper limit, this is the normal case. A positive value of NELMAXPRE gives the maximum allowed number of simultaneously active subintervals obtained by interval splitting for a planned path segment. Zero or a negative value will result in a computed maximum number, this is the normal case. NBACKMAX is the maximum number of back steps with wasted evaluations of the integrands for each planned path segment. (Normally, each evaluation of the integrands is stored in a stack. If the allocated space for the stack becomes occupied, the computations can nevertheless be continued provided a "back step" is taken and a new effort is made with a smaller planned path segment. The computations will be stopped when the number of back steps exceeds NBACKMAX.) IXTYP, NRUS, IRUS, and JRUS control the inner extrapolation process for computation of the integral along each planned path segment. They correspond in the obvious way to the parameters IXTYP0, NRUS0PRE, IRUS0 and JRUS0 for the outer extrapolation. (No minus sign for NRUS is allowed, however.) ARGDIFLIM1, ARGDIFLIM0 and ARGDIFLIM2 are certain limits for the maximal argument differences in degrees of the dispersion function between consecutive points used at a particular extrapolation depth for a particular subinterval used in the process of adaptive integration. One should have 0.0 <= argdiflim1 < argdiflim0 < argdiflim2. These parameters can be used to bring about an early subinterval splitting (which can be economical). The logic is adapted to cases of simple poles of the integrands, i.e. zeros of the dispersion function, close to the integration path. Early subinterval splitting will be enforced if the maximal argument differences in both halves of the subinterval exceed ARGDIFLIM0 and also if the maximal argument difference in one half exceeds ARGDIFLIM2 whereas no argument difference in the other half exceeds ARGDIFLIM1. In the computations, the field is built up from different contributions or parts. There will be "mode-path-branch parts" and "tail integrals". NPARTS is the number of path-mode-branch parts. For each such part, values must be assigned to some parameters. IBESSARR() specifies the type of Bessel/Hankel function to be used. The value 1 or 2 will prescribe J0 whereas -1 or -2 will prescribe 0.5*H0(1). The tables from BESSFIL will be used for the values +-1, whereas +-2 will imply very accurate computations using subroutines provided by Anders Bostrom. The value -3 for IBESSARR() is a special option implying quick computations of 0.5*H0(1) using first-order asymptotics and simple multiplications to move out in range. This option can only be used for vertical receiver arrays that are uniformly distributed in range. NPUARR() is positive to indicate a "path part", zero to indicate a "mode part", and negative to indicate a "branch part". When NPUARR() is positive, there is a path part and NPUARR() will give the number of corner points for the corresponding polygonal integration path. The complex horizontal slowness coordinates for the corner points, in their correct order, will be contained in UPHOERN(,). The P and S branch-points of a lower homogeneous half-space are specified by +-991 and +-992, respectively. Here, the plus sign indicates the first quadrant and the minus sign indicates the third quadrant. When NPUARR() is zero, there is a mode part and the string MAINMOD(), character*32, is the prefix for infiles containing modal horizontal slownesses and related information. The extensions are the frequencies in the same way as for the outfile prefix MAINAM. When XYONENORM() is nonzero, these infiles are obtained as the "mainmodu " outfiles from the program rmodfaut. When XYONENORM() vanishes, they are obtained as the "mainmodp " outfiles from the program rmodpfw run with the "iuttyp=0" option. The mode parts (residue contributions) can be computed in three different ways, depending on the value of XYONENORM(). In each case, a positive orientation of a surrounding integration path is assumed. -> When XYONENORM() is positive, the mode parts are computed by integration along a circular path around each modal horizontal slowness as specified in the MAINMOD() infile, see Appendix B of Paper IV in Ref. 0. The planned radius in s/km of the integration circle in the complex horizontal slowness plane is given by XYONENORM(). However, for each particular mode, the radius may be modified in order to ascertain that one and only one modal horizontal slowness is encircled. -> When XYONENORM() is zero, the appropriately excited mode parts as functions of depth are directly read from the MAINMOD() infile. Very little remains to be done, essentially only scaling with the appropriate range-dependent Hankel function. -> When XYONENORM() is negative, the mode parts are computed according to Theorem 1 of Ref. 7. This option should be used with caution, however, since the particular choice of matching depth can be critical. (In fact, different matching depths may be needed for different modes.) When NPUARR() is negative, finally, there is a branch part and branch-cut integrals will be computed. The value -1 indicates the P-cut of a lower homogeneous half-space, -2 indicates the S-cut and -12 indicates both cuts. (It is required that the P- and S-cuts are separated. Branch-cut integrals for an upper homogeneous half-space have not been implemented.) The integration will be performed along the cut(s), away from the real horizontal slowness axis, until the imaginary horizontal slowness part becomes AIMEND() s/km (AIMEND() is assumed to be positive). A change of integration variables is made, and this is essential close to a branch-point (see Appendix A of Paper IV in Ref. 0). HMAXDH and HMAXRELWP determine the maximum allowed (final) integration step-size in s/km along the integration path. This step-size is obtained as the minimum of HMAXDH and HMAXRELWP/(f*OR(LI)), where f is the current frequency in Hz. Hence, HMAXRELWP provides a dimension-less specification in terms of the minimum asymptotic Hankel-function period. (When OR(LI) is zero, HMAXRELWP will be irrelevant and HMAXDH will be taken.) To be cautious, cf. the Appendix to Ref. 1, HMAXRELWP should be taken somewhat smaller than 1.0. When adaptive integration with error control is performed, the available error tolerance is uniformly distributed along the integration path. However, even for a very small subinterval, over which integration is to be performed, the share of the available error tolerance will be at least EPSH01Q. The example in Remark 4 (ii) of section II in Ref. 1 corresponds to setting EPSH01Q to 0.01. PREEPS is a relative error tolerance for adaptive integration. It corresponds in the obvious way to the parameter PREEPS0 for the outer extrapolation process. If HKVARPART is positive, then the option for adaptively changing the integration path will be put into action; see the final paragraph of section III in Ref. 1. When the current integration step-size for the trapezoidal-sum estimates falls below HKVARPART times the "maximum allowed (final) integration step-size" (as determined from HMAXDH and HMAXRELWP), a path segment along the positive real axis will be moved into the fourth quadrant. If positive, HKVARPART might typically be 0.01. The absolute value of FAKTKVAR (typically 1.0) determines the relation between the depth of the roundabout path and the width of the original path segment. A minus sign for FAKTKVAR will indicate that the search for poles of the integrand within the enclosed region be omitted. The path segments to be moved will be temporarily stored on the file ftn???, where ??? is to be replaced by the number LUNKVAR. However, particular caution is needed when computing the field in a medium without anelastic attenuation but with some solid region (Ivansson and Karasalo, 1998, manuscript in preparation). A number of additional parameters are required if the check for poles of the integrand (appearing among the zeros of the dispersion function) is to be performed. This check is made using the argument principle for analytic functions (see Ref. 2) to determine the number of enclosed zeros. The parameters in question are DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX, DRELMAX,DRELGOALFAKT and ARGDIFSTORE. Since these parameters are taken from the program rmodfaut, we refer to the information file rmodfaut.txt (in this directory rpress) for the explanation. As mentioned, the computed field is built up from "mode-path-branch parts" and "tail integrals". Tail integrals will be present if ICONTTYP is nonzero. They can be computed in two different ways, depending on whether ICONTTYP is 1 or 2. There can be two tails, the first one (present if ITAIL(1) is nonzero) in the first quadrant and the second one (present if ITAIL(2) is nonzero) in the second quadrant; see equation (16) and Fig. 3 in Ref. 1. A nonzero ITAIL() can be equal to 1 or 2. For ITAIL() equal to 1, an asymptotic form of the Green's function will be used. When ITAIL() equals 2, the Green's function will be accurately computed. Each tail is specified by IBESSTAIL(), ARGRKTAIL() and UP0TAIL(). IBESSTAIL() specifies the type of Bessel/Hankel function to be used, in the same way as IBESSARR() for the mode-path-branch parts. ARGRKTAIL() gives the direction of the tail path, in degrees relative to the direction of the positive real horizontal slowness axis. Typical values are ARGRKTAIL(1)=+30.0 degrees and ARGRKTAIL(2)=-30.0 degrees. (Note: Compared to equation (16) and Fig. 3 in Ref. 1, there is a sign change and a change of direction for the second tail.) UP0TAIL(1) gives the start point of the first tail, whereas UP0TAIL(2) gives the end point of the second tail. Both these parameters are complex horizontal slowness values in s/km. The tail integrals are computed successively along asymptotic quarter-period segments for the Hankel functions. The maximum allowed (final) integration step-size in s/km along the integration path is the length of the current quarter-period segment divided by ANPHPMINTERM (typically 7.0). For adaptive integration, PREEPSTERM determines the accuracy requirements for each quarter-period computation in the same way as does PREEPS for the mode-path-branch parts. For ICONTTYP equal to 1, the chosen quarter-periods are those for the IRCONT:th vertical receiver array. (Typically, IRCONT is set to 1.) In this case, the tail computations will continue with quarter-period segments farther and farther away from UP0TAIL() until the contribution from the last quarter-period falls below a limit determined from PREEPSCONT in the same way as accuracy requirements are determined from PREEPSTERM. For ICONTTYP equal to 2, separate tail computations are performed for each vertical receiver array. (This will slow down the computations for large values of LI.) In this case, Euler's transformation is used for convergence acceleration (e.g. Dahlquist, Bjorck, Anderson 1974, Problem 3 in section 7.6). KMAXEU will denote the maximum number of columns in the "extrapolation table" whereas PREEPSEU will determine the error tolerance in the same way as accuracy requirements are determined from PREEPS. Format for the medium infile FILM --------------------------------- ZSOURC0 ZSOURC1 NS NSHC 0 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 ---------------------------------------------------- The infile FILM is read by the subroutine rmodel, and the parameters of the fluid-solid medium are set up along with the matching and receiver depths. Concerning the matching depth, see section 5 of Ref. 6. For reasons of simplicity, it is restricted in the program to be in the (interior) fluid region (see below). For rpress, the source depth will be at the matching depth. Hence, it will be in the (interior) fluid region, and the outgoing pressure field in the vicinity of the source appears as -exp(i*omeg*R/alpha) * (omeg**2)/(4*pi*(alpha**2)*R) where omeg is the angular frequency and R is the spherical distance to the source, while alpha is the medium velocity at the source. Concerning computation of transients, the factor omeg**2 will of course imply a sign-changed second derivative of the primarily given moment-tensor wavelet s(t). The initial "AA BB RRHO DBLA DBLB" data concern the upper boundary. A free or rigid boundary is specified by a vanishing or negative RRHO, respectively. An upper homogeneous half-space is specified by a positive RRHO, and the medium parameters are obtained from the data in the same way as described described below. ZZ is the depth to which the medium parameters in a data block refer. The ZZ value of a subsequent data block must be greater than (or equal to) that of the preceding data block. IIFSMTYP should be zero. (A nonzero value is possible to indicate a layer with continuous parameter variation, but we do not describe that option here.) IIMHQANI should be zero or one, the value one indicating frequency-dependent division into sublayers using the parameters MHQ0 and FREQMHQ0. (Preparation for transversely isotropic layers has been made, using other values of IIMHQANI as an indicator.) AA is the P-velocity and BB is the S-velocity. A vanishing value of BB indicates fluid. RRHO is the density. DBLA and DBLB are values for absorption (anelastic attenuation) of P- and S-waves, respectively, given in "dB per wavelength". NSH specifies the primary number of layers into which the depth interval starting from the ZZ value of the previous data block shall be divided. (As mentioned, further division may take place if IIMHQANI is one.) A vanishing NSH is used to indicate a discontinuity with a "new start". In the fluid case (i.e. when BB vanishes), NSH can be negative which indicates a nonhomogeneous region starting from the ZZ value of the previous data block. In this nonhomogeneous region, the square of the inverse P-velocity will be a linear function of depth. When abs(NSH) is greater than one, division into homogeneous (or nonhomogeneous with the linear inverse-square P-velocity) layers is done by interpolation for the medium parameters. NH is the number of receiver depths within the depth interval starting from the ZZ value of the preceding data block. If NH does not vanish, NSH must be divisible by NH. The receiver depths will be written in the outfile FILG. For rpress, receivers are only allowed in the (interior) fluid region. The last data block specifies the lower boundary. There is a lower homogeneous half-space if RRHO is positive. If RRHO of the last data block is negative or vanishes, on the other hand, the lower boundary will be rigid or free, respectively. The matching depth is specified by the parameters ZSOURC0, ZSOURC1, NS and NSHC. ZSOURC0 and ZSOURC1 must coincide with the ZZ values of two consecutive data blocks for the fluid region. NSHC must be equal to abs(NSH) of the latter data block. NS, finally, can take the values 0,1,..,NSHC and indicates precisely at which interface the matching depth is to be placed. Format for the outfiles ----------------------- Let NHYD be the number of receiver depths (i.e. the sum of the NH of the infile FILM). The output in a file with prefix MAINAM or MAINPH is given as a matrix of size NHYD x LI which is written row by row. Before the matrix, 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 integers NHYD and LI are written. Outfiles may also be created by the subroutine utspec. The idea is that the user could edit this routine to get his own desired output format. Concerning computation of transients, the outfiles with prefix MAINTR are "unformatted" in the Fortran sense (in order to reduce their size). Their structure can hopefully be understood from the Fortran code in the file rpressmain.f in the directory rpress/cpm. 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/bess (Bessel- and Hankel-function routines) rpress/slatec (slatec routines for complex Airy functions) rpress/ik (routines to allocate space for arrays) To conveniently produce the Green's function curves (essentially the Hankel-transform integrand) in Fig. 23 of Paper IV in Ref. 0, trapezoidal integration was used for simplicity (IADINTXFI set to zero) and special output was written to disk. Normally, trapezoidal integration is performed by the routine in the file rpress/int/trapezfi.f, but in this case the modified routine in the file rpress/int/trapezfiSPEC.fORIG was used. Example ------- The following input to rpress was used for computing the pressure curves in Fig. 1a of Ref. 4: (Note: The dB output from rpress is "re 1 m", whereas the results in the figure referred to were "re 1 km" and hence 60 dB higher.) 0 80 200000 ! lx,li(top),icworkdim ref4fig1a.m ! film 1 1.0 0 ! mhq0,freqmhq0, iazimi 2 ! icuttyp 1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps.. 1D-7 0.25 ! epsfsub,hstepfsub 0.00001 2 0 2 5 ! preeps0, ixtyp0,nrus0pre,irus0,jrus0 1.0 80.0 80 ! freqmin,freqint,nfr -1.0 0.0 ! rmin,rint 0.095 0.097 0.100 0.104 0.106 0.190 0.194 0.200 0.208 0.212 0.380 0.388 0.400 0.416 0.424 0.760 0.776 0.800 0.832 0.848 1.520 1.552 1.600 1.664 1.696 2.340 2.404 2.500 2.628 2.692 3.040 3.104 3.200 3.328 3.392 4.680 4.808 5.000 5.256 5.384 6.080 6.208 6.400 6.656 6.784 7.680 7.808 8.000 8.256 8.384 9.680 9.808 10.000 10.256 10.384 12.160 12.416 12.800 13.312 13.568 14.360 14.616 15.000 15.512 15.768 17.360 17.616 18.000 18.512 18.768 20.360 20.616 21.000 21.512 21.768 24.320 24.832 25.600 26.624 27.136 ref4fig1a.dat ! filo ref4fig1a.g ! filg 1 ! idb ref4fig1aa ! mainam ref4fig1ap ! mainph 0 ! intrest 1 ! ifilonpre % ! bessfil -1 ! iadintxfi 1.0 0 8 ! smaxpre,nelmaxpre,nbackmax 2 9 2 7 ! ixtyp,nrus,irus,jrus 0.0 10.0 430.0 ! argdiflim1,argdiflim0,argdiflim2 1 ! nparts 2 2 ! ibessarr(),npuarr() (0.0,0.0) (4.0,0.0) 0.01 0.95 ! hmaxdh,hmaxrelwp 0.01 0.000001 ! epsh01q,preeps 0.005 111 -1.0 ! hkvarpart,lunkvar,faktkvar 1 ! iconttyp 2 2 ! itail(1),itail(2) -2 30.0 ( 4.0,0.0) ! ibesstail(1) argrktail(1) up0tail(1) -2 -30.0 (-4.0,0.0) ! ibesstail(2) argrktail(2) up0tail(2) 7.0 ! anphpminterm 0.0000001 ! preepsterm 1 0.000001 ! ircont,preepscont The medium infile ref4fig1a.m had the following appearance: 0.000 0.050 1 1 ! matching depth = source depth = 0.050 km 0 0.000 0.000 0.000 0.0 0.0 0.000 0 0 1.450 0.000 1.000 0.0 0.0 0 0 0.050 0 0 1.450 0.000 1.000 0.0 0.0 1 1 0.100 0 0 1.450 0.000 1.000 0.0 0.0 1 0 0.100 0 0 1.550 0.120 1.600 0.4 0.6 0 0 0.103 0 0 1.550 0.120 1.600 0.4 0.6 1 0 0.103 0 0 4.000 2.300 2.500 0.2 0.3 0 0 Note the specification of the integration path with consistent "uphoern(,)" and "up0tail(1:2)"; cf. equation (16) and Fig. 3 in Ref. 1, although there is a sign change and a change of direction for the second tail there. The breakaway points +-4.0 s/km are not far enough out to include the Scholte wave, which may have a modal slowness of more than 1/0.120 s/km in this example. However, the Scholte wave should be negligible because the sediment layer is thin and the common source-receiver depth is not close to the water-sediment interface. References ---------- See the information file read.me in this directory rpress.