C
C Program GSECAL calculates a linear combination of the given
C seismograms stored in the GSE data exchange format.
C
C Version: 6.50
C Date: 2011, May 19
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 GSE files and corresponding coefficients:
C     Here, # represents the value of integer constant MFILSS specified
C             in the this program.
C             E.g., if MFILSS=12, SS# stands for SS12, COEF# stands for
C             COEF12, etc.
C     SS1='string', SS2='string', ..., SS#='string'... Strings with the
C             names of the input data files in the GSE data exchange
C             format, containing the given seismograms.
C             For the GSE data exchange format, refer to the description
C             in file 'gse.for'.
C             Just nonblank filenames are considered.
C             Defaults: SS1=' ', SS2=' ', SS3=' ', ..., SS#=' '
C     COEF1=real, COEF2=real, COEF3=real, ..., COEF#=real:
C             Coefficients of the linear combination to multiply
C             seismograms from the corresponding input GSE files.
C             Defaults: COEF1=1., COEF2=1., COEF3=1., ..., COEF#=1.
C Name of the output GSE file:
C     SS='string'... String with the name of the output file in the GSE
C             data exchange format, containing the specified linear
C             combination of the input seismograms.
C             For the GSE data exchange format, refer to the description
C             in file 'gse.for'.
C             Default: SS='ss.gse'
C Optional parameters:
C     TMIN=real... Times of the samples in the output seismograms
C             will differ from TMIN by integer multiples of DT.
C             Default: TMIN=0. (mostly sufficient)
C     DT=real... Time interval to digitize the output seismograms.
C             If DT=0. (default), time TMIN and time interval DT are
C               taken, for each component at each receiver, from the
C               first nonzero seismogram
C             Default: DT=0.
C
C=======================================================================
C
C     External functions and subroutines:
      EXTERNAL ERROR
C     ERROR
      EXTERNAL RSEP1,RSEP3T,RSEP3R
C     RSEP1,RSEP3T,RSEP3R... Subroutines of the Fortran 77 file
C             sep.for (package FORMS).
      EXTERNAL RGSE1,RGSE2,WGSE1,WGSE2,WGSE3
C     RGSE1,RGSE2,WGSE1,WGSE2,WGSE3... Subroutines of the Fortran 77
C             file gse.for (package FORMS)
C             to store seismograms in the GSE data exchange format.
      EXTERNAL LENGTH
      INTEGER  LENGTH
C     LENGTH
      EXTERNAL STRIND
C     STRIND
C     The length of character function STRIND is declared later on.
C
C.......................................................................
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C.......................................................................
C
C     Parameters:
      INTEGER LU,MFILSS,MDIGSS,MREC
      PARAMETER (LU=1,MFILSS=29,MDIGSS=2,MREC=2048)
C     LU...   Logical unit number
C     MFILSS..Maximum number of input files with synthetic seismograms.
C     MDIGSS..Number of digits of MFILSS.
C     MREC... Maximum number of receivers.
C
C     Length of the external character function:
      CHARACTER*(4+MDIGSS) STRIND
C
C     Filenames:
      CHARACTER*80 FSEP,FILESS
C
C     Coefficient of the linear combination corresponding to the file:
      REAL COEF
C
C     Storing seismograms in memory
      CHARACTER*6 REC(MREC)
      COMMON/RECC/REC
      INTEGER KSS(0:3*MREC+1)
      REAL TMIN(0:3*MREC),DT(0:3*MREC)
      REAL X1(MREC),X2(MREC),X3(MREC)
      EQUIVALENCE (KSS(0),IRAM(1))
      EQUIVALENCE (TMIN(0),IRAM(3*MREC+3))
      EQUIVALENCE (DT(0),IRAM(6*MREC+4))
      EQUIVALENCE (X1(1),IRAM(9*MREC+5))
      EQUIVALENCE (X2(1),IRAM(10*MREC+5))
      EQUIVALENCE (X3(1),IRAM(11*MREC+5))
C     REC(I)... Name of the I-th receiver.
C     KSS(3*I-3+K)... Index of the last sample of K-th seismogram
C             component at the I-th receiver.
C     KSS(3*MREC+1)... Index of the last sample of the current input
C             seismogram
C     TMIN(3*I-3+K)... Time of the first seismogram sample at the I-th
C             receiver.
C     DT(3*I-3+K)... Digitization interval of the seismogram.
C     X1(I),X2(I),X3(I)... Coordinates of the I-th receiver.
C
C     Parameters of the seismogram:
      CHARACTER*1  CHAN,TEXT
      CHARACTER*7  NAMREC
      INTEGER KOMP,NSS
      REAL X1R,X2R,X3R,T0,TD
