C 
C Subroutine file 'mttq.for' to deal with the user-defined quantities
C to be interpolated by program MTT.
C The user is welcome to edit this file according to the instructions
C provided.  He will be able to read his input parameters, and to
C define quantities in the points of the rays.  The quantities, which
C are assumed to be real-valued, are then bilinearly interpolated by
C program MTT.
C
C Similarly, the quantities may be interpolated to the points on
C wavefronts using program WFSRF.
C
C Version: 5.90
C Date: 2005, January 12
C                                                  
C Coded by Petr Bulant and Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mails: bulant@seis.karlov.mff.cuni.cz,
C              klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C A list of the quantities which may already be interpolated by MTT can
C be found in the file mtt.for.
C
C If the user wishes to interpolate other quantities, he may add them by
C editting this file. Any user-defined quantity may be interpolated.
C The quantities stored in common block POINTC
C are at user's disposal for the calculation of the user-defined
C quantity at a point of ray. The user may define the value of his
C quantity at the initial point of a ray in subprogram entry MTTQRI.
C Similarly, the user may define the value of his quantity at the other
C points of a ray in subprogram entry MTTQRA.
C To let program 'mtt.for' (or 'wfsrf.for') know about the new quantity
C for interpolation, the quantity should be registered in subroutine
C MTTQ.
C At the subprogram entry MTTQD, the user may also wish to supplement
C reading some other input parameters to be used in entries MTTQRI and
C MTTQRA for the calculation of the new quantity.  The supplemented
C input parameters must be declared in the SAVE statement in the
C beginning of subroutine MTTQ.
C The filename (and possible other input parameters) should be read by
C individual invocations of the subroutines RSEP3T, RSEP3I and RSEP3R of
C file sep.for.
C
C Places to update when defining a new quantity for interpolation:
C     SAVE statement for possible input data
C     Registering a new quantity
C     Reading possible input data
C     Calculating the value of a new quantity at the
C             initial point of a ray
C     Calculating the value of a new quantity at the
C             other points of a ray
C     You may also wish to add your description of a new quantity to the
C             list of the names of output
C             files with interpolated quantities in program MTT, and
C             to the list of the quantities
C             which may be written to the columns of the output file of
C             program WFSRF.
C
C Please, send your updated files to the author
C to include your improvements in the next version of 'mttq.for'.
C
C=======================================================================
C
C     
C
      SUBROUTINE MTTQ
C
C Subroutine designed to register the quantities, which may be
C interpolated by programs MTT or WFSRF.
C.......................................................................
C
C Dummy arguments of entries:
      CHARACTER*5  QNAME
      REAL QVALUE
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP3T,RSEP3R,RSEP3I
C     ERROR ... File error.for.
C     RSEP3T,RSEP3R,RSEP3I ...
C             File sep.for.
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C pointc.inc.
C
C Common block /MTTWFC/:
      INCLUDE 'mttwf.inc'
C mttwf.inc
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C-----------------------------------------------------------------------
C                                                    
C     Logical unit number:
      INTEGER LU
      PARAMETER (LU=13)
C
C     Input data read in entry MTTQD:
      REAL GBE11,GBE21,GBE31,GBE12,GBE22,GBE32
      REAL GBR11,GBR12,GBR22,GBY11,GBY22
      REAL AMPMAX,AMP2D1,AMP2D2,AMP2D3,AMPKI1,AMPKI2,AMPKI3
      SAVE GBE11,GBE21,GBE31,GBE12,GBE22,GBE32
      SAVE GBR11,GBR12,GBR22,GBY11,GBY22
      SAVE AMPMAX,AMP2D1,AMP2D2,AMP2D3,AMPKI1,AMPKI2,AMPKI3
C     Input data read in entry MTTQD in order to be used in entries
C     MTTQRI and MTTQRA to calculate the quantities for interpolation
C     must be declared in the SAVE statement.
C
C     Values calculated in entry MTTQD and used in entries MTTQR*:
      INTEGER IGBE11,IGBE21,IGBE31,IGBE12,IGBE22,IGBE32
      INTEGER IGBR11,IGBR12,IGBR22,IGBY11,IGBY22,ITIME
      SAVE    IGBE11,IGBE21,IGBE31,IGBE12,IGBE22,IGBE32
      SAVE    IGBR11,IGBR12,IGBR22,IGBY11,IGBY22,ITIME
C
C     Values calculated in entry MTTQRI and used in entry MTTQRA:
      REAL C11,C21,C12,C22,CDET,RU11,RU12,RU22
      SAVE C11,C21,C12,C22,CDET,RU11,RU12,RU22
