C
C Program GSECAL calculates a linear combination of the given
C seismograms stored in the GSE data exchange format.
C
C Version: 7.40
C Date: 2017, 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 zero or nonzero input seismogram.
C             Default: DT=0.
C     GSEWIDTH=positive integer... Width of the output field reserved
C             for one integer value of the seismogram. Refer to the
C             description in file 'gse.for'.
C Value of undefined quantities:
C     UNDEF=real... The value to be used for undefined real quantities.
C             Default: UNDEF=undefined value used in forms.for
C
C=======================================================================
C
C     External functions and subroutines:
      EXTERNAL ERROR
C     ERROR
      EXTERNAL RSEP1,RSEP2,RSEP3T,RSEP3R,WSEPT,WSEPR,SSEP
C     RSEP1,RSEP2,RSEP3T,RSEP3R,WSEPT,WSEPR,SSEP... Subroutines of the
C             Fortran 77 file sep.for (package FORMS).
      EXTERNAL RGSE1,RGSE2,RGSE2C,WGSE1,WGSE2D,WGSE2C,WGSE2,WGSE3
C     RGSE1,RGSE2,RGSE2C,WGSE1,WGSE2D,WGSE2C,WGSE2,WGSE3... Subroutines
C             of the Fortran 77 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.
      EXTERNAL UARRAY
      REAL UARRAY
