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
C WFSRF.
C
C Version: 7.00
C Date: 2013, May 29
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     http://sw3d.cz/staff/bulant.htm
C     http://sw3d.cz/staff/klimes.htm
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 initial
C     point of a ray.
C     Calculating the value of a new quantity at the other
C     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 files with interpolated quantities
C     in program MTT, and to the list of the quantities which may
C     be written to the columns of the output file of 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*6  QNAME
      REAL QVALUE
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP3T,RSEP3R,RSEP3I
C     ERROR... File error.for.
C     RSEP3T,RSEP3R,RSEP3I... 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 TTDMAX,AMPMAX,AMP2D1,AMP2D2,AMP2D3,AMPKI1,AMPKI2,AMPKI3
      INTEGER KQIPV
      SAVE GBE11,GBE21,GBE31,GBE12,GBE22,GBE32
      SAVE GBR11,GBR12,GBR22,GBY11,GBY22
      SAVE TTDMAX,AMPMAX,AMP2D1,AMP2D2,AMP2D3,AMPKI1,AMPKI2,AMPKI3
      SAVE KQIPV
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,UU
      SAVE C11,C21,C12,C22,CDET,RU11,RU12,RU22,UU
C
C     Values set in entry MTTQRA and used in next invocation of MTTQRA:
      INTEGER ISRC0A,IWAV0A,IRAY0A,IPT0A,ISRC0D,IWAV0D,IRAY0D,IPT0D
      REAL GREEN(32),RN(6)
      SAVE    ISRC0A,IWAV0A,IRAY0A,IPT0A,ISRC0D,IWAV0D,IRAY0D,IPT0D,
     *     GREEN,RN
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,U,U1,U2,U3
      REAL Q111,Q211,Q121,Q221,Q112,Q212,Q122,Q222
      REAL HI(18),H(18),HUI(9),RM(6)
      REAL AUX
C
C.......................................................................
C
C     Memory management:
      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'
      QNAMES(23)='ampr11'
      QNAMES(24)='ampr21'
      QNAMES(25)='ampr31'
      QNAMES(26)='ampr12'
      QNAMES(27)='ampr22'
      QNAMES(28)='ampr32'
      QNAMES(29)='ampr13'
      QNAMES(30)='ampr23'
      QNAMES(31)='ampr33'
      QNAMES(32)='ampi11'
      QNAMES(33)='ampi21'
      QNAMES(34)='ampi31'
      QNAMES(35)='ampi12'
      QNAMES(36)='ampi22'
      QNAMES(37)='ampi32'
      QNAMES(38)='ampi13'
      QNAMES(39)='ampi23'
      QNAMES(40)='ampi33'
      QNAMES(41)='mti'
      QNAMES(42)='mx4'
      QNAMES(43)='mx5'
      QNAMES(44)='mx6'
      QNAMES(45)='mp4'
      QNAMES(46)='mp5'
      QNAMES(47)='mp6'
      QNAMES(48)='mtt11'
      QNAMES(49)='mtt12'
      QNAMES(50)='mtt22'
      QNAMES(51)='mtt13'
      QNAMES(52)='mtt23'
      QNAMES(53)='mtt33'
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(54)='name'
C     If you are including more than one quantity, your second quantity
C     named 'name2' may be registered as
*     QNAMES(55)='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.)
        ISRC0A=0
        IWAV0A=0
        IRAY0A=0
        IPT0A=0
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.)
C       Optional quasi-isotropic approximation of polarization vectors
        CALL RSEP3I('QIPV',KQIPV,0)
      ELSEIF (QNAME(1:3).EQ.'mtt') THEN
C       Travel-time derivatives:
        CALL RSEP3R('TTDMAX',TTDMAX,999999.)
        ISRC0D=0
        IWAV0D=0
        IRAY0D=0
        IPT0D=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
C use all the quantities from the common block and all his input
C parameters to compute the value VALUE of his quantity at the initial
C point of the ray. This value must be then assigned to the output
C quantity QVALUE.
C
C-----------------------------------------------------------------------
C
C     Imaginary travel time:
      IF(QNAME(1:3).EQ.'mti') THEN
        QVALUE=YI(2)
        RETURN
      ENDIF
C
C     Second travel-time derivatives:
      IF(QNAME(1:3).EQ.'mtt') THEN
        QVALUE=TTDMAX
        RETURN
      ENDIF
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     Source coordinates:
      IF    (QNAME(1:3).EQ.'mx4') THEN
        QVALUE=YI(3)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mx5') THEN
        QVALUE=YI(4)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mx6') THEN
        QVALUE=YI(5)
        RETURN
      ENDIF
C
C     Derivatives of travel time with respect to source coordinates:
      IF    (QNAME(1:3).EQ.'mp4') THEN
        QVALUE=-YI(6)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mp5') THEN
        QVALUE=-YI(7)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mp6') THEN
        QVALUE=-YI(8)
        RETURN
      ENDIF
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
C values at 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     Imaginary travel time:
      IF(QNAME(1:3).EQ.'mti') THEN
        QVALUE=Y(2)
        RETURN
      ENDIF