C
C     Temporary storage locations:
      INTEGER I
      REAL E11,E21,E31,E12,E22,E32,R11,R12,R22,Y11,Y22,TIME
      REAL C31,C32,H11,H21,H31,H12,H22,H32,H13,H23,H33,UU,U,U1,U2,U3
      REAL Q111,Q211,Q121,Q221,Q112,Q212,Q122,Q222
      REAL GREEN(32),PV(9),HI(18),H(18),HUI(9)
C
C.......................................................................
C
C     Memory management (see ram.inc):
      MAXRAM=MRAM
C
C     List of the character strings identifying the quantities which may
C     be interpolated:
      DO 5, I=1,MOUT
        QNAMES(I)=' '
  5   CONTINUE
      QNAMES( 1)='MPQ11'
      QNAMES( 2)='MPQ21'
      QNAMES( 3)='MPQ31'
      QNAMES( 4)='MPQ41'
      QNAMES( 5)='MPQ12'
      QNAMES( 6)='MPQ22'
      QNAMES( 7)='MPQ32'
      QNAMES( 8)='MPQ42'
      QNAMES( 9)='MPQ13'
      QNAMES(10)='MPQ23'
      QNAMES(11)='MPQ33'
      QNAMES(12)='MPQ43'
      QNAMES(13)='MPQ14'
      QNAMES(14)='MPQ24'
      QNAMES(15)='MPQ34'
      QNAMES(16)='MPQ44'
      QNAMES(17)='GBW11'
      QNAMES(18)='GBW1'
      QNAMES(19)='GBW22'
      QNAMES(20)='GBW2'
      QNAMES(21)='AMP'
      QNAMES(22)='AMPKI'
C                                                     
C     Registering a user-defined quantity:
C     Select string NAME identifying a new quantity and, simultaneously,
C     an input SEP parameter describing the name of the output file with
C     the new quantity. Then uncomment and modify the following line
C     designed to include the string into the list of quantities, which
C     may be interpolated.
*     QNAMES(23)='NAME'
C     If you are including more than one quantity, your second quantity
C     named NAME2 may be registered as
*     QNAMES(24)='NAME2'
C     and so on. The dimension MOUT of the array QNAMES is defined in
C     include file mttwf.inc.
C     End of subroutine MTTQ:
      RETURN
C
C=======================================================================
C
C     
C
      ENTRY MTTQD(QNAME)
C
C Entry designed to read input data necessary for computation of
C the interpolated quantities. This entry is called for each quantity
C to be interpolated.
C
      IF     (QNAME(1:3).EQ.'AMP') THEN
C       Amplitudes:
        CALL RSEP3R('AMPMAX',AMPMAX,999999.)
C       Amplitudes of 2-D Green function instead of 3-D Green function
        CALL RSEP3R('AMP2D1',AMP2D1,0.)
        CALL RSEP3R('AMP2D2',AMP2D2,0.)
        CALL RSEP3R('AMP2D3',AMP2D3,0.)
C       Modification of amplitudes for Kirchhoff integral
        CALL RSEP3R('AMPKI1',AMPKI1,0.)
        CALL RSEP3R('AMPKI2',AMPKI2,0.)
        CALL RSEP3R('AMPKI3',AMPKI3,0.)
      ELSEIF (QNAME(1:3).EQ.'GBW') THEN
C       Widths of Gaussian beams:
        CALL RSEP3R('GBE11',GBE11,1.)
        CALL RSEP3R('GBE21',GBE21,0.)
        CALL RSEP3R('GBE31',GBE31,0.)
        CALL RSEP3R('GBE12',GBE12,0.)
        CALL RSEP3R('GBE22',GBE22,1.)
        CALL RSEP3R('GBE32',GBE32,0.)
        CALL RSEP3R('GBR11',GBR11,0.)
        CALL RSEP3R('GBR12',GBR12,0.)
        CALL RSEP3R('GBR22',GBR22,0.)
        CALL RSEP3R('GBY11',GBY11,0.)
        CALL RSEP3R('GBY22',GBY22,0.)
        CALL RPGRD0(LU,1,'FTIME' ,ITIME)
        CALL RPGRD0(LU,2,'FGBE11',IGBE11)
        CALL RPGRD0(LU,2,'FGBE21',IGBE21)
        CALL RPGRD0(LU,2,'FGBE31',IGBE31)
        CALL RPGRD0(LU,2,'FGBE12',IGBE12)
        CALL RPGRD0(LU,2,'FGBE22',IGBE22)
        CALL RPGRD0(LU,2,'FGBE32',IGBE32)
        CALL RPGRD0(LU,2,'FGBR11',IGBR11)
        CALL RPGRD0(LU,2,'FGBR12',IGBR12)
        CALL RPGRD0(LU,2,'FGBR22',IGBR22)
        CALL RPGRD0(LU,2,'FGBY11',IGBY11)
        CALL RPGRD0(LU,2,'FGBY22',IGBY22)
      ENDIF
