C
C Program INV1TT to evaluate the derivatives of the travel time with C respect to the model coefficients. C C Version: 5.10 C Date: 1997, September 30 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: klimes@seis.karlov.mff.cuni.cz C C....................................................................... C C Just a preliminary demo version, illustrating the usage of the C routines designed to calculate the variations of travel times with C respect to model parameters (FORTRAN77 source code files 'var.for', C 'spsp.for', 'ap.for'). C C Program INV1TT assumes all model parameters (coefficients) stored in C the common block /VALC/ as in the submitted versions of user-defined C model specification FORTRAN77 source code files 'srfc.for', 'parm.for' C and 'val.for'. Thus, unlike the other parts of the complete ray C tracing, the INV1TT program cannot work with user's modifications of C the subroutines SRFC1, SRFC2, PARM1, and PARM2. C C....................................................................... C C C Description of data files: C C Main input data file read from the interactive device (*): C (1) 'INV1TT' C 'INV1TT' name of the input file described below. C Default: 'INV1TT'='inv1tt.dat'. C C C Input data INV1TT: C This main data file contains the names of other input and output C files. The data are read in by the list directed input (free C format). In the list of input data below, each numbered paragraph C indicates the beginning of a new input operation (new READ C statement). All input variables are of the type CHARACTER. Only C the first 80 characters of the strings are significant. C (1) 'MODEL','POINTS','FTT','DATA','INV1LOG',/ C 'MODEL'... Input data file containing the model parameters. C See the file 'model.for' of the package 'MODEL'. C Description of MODEL C 'POINTS'... Input data file containing the coordinates of shot and C receiver points. Its structure is described below. C Ignored (not opened) if the DATA filename below is C blank. C Description of file POINTS C 'FTT'... Input data file containing the travel times from the C field measurements. Its structure is described below. C Ignored (not opened) if the DATA filename below is C blank. C Description of file FTT C 'DATA'..Output file containing the computed travel times and their C derivatives with respect to the model coefficients. C Its structure is described below. C If the filename is blank, no file with objective prior C information is generated. C Description of file DATA C 'INV1LOG'... Output log file. Just for your information. C Not opened and generated if the DATA filename above is C blank. C Description of file IN1LOG C /... An obligatory slash at the end of line. C Default: 'MODEL'='model.dat', 'DATA'='data.out', C 'INV1LOG'='inv1log.out'. C (2) DIST,VPOWER,/ C DIST... Maximum distance between the source or receiver point and C the initial or end point of a synthetic ray. C VPOWER... Velocity power for the computation of the travel-time C check sum. If the VPOWER-th velocity power is expressed, C in all blocks of the model, in terms of explicit functions C of model coordinates, linearly homogeneous in their C coefficients (e.g. B-splines), the travel time minus its C check sum (see the log output file) should be zero within C rounding errors. Otherwise, the check sum may have no C sense. C VPOWER=0: no check sum is evaluated and printed. C /... An obligatory slash at the end of line. C Default: VPOWER=0.0. C (3) Any times the following data: C 'CRT-R','CRT-S','CRT-I' C 'CRT-R'... File with the quantities stored along rays (see C C.R.T.5.5.1). C Description of file CRT-R C 'CRT-S'... File with the quantities at the points of C intersection of rays with the specified surface at which C the receivers are situated for the case of two-point ray C tracing (see C.R.T.5.5.2). C If this filename is not blank, just the two-point rays C with minimum travel time at each receiver are considered. C If this filename is blank, all rays are taken into C account. C Attention: All rays taken into account must start in some C of the specified sources and terminate in some of the C specified receivers, see the input file FTT. C Description of file CRT-S C 'CRT-I'... File with the quantities at the initial points C of rays, corresponding to the above file rays or points of C intersection (see C.R.T.6.1). C Description of file CRT-I C (4) / (a slash). C Example of data INV1TT C C C Input data file POINTS: C This data file contains the coordinates of shot and receiver C points. The data are read in by the list directed input C (free format). In the list of input data below, each numbered C paragraph indicates the beginning of a new input operation (new C READ statement). The CHARACTER strings are explicitly mentioned C in this description. Otherwise, if the first letter of the C symbolic name of the input variable is I-N, the corresponding C value in input data must be of the type INTEGER. Otherwise, the C input parameter is of the type REAL. C (1) Several strings terminated by / (a slash). C (2) List of the sources and receivers: Any times the following data: C (2.1) POINT,X1,X2,X3 C POINT...CHARACTER*11 string containing the name of the source or C receiver point. C X1,X2,X3... Coordinates of the source or receiver point. C (3) / (a slash) or the end of file. C Example of data POINTS C C C Input data file FTT: C This data file contains the travel times from the field C measurements. The data are read in by the list directed input C (free format). In the list of input data below, each numbered C paragraph indicates the beginning of a new input operation (new C READ statement). The CHARACTER strings are explicitly mentioned C in this description. Otherwise, if the first letter of the C symbolic name of the input variable is I-N, the corresponding C value in input data must be of the type INTEGER. Otherwise, the C input parameter is of the type REAL. C (1) Several strings terminated by / (a slash). C (2) List of the travel times: Any times the following data (2.1): C (2.1) SOURCE,RECEIVER,TFIELD,TERR C SOURCE..CHARACTER*11 string containing the name of the source C point. C RECEIVER... CHARACTER*11 string containing the name of the C receiver point. The source and receiver points may be C mutually interchanged. C TFIELD..Travel time from a field measurement. C TERR... Error of the travel time from a field measurement. C (3) / (a slash) or the end of file. C Example of data FTT C C C Output file DATA: C (1) ND,NM,(INDM(I),I=1,NM) C ND... The number of data (i.e. the number of equations). C NM... Number of unknown model coefficients. C INDM... Indices of the model coefficients influencing the travel C times. The indices correspond to the relative location in C the memory. C B-spline coefficients are listed in the same order as the C grid velocities in the file 'MODEL'. C (2) ND-times the following data (2.1): C (2.1) KD,RD,ED,WD,(GD(I),I=1,NM) C KD... Index of the field travel time within the file 'ftt'. C The field travel times are indexed consecutively 1,2,3,... C RD... Field travel time minus the computed synthetic travel C time. In the case of multiple two-point rays, the first C arrival of them is considered. C ED... Prior error of the above travel time difference. C Identical to TERR, file 'FTT', part (3). I.e. the C synthetic travel time is assumed to be sufficiently C accurate. C WD... Field travel time (stored for the purpose to assess a C posterior data misfit covariance weighting matrix). C GD... Derivatives of the synthetic travel time with respect to C the model coefficients INDM. C C C Output log file INV1LOG: C (1) For each considered ray: C (1.1) SOURCE,RECEIVER,TFIELD,TDIF,SDIST,RDIST,CHECK C SOURCE..Name of the source point. C RECEIVER... Name of the receiver point. C TFIELD..Travel time from a field measurement. C TDIF... Field travel time minus the minimum synthetic travel time. C SDIST...Distance between the source and the initial point of the C synthetic ray. C RDIST...Distance between the receiver and the end point of the C synthetic ray. C CHECK...Synthetic travel time minus the travel time resulting from C the derivatives of the theoretical travel time with C respect to the model coefficients. This quantity should C not exceed in order the numerical error of the synthetic C travel time. C In this version defined just for the models described in C the terms of velocity. C C----------------------------------------------------------------------- C C Common block /VALC/: INCLUDE 'val.inc' C val.inc C None of the storage locations of the common block are altered. C C Common block /POINTC/: INCLUDE 'pointc.inc' C pointc.inc C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of working arrays: C INTEGER MP,MT PARAMETER (MP=1000) PARAMETER (MT=4000) CHARACTER*11 POINT(MP) INTEGER KS(MT),KR(MT) REAL COOR(3,MP),TFIELD(MT),TERR(MT),TTMIN(MT) EQUIVALENCE (KS ,RAM ) EQUIVALENCE (KR ,RAM( MT+1)) EQUIVALENCE (TFIELD,RAM(2*MT+1)) EQUIVALENCE (TERR ,RAM(3*MT+1)) EQUIVALENCE (TTMIN ,RAM(4*MT+1)) EQUIVALENCE (COOR ,RAM(5*MT+1)) C INTEGER INDM(NPAR) REAL SUM(NPAR) C INDM... Indices of the unknown model parameters. EQUIVALENCE (SUM , RAM(5*MT+3*MP+1)) C C Output data: subjective prior information: REAL CS(1) C C----------------------------------------------------------------------- C C Filenames: CHARACTER*80 FILE2,FILE0,FILE1,FILE3,FILE4,FILE6 C C Logical unit numbers: INTEGER LU0,LU1,LU2,LU3,LU4,LU5,LU7,IU1,IU2 PARAMETER (LU0=10) PARAMETER (LU1=11) PARAMETER (LU2=12) PARAMETER (LU3=13) PARAMETER (LU4=14) PARAMETER (LU5=15) PARAMETER (LU7=17) C C Input data: CHARACTER*1 TEXT(10) CHARACTER*11 SRC,REC INTEGER NP,NT REAL DIST,DIST2,VPOWER C POINT...Names of the source and receiver points. C NP... Number of source and receiver points. C NT... Number of field travel times. C KS(I)...Index of the source point corresponding to the I-th field C travel time. C KR(I)...Index of the receiver point corresponding to the I-th C field travel time. C DIST... Maximum distance between the source or receiver point and C the initial or end point of a synthetic ray. C DIST2...DIST**2 C VPOWER... Velocity power for the computation of the travel-time C check SUM. C COOR... Coordinates of the source or receiver points. C TFIELD..Field travel times. C TERR... Field travel time errors. C C Output data: variations of the synthetic travel time: INTEGER NSUM,NM C NM... Number of the unknown model parameters. C C Auxiliary storage locations: INTEGER IS,IR,IT,ND,IRAYTT,I,J,K REAL TT,TTAUX,TDIF,XI1,XI2,XI3,XE1,XE2,XE3,PI(6),PE(6) REAL SI,SE,RI,RE,SDIST,RDIST C IS... Index of the source point. C IR... Index of the receiver point. C IT... Index of the field travel time. C ND... Number of synthetic travel times corresponding to the C field travel times. C IRAYTT..Index of the last processed ray. C I,J,K...Temporary storage locations. C TTMIN...Minimum synthetic travel times corresponding to the C individual field travel times. C TT... Synthetic travel time. C TTAUX...Temporary storage location. C TDIF... Field travel time minus the minimum synthetic travel time. C XI1,XI2,XI3,XE1,XE2,XE3... Coordinates of the initial and end C points of the last processed ray. C PI,PE...Slowness vectors at the initial and end points of the C last processed ray. C SI,SE,RI,RE... Squares of the distances between the source or C receiver points and the initial or end points of the ray. C SDIST...Distance between the source and the initial point of the C synthetic ray. C RDIST...Distance between the receiver and the end point of the C synthetic ray. C C....................................................................... C C Opening data files and reading the input data: C C Main input data file read from the interactive device (*): WRITE(*,'(A)') ' Enter the name of the main input data file: ' FILE0='inv1.dat' READ(*,*) FILE0 WRITE(*,'(A)') '+ ' C C Input data INV1TT: OPEN(LU0,FILE=FILE0,STATUS='OLD') FILE2='model.dat' FILE4='data.out' FILE6='inv1log.out' READ(LU0,*) FILE2,FILE0,FILE1,FILE4,FILE6 VPOWER=0. READ(LU0,*) DIST,VPOWER DIST2=DIST*DIST OPEN(LU4,FILE=FILE2,STATUS='OLD') CALL MODEL1(LU4) CLOSE(LU4) C C Number of unknown model parameters: CALL SOFT(2,0,0,0,0,0,0,0.,NM,INDM,CS) WRITE(*,'(A,I4,A)') '+',NM,' model parameters' C C....................................................................... C C Reading source and receiver points: C OPEN(LU4,FILE=FILE0,STATUS='OLD') NP=-1 READ(LU4,*,END=2) TEXT 1 CONTINUE NP=NP+1 IF(NP.GE.MP) THEN C INV1TT-01 PAUSE 'Error INV1TT-01: Too many source and receiver points' STOP END IF POINT(NP+1)=' ' READ(LU4,*,END=2) POINT(NP+1),(COOR(I,NP+1),I=1,3) IF(POINT(NP+1).NE.' ') GO TO 1 2 CONTINUE CLOSE(LU4) C C Reading field travel times: C OPEN(LU4,FILE=FILE1,STATUS='OLD') NT=0 READ(LU4,*,END=2) TEXT 3 CONTINUE NT=NT+1 IF(NT.GT.MT) THEN C INV1TT-02 PAUSE 'Error INV1TT-02: Too many field travel times' STOP END IF SRC=' ' READ(LU4,*,END=9) SRC,REC,TFIELD(NT),TERR(NT) IF(SRC.EQ.' ') THEN GO TO 9 END IF DO 4 I=1,NP IF(SRC.EQ.POINT(I)) THEN KS(NT)=I GO TO 5 END IF 4 CONTINUE C INV1TT-03 PAUSE 'Error INV1TT-03: Source name not recognized' STOP 5 CONTINUE DO 6 I=1,NP IF(REC.EQ.POINT(I)) THEN KR(NT)=I GO TO 7 END IF 6 CONTINUE C INV1TT-04 PAUSE 'Error INV1TT-04: Receiver name not recognized' STOP 7 CONTINUE GO TO 3 9 CONTINUE NT=NT-1 CLOSE(LU4) C C....................................................................... C C Computing quantities describing objective prior information: C OPEN(LU5,STATUS='SCRATCH',FORM='UNFORMATTED') OPEN(LU7,FILE=FILE6) WRITE(*,*) C KS(NT+1)=NP+1 KR(NT+1)=NP+1 TFIELD(NT+1)=0. IRAY=0 IWAVE=0 NSUM=IPAR(IPAR(IPAR(2))) DO 12 I=1,NT TTMIN(I)=999999. 12 CONTINUE C C Loop for the files with computed rays 20 CONTINUE FILE1=' ' READ(LU0,*,END=70) FILE1,FILE2,FILE3 IF(FILE1.EQ.' ') THEN GO TO 70 END IF I=INDEX(FILE1,' ') J=INDEX(FILE2,' ') K=INDEX(FILE3,' ') WRITE(*,'(''+Processing: '',3A)') FILE1(1:I),FILE2(1:J), * FILE3(1:K) WRITE(*,*) OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') IF(FILE2.EQ.' ') THEN IU1=0 IU2=LU1 ELSE IU1=LU1 IU2=LU2 OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') END IF OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED',STATUS='OLD') IRAYTT=0 C C Loop for the points of intersection of rays with the surface 30 CONTINUE C Reading the results of the complete ray tracing CALL AP00(IU1,IU2,LU3) IF(IPT.LE.1.AND.IRAYTT.NE.0)THEN C New ray - recording results for the last ray iraytt C loop for field travel times - searching for two-point ray DO 39 IT=1,NT IS=KS(IT) IR=KR(IT) SI=(COOR(1,IS)-XI1)**2+(COOR(2,IS)-XI2)**2 * +(COOR(3,IS)-XI3)**2 RI=(COOR(1,IR)-XI1)**2+(COOR(2,IR)-XI2)**2 * +(COOR(3,IR)-XI3)**2 SE=(COOR(1,IS)-XE1)**2+(COOR(2,IS)-XE2)**2 * +(COOR(3,IS)-XE3)**2 RE=(COOR(1,IR)-XE1)**2+(COOR(2,IR)-XE2)**2 * +(COOR(3,IR)-XE3)**2 IF(SE.LE.SI.AND.RI.LE.RE) THEN C interchanging source and receiver points SI=SE RE=RI IS=KR(IT) IR=KS(IT) END IF IF((SI.LE.DIST2.AND.RE.LE.DIST2)) THEN C Synthetic ray may correspond to the field travel time C check for ray directions near the source C (allowable angle deviation +-30 deg: cosine**2=0.750) C (allowable angle deviation +-15 deg: cosine**2=0.933) TTAUX=(COOR(1,IR)-COOR(1,IS))*(XE1-XI1) * +(COOR(2,IR)-COOR(2,IS))*(XE2-XI2) * +(COOR(3,IR)-COOR(3,IS))*(XE3-XI3) IF(TTAUX.GT.0..AND.TTAUX*TTAUX.GT. * 0.933*((COOR(1,IR)-COOR(1,IS))**2 * +(COOR(2,IR)-COOR(2,IS))**2 * +(COOR(3,IR)-COOR(3,IS))**2) * *((XE1-XI1)**2+(XE2-XI2)**2+(XE3-XI3)**2) ) THEN C Synthetic ray corresponds to the field travel time TTAUX=TT-PI(1)*(COOR(1,IS)-XI1) * -PI(2)*(COOR(2,IS)-XI2) * -PI(3)*(COOR(3,IS)-XI3) * +PE(1)*(COOR(1,IR)-XE1) * +PE(2)*(COOR(2,IR)-XE2) * +PE(3)*(COOR(3,IR)-XE3) IF(TTAUX.LT.TTMIN(IT)) THEN TTMIN(IT)=TTAUX C Possible minimum synthetic travel time SDIST=SQRT(SI) RDIST=SQRT(RE) WRITE(LU5) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM) END IF END IF END IF 39 CONTINUE IF(NT.EQ.0) THEN WRITE(LU5) NT+1,TT,TT,0.,0.,(SUM(I),I=1,NSUM) END IF IRAYTT=0 END IF IF(IWAVE.EQ.0) THEN GO TO 60 END IF C *** for future extensions (selection of two-point rays): C IF(IU1.NE.0) THEN C CALL AP30(IREC) C IF(IREC.EQ.0) THEN C IF(IPT.LE.1.) THEN C WRITE(*,'(''+WAVE:'',I3,'' RAY:'',I4,'' POINT:'', C * I4)') IWAVE,IRAY,IPT C END IF C GO TO 30 C END IF C END IF C *** IF(IPT.EQ.1.OR.MOD(IPT,10).EQ.0) THEN WRITE(*,'(''+WAVE:'',I3,'' RAY:'',I4,'' POINT:'',I4)') * IWAVE,IRAY,IPT END IF IRAYTT=IRAY XI1=YI(3) XI2=YI(4) XI3=YI(5) XE1=Y(3) XE2=Y(4) XE3=Y(5) CALL AP01(TT,TTAUX) CALL AP02(PI,PE) CALL AP29(NSUM,SUM) GO TO 30 C End of the loop for points of intersection of rays with surface C 60 CONTINUE CLOSE(LU1) CLOSE(LU2) GO TO 20 C C All minimum travel times and their derivatives are stored in the C scratch file LU5. C C....................................................................... C C Writing objective prior information: C 70 CONTINUE OPEN(LU4,FILE=FILE4) I=MAX0(INDEX(FILE4,' ')-1,11) WRITE(*,'(''+Generating the output: '',A)') FILE4(1:I) ND=0 DO 71 I=1,NT IF(TTMIN(I).LT.999999.) THEN ND=ND+1 END IF 71 CONTINUE C ND is the number of equations. C REWIND(LU5) WRITE(LU4,'(I9,5I13)') ND,NM WRITE(LU4,'(I9,5I13)') (INDM(I),I=1,NM) WRITE(LU7,'(2A)') ' SOURCE RECEIVER TFIELD ', * 'TFIELD-TT SDIST RDIST TT-CHECKSUM' 73 CONTINUE READ(LU5,END=79) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM) TTMIN(NT+1)=TTAUX IF(TTAUX.EQ.TTMIN(IT)) THEN C Minimum synthetic travel time C C System of linear equations: TDIF=TFIELD(IT)-TTAUX WRITE(LU4,'(I9,4X,6G13.6)') IT,TDIF,TERR(IT),TFIELD(IT) WRITE(LU4,'(6G13.6)') (SUM(INDM(I)),I=1,NM) C C Check sums and log output: IF(VPOWER.NE.0.) THEN TTAUX=0. DO 74 I=1,NM J=INDM(I) IF(IPAR(IPAR(IPAR(1))).LT.J) THEN IF(SUM(J).NE.0.) THEN TTAUX=TTAUX+RPAR(J)*SUM(J) END IF END IF 74 CONTINUE IS=KS(IT) IR=KR(IT) TTAUX=TT+VPOWER*TTAUX WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR), * TFIELD(IT),TDIF,SDIST,RDIST,TTAUX ELSE WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR), * TFIELD(IT),TDIF,SDIST,RDIST END IF END IF GO TO 73 C 79 CONTINUE CLOSE(LU4) CLOSE(LU5) CLOSE(LU7) STOP END C C======================================================================= C INCLUDE 'modelv.for' C modelv.for INCLUDE 'metric.for' C metric.for INCLUDE 'srfc.for' C srfc.for INCLUDE 'parmv.for' C parmv.for INCLUDE 'valv.for' C valv.for INCLUDE 'fitv.for' C fitv.for INCLUDE 'var.for' C var.for INCLUDE 'spsp.for' C spsp.for INCLUDE 'soft.for' C soft.for INCLUDE 'means.for' C means.for INCLUDE 'ap.for' C ap.for INCLUDE 'apvar.for' C apvar.for C C======================================================================= C