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.60 C Date: 2001, October 2 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 SAVE GBE11,GBE21,GBE31,GBE12,GBE22,GBE32 SAVE GBR11,GBR12,GBR22,GBY11,GBY22 SAVE AMPMAX 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) 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' 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(22)='NAME' C If you are including more than one quantity, your second quantity C named NAME2 may be registered as * QNAMES(23)='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.EQ.'AMP') THEN C Amplitudes: CALL RSEP3R('AMPMAX',AMPMAX,999999.) 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.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.EQ.'AMP') THEN CALL AP21(1,PV,PV,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 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