C
C                                                    
C     Reading input data for user-defined quantity:
C     The statement for reading the input data may be entered either
C     right here, or may be put to the above IF/ELSEIF/ENDIF
C     construction. If put to the above IF/ELSEIF/ENDIF construction,
C     the data will be read only if interpolating the corresponding
C     quantity. If put right here, the data will be read in any case.
C
C     Integer value IVALUE of integer input SEP parameter NAMEI is read
C     by statement
*     CALL RSEP3I('NAMEI',IVALUE,0)
C     where 0 is the default value of the parameter.
C
C     Real value VALUE of the integer or real-valued input SEP parameter
C     NAMER is read by statement
*     CALL RSEP3R('NAMER',VALUE,0.)
C     where 0. is the default value of the parameter.
C
C     End of subroutine MTTQD:
      RETURN
C
C=======================================================================
C
C     
C
      ENTRY MTTQRI(QNAME,QVALUE)
C
C Entry designed to compute the user-defined quantity at an initial
C point of a ray.
C
C Input:
C     QNAME...Character string, which characterizes the quantity.
C
C Output:
C     QVALUE..Value of the quantity at the initial point of a ray.
C
C All the quantities computed by CRT at the initial point of the ray
C under consideration are now stored in the common block
C pointc.inc. The user may use all the
C quantities from the common block and all his input parameters to
C compute the value VALUE of his quantity at the initial point of the
C ray. This value must be then assigned to the output quantity QVALUE.
C
C-----------------------------------------------------------------------
C
C     Components of the 4*4 paraxial ray propagator matrix in RCCS:
      IF(QNAME(1:3).EQ.'MPQ') THEN
        IF(QNAME.EQ.'MPQ11') THEN
          QVALUE=1.
        ELSE IF(QNAME.EQ.'MPQ21') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ31') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ41') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ12') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ22') THEN
          QVALUE=1.
        ELSE IF(QNAME.EQ.'MPQ32') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ42') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ13') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ23') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ33') THEN
          QVALUE=1.
        ELSE IF(QNAME.EQ.'MPQ43') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ14') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ24') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ34') THEN
          QVALUE=0.
        ELSE IF(QNAME.EQ.'MPQ44') THEN
          QVALUE=1.
        END IF
        RETURN
      END IF
C
C     Widths of Gaussian beams:
      IF (QNAME(1:3).EQ.'GBW') THEN
        CALL RPGRD1(1,YI(1))
        CALL RPGRD2(ITIME,TIME,YI(1))
        CALL RPGRD1(2,TIME)
        CALL RPGRD2(IGBE11,E11,GBE11)
        CALL RPGRD2(IGBE21,E21,GBE21)
        CALL RPGRD2(IGBE31,E31,GBE31)
        CALL RPGRD2(IGBE12,E12,GBE12)
        CALL RPGRD2(IGBE22,E22,GBE22)
        CALL RPGRD2(IGBE32,E32,GBE32)
        CALL RPGRD2(IGBR11,R11,GBR11)
        CALL RPGRD2(IGBR12,R12,GBR12)
        CALL RPGRD2(IGBR22,R22,GBR22)
        CALL RPGRD2(IGBY11,Y11,GBY11)
        CALL RPGRD2(IGBY22,Y22,GBY22)
        IF ((QNAME(1:4).EQ.'GBW1').AND.(Y11.LE.0.)) THEN
C         MTTQ-01
          CALL ERROR('MTTQ-01: Input parameter GBY11 is not positive')
C         To calculte and interpolate the sum of squares of Gaussian
C         beam widths identified by input parameters GBW11 or GBW1,
C         positive real-valued parameter GBY11 must be specified in the
C         input SEP parameter file or stored in file given by input
C         parameter FGBY11.
        END IF
        IF ((QNAME(1:4).EQ.'GBW2').AND.(Y22.LE.0.)) THEN
C         MTTQ-02
          CALL ERROR('MTTQ-02: Input parameter GBY22 is not positive')
C         To calculte and interpolate the sum of squares of Gaussian
C         beam widths identified by input parameters GBW22 or GBW2,
C         positive real-valued parameter GBY22 must be specified in the
C         input SEP parameter file or stored in file given by input
C         parameter FGBY22.
        END IF
C       Ray-centred coordinate system at the initial point
        UU=YI(6)**2+YI(7)**2+YI(8)**2
        U=SQRT(UU)
        H13=YI( 6)/U
        H23=YI( 7)/U
        H33=YI( 8)/U
        H11=YI( 9)
        H21=YI(10)
        H31=YI(11)
        H12=H23*H31-H33*H21
        H22=H33*H11-H13*H31
        H32=H13*H21-H23*H11
