C
C Program GRDMIGR to perform common-source Kirchhoff migration C C Version: 5.90 C Date: 2005, June 18 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/klimes.htm C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Names of the input and output files: C MIGR='string'... Name of the input ASCII file containing the C names of files with grid values of travel times and C amplitudes. C Default: MIGR='migr.dat' C GRDSEIS='string'... Name of the input ASCII file with seismograms. C File contains NSEIS*NREC real values. C Default: GRDSEIS='grdseis.out' C GRDMIGR='string'... Name of the output ASCII file containing the C grid values of the migrated image. C Default: GRDNEW='grdmigr.out' C For general description of the files with gridded data refer C to file forms.htm. C Data specifying dimensions of the input and output grids: C N1=positive integer... Number of gridpoints along the X1 axis C (inner loop). C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis C (intermediate loop). C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis C (outer loop). C Default: N3=1 C NSEIS=positive integer... Number of points of a seismogram. C Default: NSEIS=1 C OSEIS=real... Time of the first sample of a seismogram. C Default: OSEIS=0. C DSEIS=real... Sampling interval of a seismogram. C Default: DSEIS=1. C ISRC=positive integer... Index of the line in file MIGR C corresponding to the source. C Default: ISRC=1 C IREC=positive integer... Index of the line in file MIGR C corresponding to the first receiver. C Default: IREC=1 C NREC=positive integer... Number of receivers. C Receivers IREC to IREC+NREC-1 are thus considered. C Default: NREC=1 C TAPER=real... Cosine windowing of seismograms. The width of C cosine windows is TAPER intervals between receivers. C The first cosine window starts at receiver IREC-1, the C second cosine window ends at receiver IREC+NREC. C Default: TAPER=0. C Optional parameters specifying the form of the quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C NUMLIN=positive integer ... Number of reals or integers in one C line of the output file. See the description in file C forms.for. C C C Input formatted file MIGR: C For each surface point, the following line: C (1) 'GRDNUM','GRDTT','GRDAMP','GRDAKI',/ C 'GRDNUM'... Any string. Reserved for future extension. C 'GRDTT'... Gridded travel times from the surface point. C 'GRDAMP'... Gridded ray-theory amplitudes from the surface point. C 'GRDAKI'... Gridded amplitudes, modified for the Kirchhoff C integral, corresponding to the surface point. C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C....................................................................... C C Filenames and parameters: CHARACTER*80 FSEP,FNAMES,FSEIS,FMIGR,FILNUM,FILTT,FILAMP,FILAKI INTEGER LU1,LU2 PARAMETER (LU1=1,LU2=2) C C Input data: INTEGER N1,N2,N3,NSEIS,ISRC,IREC,NREC REAL OSEIS,DSEIS,TAPER C C Memory allocation: INTEGER ISEIS,IMIGR,ITT1,IAMP1,ITT2,IAMP2,NRAM C C Other variables: INTEGER N123,I1,I2,I3,I4,I,IT,ISEIS1 REAL WEIGHT,AMP1,AMP2,A,TT,DT C C----------------------------------------------------------------------- C C Reading input SEP parameter file: WRITE(*,'(A)') '+GRDMIGR: Enter input filename: ' FSEP=' ' READ(*,*) FSEP IF (FSEP.EQ.' ') THEN C GRDMIGR-01 CALL ERROR('GRDMIGR-01: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF CALL RSEP1(LU1,FSEP) WRITE(*,'(A)') '+GRDMIGR: Working ... ' C C Reading input parameters from the SEP file: CALL RSEP3T('MIGR',FNAMES,'migr.dat') CALL RSEP3T('GRDSEIS',FSEIS,'grdseis.out') CALL RSEP3T('GRDMIGR',FMIGR,'grdmigr.out') CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) N123=N1*N2*N3 CALL RSEP3I('NSEIS',NSEIS,1) CALL RSEP3R('OSEIS',OSEIS,0.) CALL RSEP3R('DSEIS',DSEIS,1.) CALL RSEP3I('ISRC',ISRC,1) CALL RSEP3I('IREC',IREC,1) CALL RSEP3I('NREC',NREC,1) CALL RSEP3R('TAPER',TAPER,0.) C C Memory allocation: ISEIS=0 IMIGR=ISEIS+NSEIS*NREC ITT1 =IMIGR+N123 IAMP1=ITT1 +N123 ITT2 =IAMP1+N123 IAMP2=ITT2 +N123 NRAM =IAMP2+N123 IF(NRAM.GT.MRAM) THEN C GRDMIGR-02 CALL ERROR('GRDMIGR-02: Too small array RAM(MRAM)') C Array RAM(MRAM) allocated in include file 'ram.inc' is too small C to contain both input and output grids. C You may wish to increase the dimension MRAM in file C ram.inc. END IF C C Reading the seismograms: CALL RARRAY * (LU2,FSEIS,'FORMATTED',.TRUE.,0.,NSEIS*NREC,RAM(ISEIS+1)) C C Initializing the image: I=0 DO 13 I3=1,N3 DO 12 I2=1,N2 DO 11 I1=1,N1 I=I+1 RAM(IMIGR+I)=0. 11 CONTINUE 12 CONTINUE 13 CONTINUE C C Reading the source grid: OPEN(LU1,FILE=FNAMES,FORM='formatted',STATUS='old') DO 14 I4=1,ISRC READ(LU1,*) FILNUM,FILTT,FILAMP,FILAKI 14 CONTINUE CALL RARRAY(LU2,FILTT ,'FORMATTED',.TRUE.,0.,N123,RAM(ITT1 +1)) CALL RARRAY(LU2,FILAMP,'FORMATTED',.TRUE.,0.,N123,RAM(IAMP1+1)) C C Finding the first receiver grid: REWIND(LU1) DO 15 I4=1,IREC-1 READ(LU1,*) FILNUM,FILTT,FILAMP,FILAKI 15 CONTINUE C C Loop over the receivers: ISEIS1=ISEIS-NSEIS+1 DO 24 I4=1,NREC ISEIS1=ISEIS1+NSEIS READ(LU1,*) FILNUM,FILTT,FILAMP,FILAKI CALL RARRAY(LU2,FILTT ,'FORMATTED',.TRUE.,0.,N123,RAM(ITT2 +1)) CALL RARRAY(LU2,FILAKI,'FORMATTED',.TRUE.,0.,N123,RAM(IAMP2+1)) C C Cosine window: WEIGHT=1. IF(FLOAT(I4).LT.TAPER) THEN WEIGHT=0.5-0.5*COS(3.14159*FLOAT(I4)/TAPER) END IF IF(FLOAT(NREC-I4+1).LT.TAPER) THEN WEIGHT=(0.5-0.5*COS(3.14159*FLOAT(NREC-I4+1)/TAPER))*WEIGHT END IF C C Loop over the depth points I=0 DO 23 I3=1,N3 DO 22 I2=1,N2 DO 21 I1=1,N1 I=I+1 A=0. AMP1=RAM(IAMP1+I) AMP2=RAM(IAMP2+I) IF(AMP1.GT.0..AND.AMP2.GT.0.) THEN TT=(RAM(ITT1+I)+RAM(ITT2+I)-OSEIS)/DSEIS IT=INT(TT) IF(0..LE.TT.AND.IT.LT.NSEIS) THEN DT=TT-FLOAT(IT) IT=ISEIS1+IT A=(1.-DT)*RAM(IT)+DT*RAM(IT+1) A=A*AMP2/AMP1 END IF END IF RAM(IMIGR+I)=RAM(IMIGR+I)+WEIGHT*A 21 CONTINUE 22 CONTINUE 23 CONTINUE C 24 CONTINUE C C Writing the migrated image: CALL WARRAY * (LU2,FMIGR,'FORMATTED',.FALSE.,0.,.FALSE.,0.,N123,RAM(IMIGR+1)) CLOSE(LU1) WRITE(*,'(A,I8,A)') '+GRDMIGR: Done (source',ISRC,'). ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= C