C PROGRAM INV1 TO EVALUATE THE DERIVATIVES OF THE TRAVEL TIME WITH C RESPECT TO THE MODEL COEFFICIENTS. 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 INV1 ASSUMES ALL MODEL PARAMETERS (COEFFICIENTS) STORED IN THE C COMMON BLOCK /VALC/ AS IN THE SUBMITTED VERSIONS OF USER-DEFINED MODEL C SPECIFICATION FORTRAN77 SOURCE CODE FILES 'SRFC.FOR', 'PARM.FOR', AND C 'VAL.FOR'. THUS, UNLIKE THE OTHER PARTS OF THE COMPLETE RAY TRACING, C THE INV1 PROGRAM CANNOT WORK WITH USER'S MODIFICATIONS OF THE C SUBROUTINES SRFC1, SRFC2, PARM1, AND PARM2. C C MAIN INPUT DATA FILE READ FROM THE INTERACTIVE DEVICE (*): C (1) 'INV1' C 'INV1' NAME OF THE INPUT FILE DESCRIBED BELOW. C DEFAULT: 'INV1'='INV1.OUT'. C C INPUT DATA 'INV1': 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) 'POINTS','FTT','MODEL','DATA','SOFT','INV1LOG',/ 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 '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 'MODEL'... INPUT DATA FILE CONTAINING THE MODEL PARAMETERS. C SEE THE FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'. 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 'SOFT'..OUTPUT FILE CONTAINING THE SUBJECTIVE PRIOR INFORMATION. C ITS STRUCTURE IS DESCRIBED BELOW. C IF THE FILENAME IS BLANK, NO FILE WITH SUBJECTIVE PRIOR C INFORMATION IS GENERATED. C 'INV1LOG'... OUTPUT LOG FILE. JUST FOR YOUR INFORMATION. C NOT OPENED AND GENERATED IF THE 'DATA' FILENAME ABOVE IS C BLANK. C /... AN OBLIGATORY SLASH AT THE END OF LINE. C DEFAULT: 'MODEL'='MODEL.DAT', 'DATA'='DATA.OUT', C 'SOFT'='SOFT.OUT', 'INV1LOG'='INV1LOG.OUT'. C (2) DIST,VPOWER,ERRINT,ERRVEL,/ 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 ERRINT..WEIGHTING FACTOR OF THE POSTERIOR INTERFACE B-SPLINE C COEFFICIENT ERROR. C ERRVEL..WEIGHTING FACTOR OF THE POSTERIOR VELOCITY B-SPLINE C COEFFICIENT ERROR. C /... AN OBLIGATORY SLASH AT THE END OF LINE. C DEFAULT: VPOWER=0.0, ERRINT=1.0, ERRVEL=1.0. C (3) ANY TIMES THE FOLLOWING DATA: C 'RAYS','ENDPOINTS','INITIALPOINTS' C 'RAYS'... FILE WITH THE QUANTITIES STORED ALONG RAYS (SEE C C.R.T.5.5.1). C 'ENDPOINTS'... 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 'INITIALPOINTS'... FILE WITH THE QUANTITIES AT THE INITIAL POINTS C OF RAYS, CORRESPONDING TO THE ABOVE FILE RAYS OR ENDPOINTS C (SEE C.R.T.6.1). C (4) / (A SLASH). C (5) (NW1(I),NW2(I),NW3(I),I=1,NW) C LIST OF PARTIAL DERIVATIVES INCLUDED IN THE SOBOLEV SCALAR PRODUCT C WHICH IS ASSUMED TO REPRESENT SUBJECTIVE PRIOR INFORMATION ABOUT C THE MODEL. C NW1,NW2,NW3... ORDERS OF PARTIAL DERIVATIVES WITH RESPECT TO C X1,X2,X3 COORDINATES. C (6) ((WCS(I,J),I=1,J),J=1,NW) C ELEMENTS OF THE CONSTANT SYMMETRIC WEIGHTING MATRIX OF THE SOBOLEV C SCALAR PRODUCT. C WCS(I,J)... COEFFICIENT OF THE PRODUCT OF C (NW1(I),NW2(I),NW3(I))-TH AND (NW1(J),NW2(J),NW3(J))-TH C PARTIAL DERIVATIVES OF FUNCTIONS IN THE SOBOLEV SCALAR C PRODUCT. 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 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 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 OUTPUT FILE 'SOFT': C (1) (INDM(I),I=1,NM) C INDM... INDICES OF THE MODEL COEFFICIENTS INFLUENCING THE TRAVEL C TIMES. C (2) (WM(I),I=1,NM) C RS... PARAMETERS (COEFFICIENTS) OF THE INITIAL MODEL. C (3) (WM(I),I=1,NM) C WM... RELATIVE POSTERIOR ERRORS OF THE MODEL PARAMETERS. C WM=ERRVEL FOR MATERIAL PARAMETERS (SEE (2) OF THE INPUT C DATA 'INV1'), C WM=ERRINT FOR STRUCTURAL INTERFACES (SEE (2) OF THE INPUT C DATA 'INV1'). C DEFAULT: ERRVEL=1.0, ERRINT=1.0. C (4) ((CS(I,J),I=1,J),J=1,NM): C CS... RELATIVE INVERSE SUBJECTIVE PRIOR INFORMATION COVARIANCE C MATRIX. NOTE THAT SUBJECTIVE DATA ARE ASSUMED TO BE MINUS C THE ABOVE INITIAL MODEL PARAMETERS, AND THE MATRIX GS C PROJECTING MODEL PARAMETERS ONTO THE SUBJECTIVE DATA IS C ASSUMED TO BE IDENTITY. C C OUTPUT LOG FILE 'LOG': 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 DATE: 1994, JANUARY 23 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C COMMON BLOCK /VALC/: INTEGER NPAR PARAMETER (NPAR=8191) INTEGER IPAR(0:NPAR) REAL RPAR(0:NPAR) EQUIVALENCE (IPAR,RPAR) COMMON/VALC/IPAR C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /POINTC/: INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF INTEGER ICB1I,IEND,ISHEET REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25) COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF, * NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI C C----------------------------------------------------------------------- C C FILENAMES: CHARACTER*80 FILE0,FILE1,FILE2,FILE3,FILE4,FILE5,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: INTEGER MP,MT PARAMETER (MP=1000) PARAMETER (MT=4000) CHARACTER*1 TEXT(10) CHARACTER*11 POINT(MP),SRC,REC INTEGER NP,NT,KS(MT),KR(MT) REAL DIST,DIST2,VPOWER,ERRINT,ERRVEL REAL COOR(3,MP),TFIELD(MT),TERR(MT) INTEGER MW PARAMETER (MW=12) INTEGER NW,NW1(MW+1),NW2(MW),NW3(MW) REAL WCS(MW*(MW+1)/2) 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 NW... NUMBER OF SPECIFIED PARTIAL DERIVATIVES. C NW1,NW2,NW3... ORDERS OF PARTIAL DERIVATIVES. C WCS(I,J)... COEFFICIENT MATRIX OF THE SOBOLEV SCALAR PRODUCT. C C OUTPUT DATA: VARIATIONS OF THE SYNTHETIC TRAVEL TIME: INTEGER NSUM,NM,INDM(NPAR) REAL SUM(NPAR) C NM... NUMBER OF THE UNKNOWN MODEL PARAMETERS. C INDM... INDICES OF THE UNKNOWN MODEL PARAMETERS. C C OUTPUT DATA: SUBJECTIVE PRIOR INFORMATION: INTEGER MM,MCS C MM... MAXIMUM NUMBER OF MODEL PARAMETERS. PARAMETER (MM=200) PARAMETER (MCS=MM*(MM+1)/2) REAL CS(MCS) C C AUXILIARY STORAGE LOCATIONS: INTEGER IS,IR,IT,ND,IRAYTT,I,J,K REAL TTMIN(MT) REAL TT,TTAUX,TDIF,XI1,XI2,XI3,XE1,XE2,XE3,PI(6),PE(6) REAL SI,SE,RI,RE,SDIST,RDIST,WEIGHT 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 WEIGHT..WEIGHTING FACTOR CORRESPONDING TO THE DERIVATIVES C (AVERAGE QUADRATIC VALUE OF THE DERIVATIVES ON INPUT) C C MORE THAN 23*MP+16*MT+4*MM*(MM+1) BYTES OF MEMORY REQUIRED. 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 'INV1': OPEN(LU0,FILE=FILE0,STATUS='OLD') FILE2='MODEL.DAT' FILE4='DATA.OUT' FILE5='SOFT.OUT' FILE6='INV1LOG.OUT' READ(LU0,*) FILE0,FILE1,FILE2,FILE4,FILE5,FILE6 VPOWER=0. ERRINT=1. ERRVEL=1. READ(LU0,*) DIST,VPOWER,ERRINT,ERRVEL 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 IF(FILE4.EQ.' ') THEN C NO OBJECTIVE INFORMATION FILE GENERATED. 10 CONTINUE FILE1=' ' READ(LU0,*,END=70) FILE1,FILE2,FILE3 IF(FILE1.EQ.' ') THEN C GO TO (SUBJECTIVE PRIOR INFORMATION) GO TO 80 END IF GO TO 10 END IF 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 PAUSE 'ERROR: 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 PAUSE 'ERROR: 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 PAUSE 'ERROR: 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 PAUSE 'ERROR: 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 ENDPOINTS OF RAYS 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 ENDPOINTS OF RAYS 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) C C....................................................................... C C SUBJECTIVE PRIOR INFORMATION: C 80 CONTINUE IF(FILE5.NE.' ') THEN C IF(NM*(NM+1)/2.GT.MCS) THEN PAUSE 'ERROR: TOO MANY MODEL PARAMETERS' STOP END IF C C OPENING OUTPUT FILE: OPEN(LU4,FILE=FILE5) I=MAX0(INDEX(FILE5,' ')-1,11) WRITE(*,'('' GENERATING THE OUTPUT: '',A)') FILE5(1:I) C C READING PRIOR SUBJECTIVE INFORMATION COEFFICIENTS: DO 81 I=1,MW+1 NW1(I)=-1 81 CONTINUE READ(LU0,*) (NW1(I),NW2(I),NW3(I),I=1,MW),NW1(MW+1) DO 82 I=1,MW IF(NW1(I).EQ.-1) THEN NW=I-1 GO TO 83 END IF 82 CONTINUE PAUSE 'ERROR: TOO LARGE ENERGY COEFFICIENT MATRIX' STOP 83 CONTINUE READ(LU0,*) (WCS(I),I=1,NW*(NW+1)/2) C C GENERATING PRIOR SUBJECTIVE INFORMATION COVARIANCE MATRIX: DO 84 I=1,MCS CS(I)=0. 84 CONTINUE DO 89 J=1,NW DO 88 I=1,J WEIGHT=WCS(J*(J-1)/2+I) IF(WEIGHT.NE.0.) THEN CALL SOFT(2,NW1(I),NW2(I),NW3(I),NW1(J),NW2(J),NW3(J), * WEIGHT,NM,INDM,CS) END IF 88 CONTINUE 89 CONTINUE C C (1,2) WRITING THE INITIAL MODEL PARAMETERS: WRITE(LU4,'(I9,5I13)') (INDM(I),I=1,NM) WRITE(LU4,'(6G13.6)') (RPAR(INDM(I)),I=1,NM) C C (3) WRITING RELATIVE POSTERIOR ERRORS OF THE MODEL PARAMETERS: DO 93 I=1,NM J=INDM(I) IF(IPAR(IPAR(IPAR(1))).LT.J) THEN SUM(I)=ERRVEL ELSE SUM(I)=ERRINT END IF 93 CONTINUE WRITE(LU4,'(24F3.0)') (SUM(I),I=1,NM) C C (4) WRITING THE PRIOR SUBJECTIVE INFORMATION COVARIANCE MATRIX: DO 94 J=1,NM WRITE(LU4,'(6G13.6)') (CS(I),I=J*(J-1)/2+1,J*(J+1)/2) 94 CONTINUE CLOSE(LU4) C END IF C STOP END C C======================================================================= C SUBROUTINE SOFT(NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,WEIGHT,NM,INDM,CS) C INTEGER NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,NM,INDM(*) REAL WEIGHT,CS(*) C C THIS SUBROUTINE ACCUMULATES THE PRIOR SUBJECTIVE INFORMATION C COVARIANCE MATRIX DESCRIBING THE SMOOTHNESS OF THE FUNCTIONS C INTERPOLATED BY MEANS OF SUBROUTINES OF THE FILE 'VAL.FOR'. C C INPUT: C NCLASS..NUMBER OF CLASSES REQUIRED. C NA1,NA2,NA3... ORDERS OF X1,X2,X3-PARTIAL DERIVATIVES OF THE FIRST C BASIS FUNCTION IN THE PRODUCT BEING INTEGRATED. C NB1,NB2,NB3... ORDERS OF X1,X2,X3-PARTIAL DERIVATIVES OF THE C SECOND BASIS FUNCTION IN THE PRODUCT BEING INTEGRATED. C WEIGHT..WEIGHTING FACTOR CORRESPONDING TO THE DERIVATIVES C (AVERAGE QUADRATIC VALUE OF THE DERIVATIVES ON INPUT) C CS... SYMMETRIC MATRIX TO BE REGULARIZED: C CONTRIBUTION CORRESPONDING TO THE DERIVATIVES C IS ADDED TO THE MATRIX CS C NONE OF THE INPUT PARAMETERS ARE ALTERED. C C OUTPUT: C NM... NUMBER OF THE UNKNOWN MODEL PARAMETERS. C INDM... INDICES OF THE UNKNOWN MODEL PARAMETERS. C CS... RESULTING SYMMETRIC SUBJECTIVE PRIOR INFORMATION C COVARIANCE MATRIX. C C COMMON BLOCK: INTEGER NPAR PARAMETER (NPAR=8191) INTEGER IPAR(0:NPAR) REAL RPAR(0:NPAR) EQUIVALENCE (IPAR,RPAR) COMMON/VALC/IPAR C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL SPSP C C ERROR MESSAGES: C 359A.. INCORRECT NUMBER OF CLASSES: C THE NUMBER OF CLASSES REQUIRED IS ZERO, NEGATIVE OR C GREATER THAN THE NUMBER OF CLASSES DEFINED. C 359B.. INSUFFICIENT MEMORY: C THE DIMENSION OF THE BUFFER B IS NOT SUFFICIENT. C C DATE: 1992, SEPTEMBER 7 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C INTEGER MM,MB C MM... MAXIMUM NUMBER OF MODEL PARAMETERS. PARAMETER (MM=200) PARAMETER (MB=MM*(MM+1)/2+3*(20*21/2)) REAL B(MB) INTEGER ICLASS,JGROUP,LFUNCT,MFUNCT,JFUNCT,LADR,MADR,IADR,NVAR INTEGER NX(3),NX1,NX2,NX3 EQUIVALENCE (NX(1),NX1),(NX(2),NX2),(NX(3),NX3) INTEGER JADR(7),JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7 EQUIVALENCE (JADR(1),JADR1),(JADR(2),JADR2),(JADR(3),JADR3) EQUIVALENCE (JADR(4),JADR4),(JADR(5),JADR5),(JADR(6),JADR6) EQUIVALENCE (JADR(7),JADR7) INTEGER I1,I2,I3,I,J,N REAL WDENS,SIGMA C C B... AUXILIARY STORAGE LOCATIONS FOR PARTIAL COVARIANCE C MATRICES. C ICLASS..INDEX OF THE CLASS. C JGROUP..ADDRESS OF A CURRENT GROUP. C LFUNCT,MFUNCT,JFUNCT... ADDRESSES OF THE FIRST, LAST AND ARBITRARY C FUNCTIONS OF THE GROUP. C LADR,MADR,IADR... ADDRESSES OF THE FIRST, LAST AND ARBITRARY C PARAMETERS OF THE CURRENT FUNCTION. C NVAR... NUMBER OF THE INDEPENDENT VARIABLES A1, A2, A3 OF THE C INTERPOLATED FUNCTION W. C NX=(NX1,NX2,NX3)... NUMBERS OF GRID LINES. C JADR=(JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7)... ADDRESSES OF C PARAMETERS DESCRIBING THE INTERPOLATED FUNCTION (GRID C COORDINATES, B-SPLINE COEFFICIENTS, B-SPLINE BASIS C FUNCTIONS). C I1,I2,I3,I,J,N... LOCAL AUXILIARY VARIABLES. C WDENS...WEIGHTING FACTOR DIVIDED BY VOLUME (AREA,LENGTH). C SIGMA...TENSION FACTOR. C C....................................................................... C NM=0 IF(NCLASS.LT.1.OR.IPAR(0).LT.NCLASS) THEN PAUSE 'ERROR 359A IN SOFT: INCORRECT NUMBER OF CLASSES' STOP END IF C LOOP FOR CLASSES: DO 92 ICLASS=1,NCLASS C LOOP FOR GROUPS OF THE CURRENT CLASS: DO 91 JGROUP=IPAR(ICLASS-1)+1,IPAR(ICLASS) LFUNCT=IPAR(JGROUP-1)+1 MFUNCT=IPAR(JGROUP) MADR =IPAR(LFUNCT-1) C C LOOP FOR FUNCTIONS OF THE CURRENT GROUP: DO 90 JFUNCT=LFUNCT,MFUNCT C STARTING AND END ADDRESSES OF THE PARAMETERS DESCRIBING THE C FUNCTION LADR=MADR+1 MADR=IPAR(JFUNCT) IF(LADR.LE.MADR) THEN WDENS=WEIGHT C TENSION FACTOR SIGMA=RPAR(LADR+5) C C THE NUMBER, TYPES AND VALUES OF THE INDEPENDENT VARIABLES C OF THE INTERPOLATED FUNCTION: C INITIAL ADDRESS IADR=LADR+6 C INITIAL NUMBER OF THE INDEPENDENT VARIABLES NVAR=0 NX1=1 NX2=1 NX3=1 JADR1=0 JADR2=0 JADR3=0 JADR4=0 C LOOP FOR THE POSSIBLE INDEPENDENT VARIABLES: DO 20 I=LADR+2,LADR+4 C TYPE OF THE POSSIBLE INDEPENDENT VARIABLE: J=IPAR(I) IF(J.GT.0) THEN N=IABS(IPAR(IADR)) IF(N.GE.2) THEN NVAR=NVAR+1 NX(NVAR)=N ELSE IF(N.EQ.1) THEN JADR(NVAR+1)=JADR(NVAR+1)+1 END IF IADR=IADR+1 END IF 20 CONTINUE JADR5=0 JADR6=0 JADR7=0 C C INTERPOLATED FUNCTION W: JADR1=IADR+JADR1 IF(NVAR.LE.0) THEN C NO INDEPENDENT VARIABLE. JADR4=JADR1 ELSE JADR2=JADR1+NX1+JADR2 WDENS=WDENS/(RPAR(JADR2-1)-RPAR(JADR1)) IF(NVAR.EQ.1) THEN C ONE INDEPENDENT VARIABLE: JADR4=JADR2 JADR5=JADR2+NX1 ELSE JADR3=JADR2+NX2+JADR3 WDENS=WDENS/(RPAR(JADR3-1)-RPAR(JADR2)) IF(NVAR.EQ.2) THEN C TWO INDEPENDENT VARIABLES: JADR4=JADR3 JADR5=JADR3+NX1*NX2 JADR6=JADR4+5*NX1 ELSE C THREE INDEPENDENT VARIABLES: JADR4=JADR3+NX3+JADR4 JADR5=JADR4+NX1*NX2*NX3 JADR6=JADR5+5*NX1 JADR7=JADR6+5*NX2 WDENS=WDENS/(RPAR(JADR4-1)-RPAR(JADR3)) END IF END IF END IF C WDENS IS WEIGHTING FACTOR DIVIDED BY VOLUME (AREA,LENGTH). C JADR4=JADR4-1 N=NX1*NX2*NX3 IF(WEIGHT.EQ.0.) THEN DO 72 I2=1,N INDM(NM+I2)=JADR4+I2 72 CONTINUE NM=NM+N ELSE IF(NA1.EQ.NB1.AND.NA2.EQ.NB2.AND.NA3.EQ.NB3) THEN I=N*(N+1)/2 J=0 ELSE I=N*N J=1 END IF IF(NA1.EQ.NB1) THEN I1=NX1*(NX1+1)/2 ELSE I1=NX1*NX1 END IF IF(NA2.EQ.NB2) THEN I2=NX2*(NX2+1)/2 ELSE I2=NX2*NX2 END IF IF(NA3.EQ.NB3) THEN I3=NX3*(NX3+1)/2 ELSE I3=NX3*NX3 END IF I1=I1+I I2=I2+I1 I3=I3+I2 IF(I3.GT.MB) THEN PAUSE 'ERROR 359B IN SOFT: INSUFFICIENT MEMORY' STOP END IF CALL SPSP(NA1,NA2,NA3,NB1,NB2,NB3,NX1,NX2,NX3, * RPAR(JADR1),RPAR(JADR2),RPAR(JADR3), * RPAR(JADR5),RPAR(JADR6),RPAR(JADR7), * SIGMA,WDENS,B(1),B(I+1),B(I1+1),B(I2+1)) I=NM*(NM+1)/2 DO 82 I2=1,N INDM(NM+I2)=JADR4+I2 I=I+NM DO 81 I1=1,I2 I=I+1 IF(J.EQ.0) THEN CS(I)=CS(I)+B(I2*(I2-1)/2+I1) ELSE CS(I)=CS(I)+B(N*(I2-1)+I1)+B(N*(I1-1)+I2) END IF 81 CONTINUE 82 CONTINUE NM=NM+N C CONTRIBUTION CORRESPONDING TO THE DERIVATIVES C IS ADDED TO THE MATRIX CS. END IF C END IF 90 CONTINUE 91 CONTINUE 92 CONTINUE C END OF LOOPS FOR FUNCTIONS. C RETURN END C C======================================================================= C