C       Minus slowness gradient in the given vectors and RCCS
        U1=(YLI(4)*E11+YLI(5)*E21+YLI(6)*E31)*UU
        U2=(YLI(4)*E12+YLI(5)*E22+YLI(6)*E32)*UU
        U3=(YLI(4)*H13+YLI(5)*H23+YLI(6)*H33)*UU
C       Transformation matrix from the given vectors to RCCS
        C11=H11*E11+H21*E21+H31*E31
        C21=H12*E11+H22*E21+H32*E31
        C31=H13*E11+H23*E21+H33*E31
        C12=H11*E12+H21*E22+H31*E32
        C22=H12*E12+H22*E22+H32*E32
        C32=H13*E12+H23*E22+H33*E32
        CDET=C11*C22-C12*C21
C       Projection of the real part of matrix M to the given vectors
        RU11=U1*C31+C31*U1-C31*C31*U3
        RU12=U1*C32+C31*U2-C31*C32*U3
        RU22=U2*C32+C32*U2-C31*C31*U3
C       Sum of squares of Gaussian beam widths corresponding to Y11
        IF(QNAME.EQ.'GBW11') THEN
          QVALUE=(C11*C11+C21*C21)/Y11
          RETURN
        END IF
        IF(QNAME.EQ.'GBW1') THEN
          QVALUE=SQRT((C11*C11+C21*C21)/Y11)
          RETURN
        END IF
C       Sum of squares of Gaussian beam widths corresponding to Y22
        IF(QNAME.EQ.'GBW22') THEN
          QVALUE=(C12*C12+C22*C22)/Y22
          RETURN
        END IF
        IF(QNAME.EQ.'GBW2') THEN
          QVALUE=SQRT((C12*C12+C22*C22)/Y22)
          RETURN
        END IF
      END IF
C
C     Amplitudes:
      IF(QNAME(1:3).EQ.'AMP') THEN
        QVALUE=AMPMAX
        RETURN
      END IF
C                                                    
C     User-defined quantity:
C     Uncomment and modify the following lines, where NAME is the name
C     of your quantity.
*     IF(QNAME.EQ.'NAME') THEN
*       insert lines to calculate the value of a new quantity
*       QVALUE=value of a new quantity
*       RETURN
*     ENDIF
C     Repeat the above lines for each new quantity.
C
C     End of entry MTTQRI:
      RETURN
C
C=======================================================================
C
C     
C
      ENTRY MTTQRA(QNAME,QVALUE)
C
C Entry designed to compute the user-defined quantity at other than
C initial point of a ray.
C
C Input:
C     QNAME...Character string, which characterizes the quantity.
C
C Output:
C     QVALUE..Value of the quantity at the point of a ray.
C
C All the quantities computed by CRT at the point of the ray
C under consideration are now stored in the common block
C pointc.inc together with the values at
C the initial point of the ray. The user may use all the
C quantities from the common block and all his input parameters to
C compute the value VALUE of his quantity at the point under
C consideration. This value must be then assigned to the output
C quantity QVALUE.
C
C-----------------------------------------------------------------------
C
C     Components of the 4*4 paraxial ray propagator matrix in RCCS:
      IF(QNAME(1:3).EQ.'MPQ') THEN
        IF(QNAME.EQ.'MPQ11') THEN
          QVALUE=Y(12)
        ELSEIF (QNAME.EQ.'MPQ21') THEN
          QVALUE=Y(13)
        ELSEIF (QNAME.EQ.'MPQ31') THEN
          QVALUE=Y(14)
        ELSEIF (QNAME.EQ.'MPQ41') THEN
          QVALUE=Y(15)
        ELSEIF (QNAME.EQ.'MPQ12') THEN
          QVALUE=Y(16)
        ELSEIF (QNAME.EQ.'MPQ22') THEN
          QVALUE=Y(17)
        ELSEIF (QNAME.EQ.'MPQ32') THEN
          QVALUE=Y(18)
        ELSEIF (QNAME.EQ.'MPQ42') THEN
          QVALUE=Y(19)
        ELSEIF (QNAME.EQ.'MPQ13') THEN
          QVALUE=Y(20)
        ELSEIF (QNAME.EQ.'MPQ23') THEN
          QVALUE=Y(21)
        ELSEIF (QNAME.EQ.'MPQ33') THEN
          QVALUE=Y(22)
        ELSEIF (QNAME.EQ.'MPQ43') THEN
          QVALUE=Y(23)
        ELSEIF (QNAME.EQ.'MPQ14') THEN
          QVALUE=Y(24)
        ELSEIF (QNAME.EQ.'MPQ24') THEN
          QVALUE=Y(25)
        ELSEIF (QNAME.EQ.'MPQ34') THEN
          QVALUE=Y(26)
        ELSEIF (QNAME.EQ.'MPQ44') THEN
          QVALUE=Y(27)
        END IF
        RETURN
      END IF
