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     E-mail: klimes@seis.karlov.mff.cuni.cz
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
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