C     UARRAY
C
C.......................................................................
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
C     Parameters:
      INTEGER LU,MFILSS,MDIGSS,MREC
      REAL UNDEF
      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     UNDEF...Undefined value.  Determined later on.
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),SRC(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)
      REAL X1SRC(MREC),X2SRC(MREC),X3SRC(MREC)
      EQUIVALENCE (KSS(0),RAM(1))
      EQUIVALENCE (TMIN(0),RAM(3*MREC+3))
      EQUIVALENCE (DT(0),RAM(6*MREC+4))
      EQUIVALENCE (X1(1),RAM(9*MREC+5))
      EQUIVALENCE (X2(1),RAM(10*MREC+5))
      EQUIVALENCE (X3(1),RAM(11*MREC+5))
      EQUIVALENCE (X1SRC(1),RAM(12*MREC+5))
      EQUIVALENCE (X2SRC(1),RAM(13*MREC+5))
      EQUIVALENCE (X3SRC(1),RAM(14*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     X1SRC(I),X2SRC(I),X3SRC(I)... Coordinates of the point source
C             corresponding to the I-th receiver.
C
C     Parameters of the seismogram:
      CHARACTER*1  CHAN,TEXT
      CHARACTER*6  NAMREC,NAMSRC
      INTEGER KOMP,NSS
      REAL X1R,X2R,X3R,X1S,X2S,X3S,T0,TD
C
C     Line of optional SEP parameter file or comments in the GSE file
      CHARACTER*80 LINE
C     Indices of the sets of SEP parameters
      INTEGER ISEP,IOLD
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,ITRAM,I
      REAL T,RT,S
C
C     Occupied array RAM:
      KSS(0)=15*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     Undefined value:
      UNDEF=UARRAY()
C
C     Defining index ISEP of the working set of SEP parameters:
      ISEP=0
      CALL SSEP(ISEP,IOLD)
      CALL SSEP(IOLD,ISEP)
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     Zero seismograms
      DO 11 I=1,3*MREC+1
        KSS(I)=KSS(0)
   11 CONTINUE
C     No seismogram defined
      DO 12 I=1,3*MREC
        DT(I)=0.
   12 CONTINUE
C     No parameters from the seismogram comment lines
      DO 13 I=1,MREC
        SRC(I)=' '
        X1SRC(I)=UNDEF
        X2SRC(I)=UNDEF
        X3SRC(I)=UNDEF
   13 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
C
C           Determining the index of the receiver
            DO 31 IREC=1,NREC
              IF(REC(IREC).EQ.NAMREC) THEN
C               Already identified receiver
                IF(ABS(X1(IREC)+999.).LT.0.001) THEN
                  X1(IREC)=X1R
                END IF
                IF(ABS(X2(IREC)+999.).LT.0.001) THEN
                  X2(IREC)=X2R
                END IF
                IF(ABS(X3(IREC)+999.).LT.0.001) THEN
                  X3(IREC)=X3R
                END IF
                IF(ABS(X1R+999.).LT.0.001) THEN
                  X1R=X1(IREC)
                END IF
                IF(ABS(X2R+999.).LT.0.001) THEN
                  X2R=X2(IREC)
                END IF
                IF(ABS(X3R+999.).LT.0.001) THEN
                  X3R=X3(IREC)
                END IF
                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           Reading the comment lines corresponding to the seismogram
            ISEP=-ISEP
            CALL SSEP(ISEP,IOLD)
   21       CONTINUE
              CALL RGSE2C(LINE,*22)
              CALL RSEP2(LINE)
            GO TO 21
   22       CONTINUE
            CALL RSEP3T('NAMESRC',NAMSRC,' ')
            CALL RSEP3R('X1SRC',X1S,UNDEF)
            CALL RSEP3R('X2SRC',X2S,UNDEF)
            CALL RSEP3R('X3SRC',X3S,UNDEF)
            CALL SSEP(IOLD,ISEP)
            IF(NAMSRC.NE.' ') THEN
              IF(SRC(IREC).EQ.' ') THEN
                SRC(IREC)=NAMSRC
              ELSE IF(SRC(IREC).NE.NAMSRC) THEN
C               GSECAL-05
                CALL ERROR('GSECAL-05: Different source name')
C               The source corresponding to the same receiver has
C               different name in different input GSE files.
              END IF
            END IF
            IF(X1S.NE.UNDEF) THEN
              IF(X1SRC(IREC).EQ.UNDEF) THEN
                X1SRC(IREC)=X1S
              ELSE IF(X1SRC(IREC).NE.X1S) THEN
C               GSECAL-06
                CALL ERROR('GSECAL-06: Different source coordinates')
C               The source corresponding to the same receiver has
C               different X1 coordinate in different input GSE files.
              END IF
            END IF
            IF(X2S.NE.UNDEF) THEN
              IF(X2SRC(IREC).EQ.UNDEF) THEN
                X2SRC(IREC)=X2S
              ELSE IF(X2SRC(IREC).NE.X2S) THEN
C               GSECAL-07
                CALL ERROR('GSECAL-07: Different source coordinates')
C               The source corresponding to the same receiver has
C               different X2 coordinate in different input GSE files.
              END IF
            END IF
            IF(X3S.NE.UNDEF) THEN
              IF(X3SRC(IREC).EQ.UNDEF) THEN
                X3SRC(IREC)=X3S
              ELSE IF(X3SRC(IREC).NE.X3S) THEN
C               GSECAL-08
                CALL ERROR('GSECAL-08: Different source coordinates')
C               The source corresponding to the same receiver has
C               different X3 coordinate in different input GSE files.
              END IF
            END IF
C
C           Defining a new seismogram
            IF(DT(J).EQ.0.) THEN
              IF(DT(0).EQ.0.) THEN
                TMIN(J)=T0
                DT(J)=TD
              ELSE
                TMIN(J)=TMIN(0)
                DT(J)=DT(0)
              END IF
            END IF
C
            IF(NSS.GE.1) THEN
C
C             Determining the memory for the new seismogram
              ILEFT=NINT((TMIN(J)-T0+TD)/DT(J)-0.501)
              IF(KSS(J).EQ.KSS(J-1)) THEN
C               Zero stored seismogram
                TMIN(J)=TMIN(J)-FLOAT(ILEFT)*DT(J)
                ILEFT=0
                IRIGHT=NINT((T0+FLOAT(NSS)*TD-TMIN(J))/DT(J)+0.499)
              ELSE
C               Non-zero stored seismogram
                IRIGHT=NINT((T0+FLOAT(NSS)*TD-TMIN(J))/DT(J)+0.499)
     *                                                -(KSS(J)-KSS(J-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
              TMIN(J)=TMIN(J)-FLOAT(ILEFT)*DT(J)
              ISHIFT=ILEFT+IRIGHT
              IF(KSS(3*MREC+1)+ISHIFT.GT.MRAM) THEN
C               GSECAL-09
                CALL ERROR('GSECAL-09: 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)
C               Index of the left-hand sample of the interval
                IT=IT+1
C               Index of the left-hand sample of the interval in RAM
                ITRAM=KSS(3*MREC)+IT
                IF(0.LE.IT.AND.IT.LE.NSS) THEN
C                 Linear interpolation of the input seismogram
                  IF(IT.EQ.0) THEN
                    S=RT*RAM(ITRAM+1)
                  ELSE IF(IT.EQ.NSS) THEN
                    S=(1.-RT)*RAM(ITRAM)
                  ELSE
                    S=(1.-RT)*RAM(ITRAM)+RT*RAM(ITRAM+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
          J=3*IREC-3+KOMP
          IF(DT(J).NE.0.) THEN
            CALL WGSE2D
            IF(SRC(IREC).NE.' ') THEN
              CALL WSEPT(LINE,'NAMESRC',SRC(IREC))
              CALL WGSE2C(LINE)
            END IF
            IF(X1SRC(IREC).NE.UNDEF) THEN
              CALL WSEPR(LINE,'X1SRC',X1SRC(IREC))
              CALL WGSE2C(LINE)
            END IF
            IF(X2SRC(IREC).NE.UNDEF) THEN
              CALL WSEPR(LINE,'X2SRC',X2SRC(IREC))
              CALL WGSE2C(LINE)
            END IF
            IF(X3SRC(IREC).NE.UNDEF) THEN
              CALL WSEPR(LINE,'X3SRC',X3SRC(IREC))
              CALL WGSE2C(LINE)
            END IF
            CALL WGSE2(LU,REC(IREC),' ',KOMP,X1(IREC),X2(IREC),X3(IREC),
     *                    TMIN(J),DT(J),KSS(J)-KSS(J-1),RAM(KSS(J-1)+1))
          END IF
   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
      INCLUDE 'forms.for'
C     forms.for
C
C=======================================================================
C