C
C     Widths of Gaussian beams:
      IF (QNAME(1:3).EQ.'GBW') THEN
        CALL RPGRD1(1,Y(1))
        CALL RPGRD2(ITIME,TIME,Y(1))
        CALL RPGRD1(2,TIME)
        IF (IGBE11.NE.0.OR.IGBE21.NE.0.OR.IGBE31.NE.0.OR.
     *      IGBE12.NE.0.OR.IGBE22.NE.0.OR.IGBE32.NE.0.) THEN
C         Recalculating matrices C and RU:
          CALL RPGRD2(IGBE11,E11,GBE11)
          CALL RPGRD2(IGBE21,E21,GBE21)
          CALL RPGRD2(IGBE31,E31,GBE31)
          CALL RPGRD2(IGBE12,E12,GBE12)
          CALL RPGRD2(IGBE22,E22,GBE22)
          CALL RPGRD2(IGBE32,E32,GBE32)
C         Minus slowness gradient in the given vectors and RCCS
          U1=(YLI(4)*E11+YLI(5)*E21+YLI(6)*E31)*UU
          U2=(YLI(4)*E12+YLI(5)*E22+YLI(6)*E32)*UU
          U3=(YLI(4)*H13+YLI(5)*H23+YLI(6)*H33)*UU
C         Transformation matrix from the given vectors to RCCS
          C11=H11*E11+H21*E21+H31*E31
          C21=H12*E11+H22*E21+H32*E31
          C31=H13*E11+H23*E21+H33*E31
          C12=H11*E12+H21*E22+H31*E32
          C22=H12*E12+H22*E22+H32*E32
          C32=H13*E12+H23*E22+H33*E32
          CDET=C11*C22-C12*C21
C         Projection of the real part of matrix M to the given vectors
          RU11=U1*C31+C31*U1-C31*C31*U3
          RU12=U1*C32+C31*U2-C31*C32*U3
          RU22=U2*C32+C32*U2-C31*C31*U3
        END IF
C       Transforming matrices of geometrical spreading to given vectors
        Q111= Y(12)*C11+Y(16)*C21
        Q211= Y(13)*C11+Y(17)*C21
        Q121= Y(12)*C12+Y(16)*C22
        Q221= Y(13)*C12+Y(17)*C22
        Q112=( Y(20)*C22-Y(24)*C12)/CDET
        Q212=( Y(21)*C22-Y(25)*C12)/CDET
        Q122=(-Y(20)*C21+Y(24)*C11)/CDET
        Q222=(-Y(21)*C21+Y(25)*C11)/CDET
C       Dependence of the widths on real part R11,R12,R22 of matrix M
        CALL RPGRD2(IGBR11,R11,GBR11)
        CALL RPGRD2(IGBR12,R12,GBR12)
        CALL RPGRD2(IGBR22,R22,GBR22)
        R11=RU11+R11
        R12=RU12+R12
        R22=RU22+R22
        Q111=Q111+Q112*R11+Q122*R12
        Q211=Q211+Q212*R11+Q222*R12
        Q121=Q121+Q112*R12+Q122*R22
        Q221=Q221+Q212*R12+Q222*R22
C       Sum of squares of Gaussian beam widths corresponding to GBY11
        CALL RPGRD2(IGBY11,Y11,GBY11)
        CALL RPGRD2(IGBY22,Y22,GBY22)
        IF ((QNAME(1:4).EQ.'GBW1').AND.(Y11.LE.0.)) THEN
C         MTTQ-03
          CALL ERROR('MTTQ-03: Input parameter GBY11 is not positive')
C         To calculte and interpolate the sum of squares of Gaussian
C         beam widths identified by input parameters GBW11 or GBW1,
C         positive real-valued parameter GBY11 must be specified in the
C         input SEP parameter file or stored in file given by input
C         parameter FGBY11.
        END IF
        IF ((QNAME(1:4).EQ.'GBW2').AND.(Y22.LE.0.)) THEN
C         MTTQ-04
          CALL ERROR('MTTQ-04: Input parameter GBY22 is not positive')
C         To calculte and interpolate the sum of squares of Gaussian
C         beam widths identified by input parameters GBW22 or GBW2,
C         positive real-valued parameter GBY22 must be specified in the
C         input SEP parameter file or stored in file given by input
C         parameter FGBY22.
        END IF
        IF(QNAME.EQ.'GBW11') THEN
          QVALUE=(Q111*Q111+Q211*Q211)/Y11+(Q112*Q112+Q212*Q212)*Y11
          RETURN
        END IF
        IF(QNAME.EQ.'GBW1') THEN
          QVALUE=SQRT(
     *          (Q111*Q111+Q211*Q211)/Y11+(Q112*Q112+Q212*Q212)*Y11)
          RETURN
        END IF