C
C     Second travel-time derivatives:
      IF(QNAME(1:3).EQ.'mtt') THEN
        IF ((ISRC.NE.ISRC0D).OR.(IWAVE.NE.IWAV0D).OR.(IRAY.NE.IRAY0D)
     *      .OR.(IPT.NE.IPT0D)) THEN
          ISRC0D=ISRC
          IWAV0D=IWAVE
          IRAY0D=IRAY
          IPT0D=IPT
          CALL AP03(0,HI,H,HUI)
          CALL AP08(0,H,HUI,RM,RN)
        ENDIF
        IF    (QNAME(1:5).EQ.'mtt11') THEN
          QVALUE=RN(1)
        ELSEIF(QNAME(1:5).EQ.'mtt12') THEN
          QVALUE=RN(2)
        ELSEIF(QNAME(1:5).EQ.'mtt22') THEN
          QVALUE=RN(3)
        ELSEIF(QNAME(1:5).EQ.'mtt13') THEN
          QVALUE=RN(4)
        ELSEIF(QNAME(1:5).EQ.'mtt23') THEN
          QVALUE=RN(5)
        ELSEIF(QNAME(1:5).EQ.'mtt33') THEN
          QVALUE=RN(6)
        ENDIF
        QVALUE=AMIN1(QVALUE,TTDMAX)
        QVALUE=AMAX1(QVALUE,-TTDMAX)
        RETURN
      ENDIF
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-C32*C32*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
        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
C       Sum of squares of Gaussian beam widths corresponding to GBY11
        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
        IF ((ISRC.NE.ISRC0A).OR.(IWAVE.NE.IWAV0A).OR.(IRAY.NE.IRAY0A)
     *      .OR.(IPT.NE.IPT0A)) THEN
          CALL AP21(KQIPV,GREEN)
          ISRC0A=ISRC
          IWAV0A=IWAVE
          IRAY0A=IRAY
          IPT0A=IPT
          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
            AUX=6.2831853*ABS(Y(20)*Y(25)-Y(21)*Y(24))
            AUX=AUX*(E11*E11+E21*E21)/ABS(Q111*Q221-Q121*Q211)
            AUX=SQRT(AUX)
            DO 7, I=15,32
              GREEN(I)=GREEN(I)*AUX
   7        CONTINUE
          ENDIF
        ENDIF
        IF     (QNAME(1:6).EQ.'ampr11') THEN
          QVALUE=GREEN(15)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi11') THEN
          QVALUE=GREEN(16)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr21') THEN
          QVALUE=GREEN(17)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi21') THEN
          QVALUE=GREEN(18)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr31') THEN
          QVALUE=GREEN(19)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi31') THEN
          QVALUE=GREEN(20)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr12') THEN
          QVALUE=GREEN(21)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi12') THEN
          QVALUE=GREEN(22)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr22') THEN
          QVALUE=GREEN(23)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi22') THEN
          QVALUE=GREEN(24)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr32') THEN
          QVALUE=GREEN(25)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi32') THEN
          QVALUE=GREEN(26)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr13') THEN
          QVALUE=GREEN(27)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi13') THEN
          QVALUE=GREEN(28)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr23') THEN
          QVALUE=GREEN(29)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi23') THEN
          QVALUE=GREEN(30)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampr33') THEN
          QVALUE=GREEN(31)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ELSEIF (QNAME(1:6).EQ.'ampi33') THEN
          QVALUE=GREEN(32)
          QVALUE=AMIN1(QVALUE,AMPMAX)
          QVALUE=AMAX1(QVALUE,-AMPMAX)
          RETURN
        ENDIF
C
        QVALUE=0.
        DO 10 I=15,32
          AUX=GREEN(I)
          AUX=AMIN1(AUX,AMPMAX)
          AUX=AMAX1(AUX,-AMPMAX)
          QVALUE=QVALUE+AUX**2
   10   CONTINUE
        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     Source coordinates:
      IF    (QNAME(1:3).EQ.'mx4') THEN
        QVALUE=YI(3)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mx5') THEN
        QVALUE=YI(4)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mx6') THEN
        QVALUE=YI(5)
        RETURN
      ENDIF
C
C     Derivatives of travel time with respect to source coordinates:
      IF    (QNAME(1:3).EQ.'mp4') THEN
        QVALUE=-YI(6)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mp5') THEN
        QVALUE=-YI(7)
        RETURN
      ELSEIF(QNAME(1:3).EQ.'mp6') THEN
        QVALUE=-YI(8)
        RETURN
      ENDIF
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
C should be 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             IGRID0=0 indicates no input grid.
C             Minus sign of IGRID0 indicates an undefined input grid
C             value in the input grid.  In this case, the index in array
C             RAM is -IGRID0.
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
      EXTERNAL UARRAY
      REAL UARRAY