C
C     Other variables:
      INTEGER IFILSS,NREC,IREC,J,ILEFT,IRIGHT,ISHIFT
C     IFILSS..Index of the input file.
C     NREC... Number of already identified receivers.
C     IREC... Index of the receiver.
C     J...    Index of the seismogram.
C     ILEFT,IRIGHT,ISHIFT... Extending memory for the new seismogram.
C
C     Temporary storage locations:
      INTEGER IT,I
      REAL T,RT,S
C
C     Occupied array RAM:
      KSS(0)=12*MREC+4
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GSECAL: Enter input filename: '
      FSEP=' '
      READ(*,*) FSEP
      WRITE(*,'(A)') '+GSECAL: Working...            '
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.NE.' ') THEN
        CALL RSEP1(LU,FSEP)
      ELSE
C       GSECAL-01
        CALL ERROR('GSECAL-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.
      END IF
C
C     Reading optional input data:
      CALL RSEP3R('TMIN',TMIN(0),0.)
      CALL RSEP3R('DT',DT(0),0.)
C
C     Initial values:
C     No receiver
      NREC=0
C     No seismogram
      DO 11 I=1,3*MREC+1
        KSS(I)=KSS(0)
   11 CONTINUE
C
C.......................................................................
C
C     Loop for input GSE files:
      DO 88 IFILSS=1,MFILSS
        CALL RSEP3T(STRIND('SS',IFILSS),FILESS,' ')
        CALL RSEP3R(STRIND('COEF',IFILSS),COEF,1.)
        IF(FILESS.NE.' ') THEN
C
C         Opening the input GSE file
          WRITE(*,'(2A)') '+GSECAL: Reading ',FILESS(1:63)
          OPEN(LU,FILE=FILESS,STATUS='OLD')
          CALL RGSE1(LU,TEXT)
C
C         Loop for input seismograms
   20     CONTINUE
            I=KSS(3*MREC)
            CALL RGSE2(LU,NAMREC,CHAN,KOMP,X1R,X2R,X3R,T0,TD,
     *                                              NSS,MRAM-I,RAM(I+1))
            IF(NSS.LE.-1) THEN
C             End of the GSE file
              GO TO 87
            END IF
            KSS(3*MREC+1)=KSS(3*MREC)+NSS
            IF(KOMP.LT.1.OR.KOMP.GT.3) THEN
C             GSECAL-02
              CALL ERROR('GSECAL-02: Component other than 1, 2 or 3')
C             Input GSE file contains seismogram component other than
C             1, 2 or 3, which is not allowed by this version of this
C             program.
            END IF
            IF(NSS.GE.1) THEN
C
C             Determining the index of the receiver
              DO 31 IREC=1,NREC
                IF(REC(IREC).EQ.NAMREC) THEN
C                 Already identified receiver
                  IF(X1(IREC).NE.X1R.OR.
     *               X2(IREC).NE.X2R.OR.
     *               X3(IREC).NE.X3R) THEN
C                   GSECAL-03
                    CALL ERROR('GSECAL-03: Different coordinates')
C                   The same receiver has different coordinates in
C                   different input GSE files.
                  END IF
                  GO TO 32
                END IF
   31         CONTINUE
                IF(NREC.GE.MREC) THEN
C                 GSECAL-04
                  CALL ERROR('GSECAL-04: Too many receivers')
C                 Number of receivers exceeds integer constant MREC
C                 specified in the this program.
C                 You may wish to increase parameter MREC.
                END IF
                NREC=NREC+1
                REC(NREC)=NAMREC
                IREC=NREC
                X1(IREC)=X1R
                X2(IREC)=X2R
                X3(IREC)=X3R
   32         CONTINUE
C             Index of the seismogram
              J=3*IREC-3+KOMP
C
C             Determining the memory for the new seismogram
              IF(KSS(J).EQ.KSS(J-1)) THEN
C               Zero stored seismogram
                IF(DT(0).EQ.0.) THEN
                  TMIN(J)=T0
                  DT(J)=TD
                  ILEFT=0
                  IRIGHT=NSS
                ELSE
                  TMIN(J)=TMIN(0)
                  DT(J)=DT(0)
                  ILEFT=NINT((TMIN(J)-T0+TD)/DT(J)-0.5)
                  TMIN(J)=TMIN(J)-FLOAT(ILEFT)*DT(J)
                  ILEFT=0
                  IRIGHT=INT((T0+FLOAT(NSS)*TD-TMIN(J))/DT(J))+1
                END IF
              ELSE
C               Non-zero stored seismogram
                ILEFT=INT((TMIN(J)-T0+TD)/DT(J))
                IRIGHT=INT((T0+FLOAT(NSS)*TD-TMIN(J))/DT(J))
     *                                              -(KSS(J)-KSS(J-1)-1)
                ILEFT=MAX0(0,ILEFT)
                IRIGHT=MAX0(0,IRIGHT)
              END IF
C             ILEFT is the number of new samples before the seismogram
C             IRIGHT is the number of new samples after the seismogram
C
C             Preparing the memory for the new seismogram
              ISHIFT=ILEFT+IRIGHT
              IF(KSS(3*MREC+1)+ISHIFT.GT.MRAM) THEN
C               GSECAL-05
                CALL ERROR('GSECAL-05: Too small array RAM(MRAM)')
C               Array RAM(MRAM) allocated in include file 'ram.inc'
C               is too small to contain the samples of the calculated
C               linear combination of the given seismograms and the
C               samples of the currently read seismogram.
C               You may wish to increase the dimension MRAM in file
C               ram.inc.
              END IF
              DO 41 I=J,3*MREC+1
                KSS(I)=KSS(I)+ISHIFT
   41         CONTINUE
              DO 42 I=KSS(3*MREC+1),KSS(J)+1,-1
                RAM(I)=RAM(I-ISHIFT)
   42         CONTINUE
              DO 43 I=KSS(J),KSS(J)-IRIGHT+1,-1
                RAM(I)=0.
   43         CONTINUE
              DO 44 I=KSS(J)-IRIGHT,KSS(J-1)+ILEFT+1,-1
                RAM(I)=RAM(I-ILEFT)
   44         CONTINUE
              DO 45 I=KSS(J-1)+ILEFT,KSS(J-1)+1,-1
                RAM(I)=0.
   45         CONTINUE
C
C             Loop over the samples of the stored seismogram
              DO 60 I=KSS(J-1)+1,KSS(J)
C               Time of the stored sample
                T=TMIN(J)+FLOAT(I-KSS(J-1)-1)*DT(J)
C               Number of intervals from the first input sample
                RT=(T-T0)/TD
C               Integer part of the number of intervals
                IT=INT(RT+2.)-2
C               Fractional part of the number of intervals
                RT=RT-FLOAT(IT)
                IT=IT+1
                IF(0.LE.IT.AND.IT.LE.NSS) THEN
C                 Linear interpolation of the input seismogram
                  IF(IT.EQ.0) THEN
                    S=(1.-RT)*RAM(KSS(J-1)+IT+1)
                  ELSE IF(IT.EQ.NSS) THEN
                    S=RT*RAM(KSS(J-1)+IT)
                  ELSE
                    S=RT*RAM(KSS(J-1)+IT)+(1.-RT)*RAM(KSS(J-1)+IT+1)
                  END IF
C                 Adding the seismogram to the memory
                  RAM(I)=RAM(I)+COEF*S
                END IF
   60         CONTINUE
C
            END IF
          GO TO 20
C
C         End of loop for input seismograms
   87     CONTINUE
          CLOSE(LU)
C
        END IF
   88 CONTINUE
C     End of loop for input GSE files
C
C.......................................................................
C
C     Writing the output GSE file:
      CALL RSEP3T('SS',FILESS,'ss.gse')
      OPEN(LU,FILE=FILESS)
      CALL WGSE1(LU,' ')
      DO 92 IREC=1,NREC
        DO 91 KOMP=1,3
          I=3*IREC-3+KOMP
          CALL WGSE2(LU,REC(IREC),' ',KOMP,X1(IREC),X2(IREC),X3(IREC),
     *                    TMIN(I),DT(I),KSS(I)-KSS(I-1),RAM(KSS(I-1)+1))
   91   CONTINUE
   92 CONTINUE
      CALL WGSE3(LU)
      CLOSE(LU)
C
      WRITE(*,'(A)') '+GSECAL: Done.                                   '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'gse.for'
C     gse.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C