C       Sum of squares of Gaussian beam widths corresponding to GBY22
        IF(QNAME.EQ.'GBW22') THEN
          QVALUE=(Q121*Q121+Q221*Q221)/Y22+(Q122*Q122+Q222*Q222)*Y22
          RETURN
        END IF
        IF(QNAME.EQ.'GBW2') THEN
          QVALUE=SQRT(
     *          (Q121*Q121+Q221*Q221)/Y22+(Q122*Q122+Q222*Q222)*Y22)
          RETURN
        END IF
      END IF
C
C     Amplitudes:
      IF(QNAME(1:3).EQ.'AMP') THEN
        CALL AP21(1,PV,PV,GREEN)
        QVALUE=0.
        DO 10 I=15,32
          QVALUE=QVALUE+GREEN(I)**2
   10   CONTINUE
        IF(AMP2D1.NE.0..OR.AMP2D2.NE.0..OR.AMP2D3.NE.0.) THEN
C         Amplitudes of 2-D Green function
          CALL AP03(0,HI,H,HUI)
          E11=HI(1)*AMP2D1+HI(2)*AMP2D2+HI(3)*AMP2D3
          E21=HI(4)*AMP2D1+HI(5)*AMP2D2+HI(6)*AMP2D3
          Q111=Y(12)*E11+Y(16)*E21
          Q211=Y(13)*E11+Y(17)*E21
          Q121=Y(20)*E21-Y(24)*E11
          Q221=Y(21)*E21-Y(25)*E11
          QVALUE=QVALUE*6.2831853*ABS(Y(20)*Y(25)-Y(21)*Y(24))
          QVALUE=QVALUE*(E11*E11+E21*E21)/ABS(Q111*Q221-Q121*Q211)
        END IF
        QVALUE=SQRT(QVALUE)
        IF(QNAME.EQ.'AMPKI') THEN
C         Modification of amplitudes for Kirchhoff integral:
          QVALUE=QVALUE*(YI(6)*AMPKI1+YI(7)*AMPKI2+YI(8)*AMPKI3)
          QVALUE=QVALUE*2.*YLI(3)/ABS(YI(6)**2+YI(7)**2+YI(8)**2)
        END IF
        QVALUE=AMIN1(QVALUE,AMPMAX)
        RETURN
      END IF
C                                                     
C     User-defined quantity:
C     Uncomment and modify the following lines, where NAME is the name
C     of your quantity.
*     IF(QNAME.EQ.'NAME') THEN
*       insert lines to calculate the value of a new quantity
*       QVALUE=value of a new quantity
*       RETURN
*     ENDIF
C     Repeat the above lines for each new quantity.
C
C     End of entry MTTQRA:
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RPGRD0(LU,KTIME,NAME,IGRID0)
C     Dummy arguments of all entries:
      CHARACTER*(*) NAME
      INTEGER LU,KTIME,IGRID0
      REAL TIME,VALUE,DEF
C
C Reading Ray-Parameter GRiD values into array RAM.
C Variable MAXRAM in file ram.inc should be
C defined before invocation of this subroutine.
C
C Input:
C     LU...   Logical unit number to be used for the input.
C     KTIME...Flag denoting time axis:
C             KTIME=1: Travel time along ray.
C             KTIME=2: Two-way travel time.
C     NAME... Name of the SEP parameter to specify the filename of the
C             ray-parameter grid of values.
C Output:
C     IGRID0..Index in array RAM corresponding to the first gridpoint
C             of the stored part of the ray-parameter grid.
C             Index IGRID0 is used as input of entry RPGRD2 identifying
C             the quantity to interpolate.
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C pointc.inc.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C-----------------------------------------------------------------------
C
C     Numerical parameter:
      REAL UNDEF
      PARAMETER (UNDEF=-999999.)
C
C     Values defined in subroutine RPGRD0 and used in entry RPGRD1:
      INTEGER NTT,NTIME,NPAR1,NPAR2,NWAVES,NGRID1,NGRID2
      SAVE    NTT,NTIME,NPAR1,NPAR2,NWAVES,NGRID1,NGRID2
      REAL    OTT,DTT,OTIME,DTIME,OPAR1,DPAR1,OPAR2,DPAR2
      SAVE    OTT,DTT,OTIME,DTIME,OPAR1,DPAR1,OPAR2,DPAR2
