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 Version: 5.40
C Date: 2000, May 31
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' know about the new quantity for
C interpolation, the quantity should be registered in subroutine MTTQ.
C The user should supplement subroutine MTTQ with reading the filename
C of the output file with the interpolated quantity. At the same place,
C the user may also wish to supplement reading some other input
C parameters to be used in entries MTTQRI and MTTQRA for the calculation
C of the new quantity.  The supplemented input parameters must be
C declared in the SAVE statement in the 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     Reading possible input data
C     Registering a new quantity
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.
C Please, send your updated files to the author
C to include your improvements in the next version of 'mttq.for'.
C
C=======================================================================
C
      SUBROUTINE MTTQ
C
C Subroutine designed to read the input data describing the names
C of the output files with the user-defined interpolated quantities,
C and also to read any other input parameters, which will be used by the
C user. All the parameters to be read are assumed to be specified in
C the SEP format in the
C MTT input file.
C
C.......................................................................
C
C Dummy arguments of entries:
      CHARACTER*80 QNAME
      REAL QVALUE
C
C Subroutines and external functions required:
      EXTERNAL CIEROR,ERROR,RSEP3T,RSEP3R,RSEP3I
C     CIEROR ... File mtt.for.
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 /MTTC/ to store the information about triangles and rays.
      INCLUDE 'mtt.inc'
C mtt.inc
C
C-----------------------------------------------------------------------
C                                                    
C     Input data:
      REAL GBE11,GBE21,GBE31,GBE12,GBE22,GBE32
      REAL GBR11,GBR12,GBR22,GBY11,GBY22
      REAL AMPMAX
      SAVE GBE11,GBE21,GBE31,GBE12,GBE22,GBE32
      SAVE GBR11,GBR12,GBR22,GBY11,GBY22
      SAVE AMPMAX
C     Input data read in this subroutine 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 MTTQRI and used in entry MTTQRA:
      REAL R11,R12,R22,C11,C21,C12,C22,CDET
      SAVE R11,R12,R22,C11,C21,C12,C22,CDET
C     The values must be declared in the SAVE statement.
C
C     Temporary storage locations:
      CHARACTER*80 TXTERR
      INTEGER I
      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)
C
C.......................................................................
C     Reading the input data:
C
      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 RSEP3R('AMPMAX',AMPMAX,999999.)