C
C     Numerical parameter:
      REAL UNDEF
      SAVE UNDEF
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     Temporary variables used in entry RPGRD2:
      REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6,AUX7,AUX8
C
C     Initial value:
      DATA NGRID1/0/,NGRID2/0/
C
      UNDEF=UARRAY()
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))
        IGRID0=MAXRAM-NGRID+1
        IF(ICRT.LT.NCRT) THEN
          ISHIFT=NGRID*(NCRT-ICRT)
          DO 11 I=MAXRAM,IGRID0,-1
            RAM(I)=RAM(I-ISHIFT)
   11     CONTINUE
        END IF
        MAXRAM=MAXRAM-NGRID
        DO 12 I=IGRID0,MAXRAM
          IF(RAM(I).EQ.UNDEF) THEN
C           Indicating an undefined input grid value
            IGRID0=-IGRID0
            GO TO 13
          END IF
   12   CONTINUE
   13   CONTINUE
      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             IGRID0=0 indicates no input grid.
C             Minus sign of IGRID0 indicates an undefined input grid
C             value in the input grid.  In this case, the index in array
C             RAM is -IGRID0.
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 IF(IGRID0.GT.0) THEN
C       No undefined value in the grid
        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)
      ELSE
C       There is an undefined value in the grid
        AUX1=RAM(IGRID1-IGRID0)
        AUX2=RAM(IGRID2-IGRID0)
        AUX3=RAM(IGRID3-IGRID0)
        AUX4=RAM(IGRID4-IGRID0)
        AUX5=RAM(IGRID5-IGRID0)
        AUX6=RAM(IGRID6-IGRID0)
        AUX7=RAM(IGRID7-IGRID0)
        AUX8=RAM(IGRID8-IGRID0)
        IF(AUX1.NE.UNDEF.AND.
     *     AUX2.NE.UNDEF.AND.
     *     AUX3.NE.UNDEF.AND.
     *     AUX4.NE.UNDEF.AND.
     *     AUX5.NE.UNDEF.AND.
     *     AUX6.NE.UNDEF.AND.
     *     AUX7.NE.UNDEF.AND.
     *     AUX8.NE.UNDEF) THEN
C         No undefined value in the current grid cell
          VALUE=AGRID1*AUX1
     *         +AGRID2*AUX2
     *         +AGRID3*AUX3
     *         +AGRID4*AUX4
     *         +AGRID5*AUX5
     *         +AGRID6*AUX6
     *         +AGRID7*AUX7
     *         +AGRID8*AUX8
        ELSE
C         There is an undefined value in the current grid cell
          VALUE=0.
          IF(AGRID1.NE.0.) THEN
            IF(AUX1.NE.UNDEF) THEN
              VALUE=VALUE+AGRID1*AUX1
            ELSE
              GO TO 90
            END IF
          END IF
          IF(AGRID2.NE.0.) THEN
            IF(AUX2.NE.UNDEF) THEN
              VALUE=VALUE+AGRID2*AUX2
            ELSE
              GO TO 90
            END IF
          END IF
          IF(AGRID3.NE.0.) THEN
            IF(AUX3.NE.UNDEF) THEN
              VALUE=VALUE+AGRID3*AUX3
            ELSE
              GO TO 90
            END IF
          END IF
          IF(AGRID4.NE.0.) THEN
            IF(AUX4.NE.UNDEF) THEN
              VALUE=VALUE+AGRID4*AUX4
            ELSE
              GO TO 90
            END IF
          END IF
          IF(AGRID5.NE.0.) THEN
            IF(AUX5.NE.UNDEF) THEN
              VALUE=VALUE+AGRID5*AUX5
            ELSE
              GO TO 90
            END IF
          END IF
          IF(AGRID6.NE.0.) THEN
            IF(AUX6.NE.UNDEF) THEN
              VALUE=VALUE+AGRID6*AUX6
            ELSE
              GO TO 90
            END IF
          END IF
          IF(AGRID7.NE.0.) THEN
            IF(AUX7.NE.UNDEF) THEN
              VALUE=VALUE+AGRID7*AUX7
            ELSE
              GO TO 90
            END IF
          END IF
          IF(AGRID8.NE.0.) THEN
            IF(AUX8.NE.UNDEF) THEN
              VALUE=VALUE+AGRID8*AUX8
            ELSE
              GO TO 90
            END IF
          END IF
        END IF
      END IF
      RETURN
C
   90 CONTINUE
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, FGBR12,
C       FGBR22, FGBY11, FGBY22 or FTIME, are allowed only at the
C       gridpoints not used during interpolation within ray cells.
C       See the input data for program
C       mtt.for.
      RETURN
      END
C
C=======================================================================
C