C
C     Values defined in entry RPGRD1 and used in entry RPGRD2:
      INTEGER NGRID
      SAVE    NGRID
      INTEGER IGRID1,IGRID2,IGRID3,IGRID4,IGRID5,IGRID6,IGRID7,IGRID8
      SAVE    IGRID1,IGRID2,IGRID3,IGRID4,IGRID5,IGRID6,IGRID7,IGRID8
      REAL    AGRID1,AGRID2,AGRID3,AGRID4,AGRID5,AGRID6,AGRID7,AGRID8
      SAVE    AGRID1,AGRID2,AGRID3,AGRID4,AGRID5,AGRID6,AGRID7,AGRID8
C
C     Temporary variables used in entry RPGRD0:
      CHARACTER*80 FILE
      INTEGER NCRT,ICRT,ISHIFT,I
C
C     Temporary variables used in entry RPGRD1:
      INTEGER ITT,IPAR1,IPAR2,JTT,JPAR1,JPAR2,N
      REAL ATT,APAR1,APAR2
C
C     Initial value:
      DATA NGRID1/0/,NGRID2/0/
C
C-----------------------------------------------------------------------
C
      CALL RSEP3T(NAME,FILE,' ')
      IF(FILE.EQ.' ') THEN
        IGRID0=0
      ELSE
        IF(NGRID1.EQ.0.AND.NGRID2.EQ.0) THEN
          CALL RSEP3I('NTT'  ,NTT  ,1)
          CALL RSEP3R('OTT'  ,OTT  ,0.)
          CALL RSEP3R('DTT'  ,DTT  ,1.)
          CALL RSEP3I('NTIME',NTIME,1)
          CALL RSEP3R('OTIME',OTIME,0.)
          CALL RSEP3R('DTIME',DTIME,1.)
          CALL RSEP3I('NPAR1',NPAR1,1)
          CALL RSEP3R('OPAR1',OPAR1,0.)
          CALL RSEP3R('DPAR1',DPAR1,1.)
          CALL RSEP3I('NPAR2',NPAR2,1)
          CALL RSEP3R('OPAR2',OPAR2,0.)
          CALL RSEP3R('DPAR2',DPAR2,1.)
          CALL RSEP3I('NWAVES',NWAVES,1)
          CALL RSEP3I('NCRT' ,NCRT ,1)
          CALL RSEP3I('ICRT' ,ICRT ,1)
          IF(ICRT.LT.1.OR.NCRT.LT.ICRT) THEN
C           MTTQ-05
            CALL ERROR('MTTQ-05: Incorrect input parameter ICRT')
          END IF
        END IF
        IF(KTIME.LE.1) THEN
          NGRID1=NTT*NPAR1*NPAR2*NWAVES
          NGRID=NGRID1
        ELSE
          NGRID2=NTIME*NPAR1*NPAR2*NWAVES
          NGRID=NGRID2
        END IF
        IF(MAXRAM-NGRID*NCRT.LT.0) THEN
C         MTTQ-06
          CALL ERROR('MTTQ-06: Too small array RAM(MRAM)')
C         Too small array RAM(MRAM) to allocate input ray-parameter
C         grid values.  If possible, increase dimension MRAM in include
C         file ram.inc.
        END IF
        CALL RARRAY(LU,FILE,'FORMATTED',.TRUE.,UNDEF,NGRID*NCRT,
     *                                    RAM(MAXRAM-NGRID*NCRT+1))
        DO 11 I=MAXRAM-NGRID*NCRT+1,MAXRAM
          IF(RAM(I).EQ.UNDEF) THEN
C           MTTQ-07
            CALL ERROR('MTTQ-07: Undefined input grid value')
C           Undefined input grid values in files, given by parameters
C           FGBE11, FGBE21, FGBE31, FGBE12, FGBE22, FGBE32, FGBR11,
C           FGBR12, FGBR22, FGBY11, FGBY22 or FTIME, are not allowed
C           in this version.  See the input data for program
C           mtt.for.
          END IF
   11   CONTINUE
        IGRID0=MAXRAM-NGRID+1
        IF(ICRT.LT.NCRT) THEN
          ISHIFT=NGRID*(NCRT-ICRT)
          DO 12 I=MAXRAM,IGRID0,-1
            RAM(I)=RAM(I-ISHIFT)
   12     CONTINUE
        END IF
        MAXRAM=MAXRAM-NGRID
      END IF
      RETURN
C
C=======================================================================
C
C     
C
      ENTRY RPGRD1(KTIME,TIME)