C
C                                                    
C     User-defined quantity:
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.......................................................................
C
C     Registering quantities for bilinear interpolation within ray cells
C     and reading the filenames of the respective output files:
C
C     Components of the 4*4 paraxial ray propagator matrix in RCCS:
      CALL RSEP3T('MPQ11',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ11'
      ENDIF
      CALL RSEP3T('MPQ21',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ21'
      ENDIF
      CALL RSEP3T('MPQ31',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ31'
      ENDIF
      CALL RSEP3T('MPQ41',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ41'
      ENDIF
      CALL RSEP3T('MPQ12',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ12'
      ENDIF
      CALL RSEP3T('MPQ22',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ22'
      ENDIF
      CALL RSEP3T('MPQ32',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ32'
      ENDIF
      CALL RSEP3T('MPQ42',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ42'
      ENDIF
      CALL RSEP3T('MPQ13',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ13'
      ENDIF
      CALL RSEP3T('MPQ23',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ23'
      ENDIF
      CALL RSEP3T('MPQ33',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ33'
      ENDIF
      CALL RSEP3T('MPQ43',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ43'
      ENDIF
      CALL RSEP3T('MPQ14',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ14'
      ENDIF
      CALL RSEP3T('MPQ24',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ24'
      ENDIF
      CALL RSEP3T('MPQ34',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ34'
      ENDIF
      CALL RSEP3T('MPQ44',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='MPQ44'
      ENDIF
C
C     Widths of Gaussian beams:
      CALL RSEP3T('GBW11',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='GBW11'
        IF(GBY11.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 parameter GBW11, positive
C         real-valued parameter GBY11 must be specified in the input
C         SEP parameter or history file.
        END IF
      ENDIF
      CALL RSEP3T('GBW22',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='GBW22'
        IF(GBY22.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 parameter GBW22, positive
C         real-valued parameter GBY22 must be specified in the input
C         SEP parameter or history file.
        END IF
      ENDIF
      CALL RSEP3T('GBW1',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='GBW1'
        IF(GBY11.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 parameter GBW1, positive
C         real-valued parameter GBY11 must be specified in the input
C         SEP parameter or history file.
        END IF
      ENDIF
      CALL RSEP3T('GBW2',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='GBW2'
        IF(GBY22.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 parameter GBW2, positive
C         real-valued parameter GBY22 must be specified in the input
C         SEP parameter or history file.
        END IF
      ENDIF
      CALL RSEP3T('AMP',FOUT(NOUT+1),' ')
      IF (FOUT(NOUT+1).NE.' ') THEN
        NOUT=NOUT+1
        IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
        CHOUT(NOUT)='AMP'
      ENDIF
C                                                     
C     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 lines
C     designed to read the filename of the output file with the new
C     quantity and to register the quantity if the output filename is
C     not blank.
*     CALL RSEP3T('NAME',FOUT(NOUT+1),' ')
*     IF (FOUT(NOUT+1).NE.' ') THEN
*       NOUT=NOUT+1
*       IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR)
*       CHOUT(NOUT)='NAME'
*       insert lines to check the values of input data
*     ENDIF
C     Repeat the above lines for each new quantity.
C
C     End of subroutine MTTQ:
      RETURN
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
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)*GBE11+YLI(5)*GBE21+YLI(6)*GBE31)*UU
        U2=(YLI(4)*GBE12+YLI(5)*GBE22+YLI(6)*GBE32)*UU
        U3=(YLI(4)*  H13+YLI(5)*  H23+YLI(6)*  H33)*UU
C       Transformation matrix from the given vectors to RCCS
        C11=H11*GBE11+H21*GBE21+H31*GBE31
        C21=H12*GBE11+H22*GBE21+H32*GBE31
        C31=H13*GBE11+H23*GBE21+H33*GBE31
        C12=H11*GBE12+H21*GBE22+H31*GBE32
        C22=H12*GBE12+H22*GBE22+H32*GBE32
        C32=H13*GBE12+H23*GBE22+H33*GBE32
        CDET=C11*C22-C12*C21
C       Projection of the real part of matrix M to the given vectors
        R11=GBR11+U1*C31+C31*U1-C31*C31*U3
        R12=GBR12+U1*C32+C31*U2-C31*C32*U3
        R22=GBR22+U2*C32+C32*U2-C31*C31*U3
C       Sum of squares of Gaussian beam widths corresponding to GBY11
        IF(QNAME.EQ.'GBW11') THEN
          QVALUE=(C11*C11+C21*C21)/GBY11
          RETURN
        END IF
        IF(QNAME.EQ.'GBW1') THEN
          QVALUE=SQRT((C11*C11+C21*C21)/GBY11)
          RETURN
        END IF
C       Sum of squares of Gaussian beam widths corresponding to GBY22
        IF(QNAME.EQ.'GBW22') THEN
          QVALUE=(C12*C12+C22*C22)/GBY22
          RETURN
        END IF
        IF(QNAME.EQ.'GBW2') THEN
          QVALUE=SQRT((C12*C12+C22*C22)/GBY22)
          RETURN
        END IF
      END IF
C
C     Amplitudes:
      IF(QNAME.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
      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
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
        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
        IF(QNAME.EQ.'GBW11') THEN
          QVALUE=(Q111*Q111+Q211*Q211)/GBY11+(Q112*Q112+Q212*Q212)*GBY11
          RETURN
        END IF
        IF(QNAME.EQ.'GBW1') THEN
          QVALUE=SQRT(
     *          (Q111*Q111+Q211*Q211)/GBY11+(Q112*Q112+Q212*Q212)*GBY11)
          RETURN
        END IF
C       Sum of squares of Gaussian beam widths corresponding to GBY22
        IF(QNAME.EQ.'GBW22') THEN
          QVALUE=(Q121*Q121+Q221*Q221)/GBY22+(Q122*Q122+Q222*Q222)*GBY22
          RETURN
        END IF
        IF(QNAME.EQ.'GBW2') THEN
          QVALUE=SQRT(
     *          (Q121*Q121+Q221*Q221)/GBY22+(Q122*Q122+Q222*Q222)*GBY22)
          RETURN
        END IF
      END IF
C
C     Amplitudes:
      IF(QNAME.EQ.'AMP') THEN
        CALL AP21(GREEN)
        QVALUE=0.
        DO 10 I=15,32
          QVALUE=QVALUE+GREEN(I)**2
   10   CONTINUE
        QVALUE=SQRT(QVALUE)
        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