C     INTEGER KTIME
C     REAL TIME
C
C Entry to be called once for each point of a ray, before entry RPGRD2.
C This entry calculates the indices and weights for trilinear
C interpolation within given ray-parameter grid.
C
C Input:
C     KTIME...Flag denoting time axis:
C             KTIME=1: Travel time along ray.
C             KTIME=2: Two-way travel time.
C     TIME... Travel time at the point of the ray.
C             TIME=YI(1) at the initial point of the ray,
C             TIME=Y(1) at other points of the ray.
C
C No output.
C
C-----------------------------------------------------------------------
C
      IF(KTIME.LE.1) THEN
        NGRID=NGRID1
      ELSE
        NGRID=NGRID2
      END IF
C
      IF(NGRID.GT.0) THEN
        IF(KTIME.LE.1) THEN
          ATT=(TIME-OTT)/DTT
          N=NTT
        ELSE
          ATT=(TIME-OTIME)/DTIME
          N=NTIME
        END IF
        ITT=NINT(ATT-0.5)
        IF(ITT.LT.0) THEN
          ITT=0
          JTT=0
          ATT=0.
        ELSE IF(ITT.GE.N-1) THEN
          ITT=N-1
          JTT=0
          ATT=0.
        ELSE
          JTT=1
          ATT=ATT-FLOAT(ITT)
        END IF
C
        APAR1=(YI(20)-OPAR1)/DPAR1
        IPAR1=NINT(APAR1-0.5)
        IF(IPAR1.LT.0) THEN
          IPAR1=0
          JPAR1=0
          APAR1=0.
        ELSE IF(IPAR1.GE.NPAR1-1) THEN
          IPAR1=NPAR1-1
          JPAR1=0
          APAR1=0.
        ELSE
          JPAR1=N
          APAR1=APAR1-FLOAT(IPAR1)
        END IF
C
        APAR2=(YI(21)-OPAR2)/DPAR2
        IPAR2=NINT(APAR2-0.5)
        IF(IPAR2.LT.0) THEN
          IPAR2=0
          JPAR2=0
          APAR2=0.
        ELSE IF(IPAR2.GE.NPAR2-1) THEN
          IPAR2=NPAR2-1
          JPAR2=0
          APAR2=0.
        ELSE
          JPAR2=N*NPAR1
          APAR2=APAR2-FLOAT(IPAR2)
        END IF
C
        IF(IWAVE.GT.NWAVES) THEN
C         MTTQ-08
          CALL ERROR('MTTQ-08: More elementary waves than NWAVES')
C         Index IWAVE of the elementary wave is greater than number
C         NWAVES of elementary waves corresponding to ray-parameter
C         grid.
        END IF
        IGRID1=ITT+N*(IPAR1+NPAR2*(IPAR2+NPAR2*(IWAVE-1)))
        IGRID2=IGRID1+JTT
        IGRID3=IGRID1+JPAR1
        IGRID4=IGRID2+JPAR1
        IGRID5=IGRID1+JPAR2
        IGRID6=IGRID2+JPAR2
        IGRID7=IGRID3+JPAR2
        IGRID8=IGRID4+JPAR2
        AGRID1=(1.-ATT)*(1.-APAR1)*(1.-APAR2)
        AGRID2=    ATT *(1.-APAR1)*(1.-APAR2)
        AGRID3=(1.-ATT)*    APAR1 *(1.-APAR2)
        AGRID4=    ATT *    APAR1 *(1.-APAR2)
        AGRID5=(1.-ATT)*(1.-APAR1)*    APAR2
        AGRID6=    ATT *(1.-APAR1)*    APAR2
        AGRID7=(1.-ATT)*    APAR1 *    APAR2
        AGRID8=    ATT *    APAR1 *    APAR2
      END IF
      RETURN
C
C=======================================================================
C
C     
C
      ENTRY RPGRD2(IGRID0,VALUE,DEF)
C     INTEGER IGRID
C     REAL VALUE
C
C Entry to calculate the given quantity at the current point of the ray
C by bilinear interpolation within the ray-parameter grid.
C
C Input:
C     IGRID0..Index in array RAM identifying the gridded quantity.
C     DEF...  Default value in case of no grid specified.
C Output:
C     VALUE...Interpolated value of the quantity under consideration.
C
C-----------------------------------------------------------------------
C
      IF(NGRID.EQ.0.OR.IGRID0.EQ.0) THEN
        VALUE=DEF
      ELSE
        VALUE=AGRID1*RAM(IGRID0+IGRID1)
     *       +AGRID2*RAM(IGRID0+IGRID2)
     *       +AGRID3*RAM(IGRID0+IGRID3)
     *       +AGRID4*RAM(IGRID0+IGRID4)
     *       +AGRID5*RAM(IGRID0+IGRID5)
     *       +AGRID6*RAM(IGRID0+IGRID6)
     *       +AGRID7*RAM(IGRID0+IGRID7)
     *       +AGRID8*RAM(IGRID0+IGRID8)
      END IF
      RETURN
      END
C
C=======================================================================
C