C 
C Program GREENSS to read a formatted file containing the ray-theory
C elastodynamic Green function and to generate ray-theory time-domain
C synthetic seismograms (without attenuation) or frequency-domain
C response functions (including causal Futterman's attenuation).
C
C Version: 5.10
C Date: 1997, September 26
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                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'GREEN','SOURCE','FRDAT','SIGNAL','SS',KOMP
C     'GREEN'... Name of the input formatted file with the Green tensor.
C             Description of file GREEN
C     'SOURCE'... Name of the input formatted file containing the
C             complex-valued seismic force or moment.
C             Description of file SOURCE
C     'FRDAT'... Name of the input formatted file containing the
C             specification of the frequency domain for the calculation
C             of response functions.
C             If blank, the time-domain synthetic seismograms are
C             generated instead of frequency-domain response functions,
C             see 'SIGNAL'.
C             Description of file FRDAT
C     'SIGNAL'... Name of the input file in the GSE format containing
C             the source time function and its Hilbert transform (for
C             seismic force), or their derivatives (for seismic moment).
C             To be submitted if 'SIGNAL'=' ', otherwise has no meaning.
C             (One of filenames 'FRDAT' and 'SIGNAL' must be nonblank.)
C             Description of GSE
C     'RF'... Name of the output file in the GSE format containing the
C             ray-theory seismograms or the frequency-domain Response
C             Function.
C             Description of file RF
C     KOMP... KOMP=0: All 3 components of the synthetic seismograms are
C               stored in the output GSE file.
C             KOMP=1: The 1st component of the synthetic seismograms is
C               stored in the output GSE file.
C             KOMP=2: The 2nd component of the synthetic seismograms is
C               stored in the output GSE file.
C             KOMP=3: The 3rd component of the synthetic seismograms is
C               stored in the output GSE file.
C     Default: 'GREEN'='green.out', 'SOURCE'='source.dat',
C             'FRDAT'='fr.dat', 'SIGNAL'='source.gse', 'RF'='rf.out'.
C             Note: File 'SIGNAL'='source.gse' need not exist (i.e. is
C               not read) if 'SOURCE' is not blank.
C
C                                                   
C Input formatted file GREEN:
C (1) / (a slash).
C (2) For each two-point ray (2.1):
C (2.1) 'R','S',(GREEN(I),I=1,32),/
C     'R'...  String in apostrophes identifying the receiver.  Only
C             the first 6 characters are written to the output GSE
C             file.  The strings corresponding to different receivers
C             thus should, if possible, differ in the first 6
C             characters.  If this is not the case, at least in the
C             first 12 characters.  All lines corresponding to the same
C             receiver must be consecutive.
C     'S'...  String in apostrophes describing the source.  Not taken
C             into account.
C     GREEN(1)... Travel time between receiver and source.
C     GREEN(2)... Imaginary part of the complex-valued travel time
C             between receiver and source due to attenuation.
C     GREEN(3:8)... Coordinates of the receiver and coordinates of the
C             source.
C     GREEN(9:14)... Derivatives of the travel time with respect to the
C             coordinates of the receiver and coordinates of the source.
C     GREEN(15:32)... 1000000 times enlarged amplitude of the Green
C             function: contravariant components of the complex-valued
C             3*3 matrix Gij in model coordinates, where the first
C             subscript corresponds to the receiver and the second
C             subscript corresponds to the source.  The components are
C             ordered as
C             Re(G11),Im(G11),Re(G21),Im(G21),Re(G31),Im(G31),
C             Re(G12),Im(G12),Re(G22),Im(G22),Re(G32),Im(G32),
C             Re(G13),Im(G13),Re(G23),Im(G23),Re(G33),Im(G33).
C     /...    An obligatory slash after at the end of line, in place
C             where the slowness vector components could be written.
C (3) / (a slash).
C
C                                                  
C Input formatted file SOURCE:
C (1) SFR1,SFI1,SFR2,SFI2,SFR3,SFI3,/
C     Components of the complex-valued vectorial seismic force.  The
C     seismic force is assumed to be the product of this vector and the
C     source time function submitted in file 'SIGNAL'.
C     Note:   The 'unit' radiation pattern corresponds to
C             SF = 4 pi rho v**2
C (2) SMR11,SMI11,SMR12,SMI12,SMR13,SMI13,
C     SMR21,SMI21,SMR22,SMI22,SMR23,SMI23,
C     SMR31,SMI31,SMR32,SMI32,SMR33,SMI33,/
C     Components of the transposed complex-valued 3*3 seismic-moment
C     tensor.  The tensor is transposed in order not to look transposed
C     in the input data.  The time derivative of the seismic moment is
C     assumed to be the product of this vector and the time function
C     submitted in file 'SIGNAL'.
C     Note:   The 'unit' radiation pattern corresponds to
C             SM = 4 pi rho v**3
C Example of data file SOURCE
C
C                                                   
C Input formatted file FRDAT:
C (1) NPTS,FMIN,FMAX,TD,TINT,TIMUL,FREF,/
C     NPTS... Number of time steps for the fast Fourier transform.
C             Will be used to convert the time step to the frequency
C             step.
C             NPTS=0: Frequency step is specified instead of time step.
C     FMIN,FMAX... Response functions are calculated for frequencies
C             from FMIN to FMAX.  FMIN and FMAX are rounded to the
C             nearest multiples of the time step.
C     TD...   NPTS=0: Frequency step.
C             Otherwise: Time step.  Frequency step=1/(NPTS*TD).
C     TINT... Maximum time interval.  The contribution of the ray to
C             the seismogram is taken into account only if the travel
C             time does not exceed TMIN+TINT, where TMIN is the minimum
C             travel time over the preceding rays arriving at the '
C             receiver.  Useful to remove alias in the time domain.
C             TINT=0: Infinite time interval.
C     TIMUL...Multiplication factor for the imaginary part TI of the
C             travel time.  It may be used to globally decrease or
C             increase attenuation in the whole model.
C     FREF... Reference frequency for the Futterman's (1962)
C             quasi-causal attenuation.  The travel times TR and TI
C             describing the Green function are assumed to correspond
C             to this frequency.  Frequency-dependent travel times are
C             then given by
C               Re TT(F)=TR-TI*ln(F/FREF)*2/pi,
C               Im TT(F)=TI.
C             FREF=0: Noncausal attenuation, just for test purposes,
C               Re TT(F)=TR,
C               Im TT(F)=TI.
C     Defaults: NPTS=0, FMIN=0, FMAX=1/(NPTS*TD) if NPTS.NE.0, TINT=0,
C             TIMUL=1, FREF=(FMIN+FMAX)/2.
C
C                                                      
C Output formatted file RF:
C If 'FRDAT'=' ':  Ray-theory time-domain synthetic seismograms (without
C     attenuation) in the GSE format.  The may be, e.g., plotted by
C     program 'sp.for'.
C     Program 'greenss.for' stores in the comment lines of the waveform
C     identification section the hypocentral coordinates identified by
C     strings 'XS1 ', 'XS2 ' and 'XS3 '.
C     Description of the GSE format
C Otherwise:  Ray-theory frequency-domain response functions (including
C     causal Futterman's attenuation), saved in the format 'RF'
C     described in 'ss.for'.
C     Description of file RF
C     Program 'greenss.for' is prepared to generate the response
C     functions also in the GSE format in future versions.  In such a
C     case program 'greenss.for' would store in the comment lines of the
C     waveform identification section the hypocentral coordinates
C     identified by strings 'XS1 ', 'XS2 ' and 'XS3 ', and times of the
C     first and last considered arrivals to each receiver, identified by
C     'TMIN' and 'TMAX'.
C     Description of the GSE format
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL WGSE1,WGSE2,WGSE3,RGSE2
C     WGSE1,WGSE2,WGSE3,RGSE2... File 'gse.for'.
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working arrays:
      INTEGER MSEIS
      PARAMETER (MSEIS=MRAM/6)
      REAL SEIS1(MSEIS),SEIS2(MSEIS),SEIS3(MSEIS)
      REAL SEIS4(MSEIS),SEIS5(MSEIS),SEIS6(MSEIS)
      EQUIVALENCE (SEIS1,RAM          )
      EQUIVALENCE (SEIS2,RAM(  MSEIS+1))
      EQUIVALENCE (SEIS3,RAM(2*MSEIS+1))
      EQUIVALENCE (SEIS4,RAM(3*MSEIS+1))
      EQUIVALENCE (SEIS5,RAM(4*MSEIS+1))
      EQUIVALENCE (SEIS6,RAM(5*MSEIS+1))
C
C-----------------------------------------------------------------------
C
C     Input and output files:
      INTEGER LU1,LU2
      PARAMETER (LU1=1,LU2=2)
      CHARACTER*80 FILE1,FILE2,FILE3,FILE4,FILE5,LINE
      CHARACTER*1  TEXT
C     Undefined value:
      REAL UNDEF
      PARAMETER (UNDEF=-999999.)
C     Seismic force:
      REAL SFR1,SFI1,SFR2,SFI2,SFR3,SFI3
C     Seismic moment:
      REAL SMR11,SMI11,SMR12,SMI12,SMR13,SMI13
      REAL SMR21,SMI21,SMR22,SMI22,SMR23,SMI23
      REAL SMR31,SMI31,SMR32,SMI32,SMR33,SMI33
C     Green function:
      CHARACTER*12 TXTOLD,TXTREC,TXTSRC
      REAL GREEN(32)
C     Data for the frequency band:
      INTEGER NPTS,NF
      REAL FMIN,FMAX,TD,TINT,TIMUL,FREF,FSTEP
C     Seismograms:
      LOGICAL LGSE,LWRITE
      INTEGER NSEIS,NSEIS4,NSEIS5
      REAL TSTAR0,TSTEP0,TSTARH,TSTART,TSTEP
C     Force at the source and displacement at the receiver:
      REAL FR1,FI1,FR2,FI2,FR3,FI3,AR1,AI1,AR2,AI2,AR3,AI3
C     Temporary storage locations:
      CHARACTER*1 STAT,CHAN
      INTEGER I01,I02,I03,I04,I05,I06,I07,I08,I09,I10,I11,I12
      INTEGER KOMP,ISHIFT,I
      REAL AMAX,AR,AI,TR,TI,F,OMEGA
      REAL XR1,XR2,XR3,XS1,XS2,XS3,TSHIFT,W0,W1
C     Source coordinates transferred through the GSE file:
      INTEGER MCOM,NCOM
      PARAMETER (MCOM=5)
      CHARACTER*4  TCOM(MCOM)
      REAL         VCOM(MCOM)
      DATA TCOM/'XS1 ','XS2 ','XS3 ','TMIN','TMAX'/
C
C.......................................................................
C
C     Format of the output response function (.TRUE.'GSE', .FALSE.'RF'):
      LGSE=.FALSE.
C
C     Names of input and output files:
      FILE1='green.out'
      FILE2='source.dat'
      FILE3='fr.dat'
      FILE4='source.gse'
      FILE5='rf.gse'
      KOMP=0
      WRITE(*,'(A)')
     *     ' Enter green.out,source.dat,fr.dat,source.gse,rf.gse,KOMP: '
      READ(*,*) FILE1,FILE2,FILE3,FILE4,FILE5,KOMP
C
C     Reading seismic force or seismic moment:
      OPEN(LU2,FILE=FILE2,STATUS='OLD')
      SFR1=0.
      SFI1=0.
      SFR2=0.
      SFI2=0.
      SFR3=0.
      SFI3=0.
      READ(LU2,*) SFR1,SFI1,SFR2,SFI2,SFR3,SFI3
      SMR11=0.
      SMI11=0.
      SMR12=0.
      SMI12=0.
      SMR22=0.
      SMI22=0.
      SMR13=0.
      SMI13=0.
      SMR23=0.
      SMI23=0.
      SMR33=0.
      SMI33=0.
      READ(LU2,*) SMR11,SMI11,SMR12,SMI12,SMR13,SMI13,
     *            SMR21,SMI21,SMR22,SMI22,SMR23,SMI23,
     *            SMR31,SMI31,SMR32,SMI32,SMR33,SMI33
      CLOSE(LU2)
C
C     Reading the data for the frequency domain:
      IF(FILE3.NE.' ') THEN
        OPEN(LU2,FILE=FILE3,STATUS='OLD')
        NPTS=0
        FMIN=0.
        FMAX=UNDEF
        TINT=0.
        TIMUL=1.
        FREF=UNDEF
        READ(LU2,*) NPTS,FMIN,FMAX,TD,TINT,TIMUL,FREF
        IF(FMAX.EQ.UNDEF) THEN
          FMAX=1/TD
        END IF
        IF(FREF.EQ.UNDEF) THEN
          FREF=(FMIN+FMAX)/2.
        END IF
        CLOSE(LU2)
        IF(NPTS.EQ.0) THEN
          FSTEP=TD
        ELSE
          FSTEP=1./(FLOAT(NPTS)*TD)
        END IF
        FMIN=FSTEP*INT(FMIN/FSTEP+.5)
        NF= INT((FMAX-FMIN)/FSTEP+.5)+1
C       Parameters for the response functions: FMIN,FSTEP,NF,TIMUL,FREF.
        NCOM=MCOM
      ELSE
        NCOM=3
C
C       Reading the source time function and its Hilbert transform:
        IF(FILE4.NE.' ') THEN
          OPEN(LU2,FILE=FILE4,STATUS='OLD')
          CALL RGSE2(LU2,STAT,CHAN,
     *                   I,XS1,XS2,XS3,TSTAR0,TSTEP0,NSEIS4,MSEIS,SEIS4)
          CALL RGSE2(LU2,STAT,CHAN,
     *                   I,XS1,XS2,XS3,TSTARH,TSTEP ,NSEIS5,MSEIS,SEIS5)
          CLOSE(LU2)
          IF(TSTEP0.NE.TSTEP) THEN
C           GREENSS-01
            PAUSE 'Error GREENSS-01: Different time steps.'
            STOP
          END IF
        ELSE
C         GREENSS-02
          PAUSE
     *    'Error GREENSS-02: One of files FRDAT or SIGNAL must be given'
          STOP
        END IF
C
      END IF
C
C.......................................................................
C
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
      OPEN(LU2,FILE=FILE5)
      IF(NCOM.NE.MCOM) THEN
        CALL WGSE1(LU2,' ')
      ELSE
        IF(LGSE) THEN
          CALL WGSE1(LU2,' ')
        ELSE
          LWRITE=.TRUE.
        END IF
      END IF
C
C     Loop over rays:
      TXTREC='$'
   30 CONTINUE
        TXTOLD=TXTREC
        TXTREC='$'
        XR1=GREEN(3)
        XR2=GREEN(4)
        XR3=GREEN(5)
        VCOM(1)=GREEN(6)
        VCOM(2)=GREEN(7)
        VCOM(3)=GREEN(8)
C       Preparing source coordinates for output in the GSE file:
        CALL WGSE2D()
        DO 31 I=1,NCOM
          CALL WSEPR(LINE,TCOM(I),VCOM(I))
          CALL WGSE2C(LINE)
   31   CONTINUE
C
        READ(LU1,*) TXTREC,TXTSRC,GREEN
        IF(TXTOLD.NE.'$'.AND.TXTREC.NE.TXTOLD) THEN
C         Writing the previous seismogram:
C         --------------------------------
          IF(NCOM.NE.MCOM) THEN
            IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN
              CALL WGSE2(LU2,TXTOLD,' ',1,XR1,XR2,XR3,
     *                             TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS1)
            END IF
            IF(KOMP.LE.0.OR.KOMP.EQ.2) THEN
              CALL WGSE2(LU2,TXTOLD,' ',2,XR1,XR2,XR3,
     *                             TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS2)
            END IF
            IF(KOMP.LE.0.OR.KOMP.EQ.3) THEN
              CALL WGSE2(LU2,TXTOLD,' ',3,XR1,XR2,XR3,
     *                             TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS3)
            END IF
            IF(NSEIS.GT.MSEIS) THEN
C             GREENSS-51
              WRITE(*,'(A,I12,2A)')
     *         ' Warning GREENSS-51:',NSEIS,
     *         ' non-zero samples at receiver ',TXTOLD
              PAUSE
            END IF
          ELSE
            IF(LGSE) THEN
              IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN
                CALL WGSE2(LU2,TXTOLD,' ', 1,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS1)
                CALL WGSE2(LU2,TXTOLD,' ',11,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS4)
              END IF
              IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN
                CALL WGSE2(LU2,TXTOLD,' ', 2,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS2)
                CALL WGSE2(LU2,TXTOLD,' ',12,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS5)
              END IF
              IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN
                CALL WGSE2(LU2,TXTOLD,' ', 3,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS3)
                CALL WGSE2(LU2,TXTOLD,' ',13,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS6)
              END IF
            ELSE
              IF(LWRITE)THEN
                WRITE(LU2,'(A)') '/'
                WRITE(LU2,'(3(F7.3,1X),A)') VCOM(1),VCOM(2),VCOM(3),'/'
                WRITE(LU2,'(2(E11.5,1X),I8,1X,A)') FMIN,FSTEP,NF,'/'
                LWRITE=.FALSE.
              END IF
              AMAX=0.
              DO 32 I=1,NF
                AMAX=AMAX1(ABS(SEIS1(I)),ABS(SEIS2(I)),ABS(SEIS3(I)),
     *                     ABS(SEIS4(I)),ABS(SEIS5(I)),ABS(SEIS6(I)),
     *                                                             AMAX)
   32         CONTINUE
              WRITE(LU2,'(3(F7.3,1X),3(E11.5,1X),A)')
     *                              XR1,XR2,XR3,VCOM(4),VCOM(5),AMAX,'/'
              IF(VCOM(4).LE.VCOM(5)) THEN
                DO 33 I=1,NF,2
                  I01=IFIX(99999.1*SEIS1(I)/AMAX)
                  I02=IFIX(99999.1*SEIS4(I)/AMAX)
                  I03=IFIX(99999.1*SEIS2(I)/AMAX)
                  I04=IFIX(99999.1*SEIS5(I)/AMAX)
                  I05=IFIX(99999.1*SEIS3(I)/AMAX)
                  I06=IFIX(99999.1*SEIS6(I)/AMAX)
                  IF(I.LT.NF) THEN
                    I07=IFIX(99999.1*SEIS1(I+1)/AMAX)
                    I08=IFIX(99999.1*SEIS4(I+1)/AMAX)
                    I09=IFIX(99999.1*SEIS2(I+1)/AMAX)
                    I10=IFIX(99999.1*SEIS5(I+1)/AMAX)
                    I11=IFIX(99999.1*SEIS3(I+1)/AMAX)
                    I12=IFIX(99999.1*SEIS6(I+1)/AMAX)
                    WRITE(LU2,'(12I6)') I01,I02,I03,I04,I05,I06,
     *                                  I07,I08,I09,I10,I11,I12
                  ELSE
                    WRITE(LU2,'(12I6)') I01,I02,I03,I04,I05,I06
                  END IF
   33           CONTINUE
              END IF
            END IF
          END IF
        END IF
        IF(TXTREC.EQ.'$') THEN
C         No more two-point rays:
C         -----------------------
          GO TO 80
        END IF
        IF(TXTREC.NE.TXTOLD) THEN
C         New receiver:
C         -------------
          DO 41 I=1,MSEIS
            SEIS1(I)=0.
            SEIS2(I)=0.
            SEIS3(I)=0.
   41     CONTINUE
          IF(NCOM.NE.MCOM) THEN
            ISHIFT=INT(GREEN(1)/TSTEP)
            TSTART=TSTAR0+FLOAT(ISHIFT)*TSTEP
            NSEIS=0
          ELSE
            VCOM(4)= 999999.
            VCOM(5)=-999999.
            DO 42 I=1,MSEIS
              SEIS4(I)=0.
              SEIS5(I)=0.
              SEIS6(I)=0.
   42       CONTINUE
          END IF
        END IF
        DO 43 I=15,32
          GREEN(I)=GREEN(I)/1000000.
   43   CONTINUE
C
C       Adding the contribution from a new two-point ray:
C       -------------------------------------------------
C
C       Complex-valued amplitude corresponding to the given source:
        FR1=SFR1-SMR11*GREEN(12)-SMR12*GREEN(13)-SMR13*GREEN(14)
        FI1=SFI1-SMI11*GREEN(12)-SMI12*GREEN(13)-SMI13*GREEN(14)
        FR2=SFR2-SMR21*GREEN(12)-SMR22*GREEN(13)-SMR23*GREEN(14)
        FI2=SFI2-SMI21*GREEN(12)-SMI22*GREEN(13)-SMI23*GREEN(14)
        FR3=SFR3-SMR31*GREEN(12)-SMR32*GREEN(13)-SMR33*GREEN(14)
        FI3=SFI3-SMI31*GREEN(12)-SMI32*GREEN(13)-SMI33*GREEN(14)
        AR1=GREEN(15)*FR1+GREEN(21)*FR2+GREEN(27)*FR3
     *     -GREEN(16)*FI1-GREEN(22)*FI2-GREEN(28)*FI3
        AR2=GREEN(17)*FR1+GREEN(23)*FR2+GREEN(29)*FR3
     *     -GREEN(18)*FI1-GREEN(24)*FI2-GREEN(30)*FI3
        AR3=GREEN(19)*FR1+GREEN(25)*FR2+GREEN(31)*FR3
     *     -GREEN(20)*FI1-GREEN(26)*FI2-GREEN(32)*FI3
        AI1=GREEN(15)*FI1+GREEN(21)*FI2+GREEN(27)*FI3
     *     +GREEN(16)*FR1+GREEN(22)*FR2+GREEN(28)*FR3
        AI2=GREEN(17)*FI1+GREEN(23)*FI2+GREEN(29)*FI3
     *     +GREEN(18)*FR1+GREEN(24)*FR2+GREEN(30)*FR3
        AI3=GREEN(19)*FI1+GREEN(25)*FI2+GREEN(31)*FI3
     *     +GREEN(20)*FR1+GREEN(26)*FR2+GREEN(32)*FR3
C
C       Time domain or frequency domain:
        IF(NCOM.NE.MCOM) THEN
C
C         Adding the multiple of the source time function:
          TSHIFT=(TSTAR0+GREEN(1)-TSTART)/TSTEP
          ISHIFT=INT(TSHIFT)
          W1=TSHIFT-FLOAT(ISHIFT)
          W0=1.-W1
C         W0,W1 are the weights of shifts ISHIFT,ISHIFT+1
          IF(ISHIFT.LT.0) THEN
C           Shifting the start time to the new position:
            TSTART=TSTART+FLOAT(ISHIFT)*TSTEP
            NSEIS=NSEIS-ISHIFT
            DO 51 I=MAX0(MSEIS,NSEIS),-ISHIFT+1,-1
              SEIS1(I)=SEIS1(I+ISHIFT)
              SEIS2(I)=SEIS2(I+ISHIFT)
              SEIS3(I)=SEIS3(I+ISHIFT)
   51       CONTINUE
            DO 52 I=MAX0(MSEIS,-ISHIFT),1,-1
              SEIS1(I)=0.
              SEIS2(I)=0.
              SEIS3(I)=0.
   52       CONTINUE
            ISHIFT=0
          END IF
          NSEIS=MAX0(NSEIS,NSEIS4+ISHIFT)
          DO 53 I=1,MIN(NSEIS4,MSEIS-ISHIFT)
            SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)+W0*AR1*SEIS4(I)
            SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)+W0*AR2*SEIS4(I)
            SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)+W0*AR3*SEIS4(I)
   53     CONTINUE
          ISHIFT=ISHIFT+1
          NSEIS=MAX0(NSEIS,NSEIS4+ISHIFT)
          DO 54 I=1,MIN(NSEIS4,MSEIS-ISHIFT)
            SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)+W1*AR1*SEIS4(I)
            SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)+W1*AR2*SEIS4(I)
            SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)+W1*AR3*SEIS4(I)
   54     CONTINUE
C
C         Adding the multiple of the Hilbert transform:
          TSHIFT=(TSTARH+GREEN(1)-TSTART)/TSTEP
          ISHIFT=INT(TSHIFT)
          W1=TSHIFT-FLOAT(ISHIFT)
          W0=1.-W1
C         W0,W1 are the weights of shifts ISHIFT,ISHIFT+1
          IF(ISHIFT.LT.0) THEN
C           Shifting the start time to the new position:
            TSTART=TSTART+FLOAT(ISHIFT)*TSTEP
            NSEIS=NSEIS-ISHIFT
            DO 61 I=MIN0(MSEIS,NSEIS),-ISHIFT+1,-1
              SEIS1(I)=SEIS1(I+ISHIFT)
              SEIS2(I)=SEIS2(I+ISHIFT)
              SEIS3(I)=SEIS3(I+ISHIFT)
   61       CONTINUE
            DO 62 I=MIN0(MSEIS,-ISHIFT),1,-1
              SEIS1(I)=0.
              SEIS2(I)=0.
              SEIS3(I)=0.
   62       CONTINUE
            ISHIFT=0
          END IF
          NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT)
          DO 63 I=1,MIN(NSEIS5,MSEIS-ISHIFT)
            SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)-W0*AI1*SEIS5(I)
            SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)-W0*AI2*SEIS5(I)
            SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)-W0*AI3*SEIS5(I)
   63     CONTINUE
          ISHIFT=ISHIFT+1
          NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT)
          DO 64 I=1,MIN(NSEIS5,MSEIS-ISHIFT)
            SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)-W1*AI1*SEIS5(I)
            SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)-W1*AI2*SEIS5(I)
            SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)-W1*AI3*SEIS5(I)
   64     CONTINUE
        ELSE
C
C         Response functions:
          VCOM(4)=AMIN1(GREEN(1),VCOM(4))
          IF(TINT.EQ.0..OR.GREEN(1).LE.VCOM(4)+TINT) THEN
            VCOM(5)=AMAX1(GREEN(1),VCOM(5))
            DO 69 I=1,NF
              F=FMIN+FSTEP*FLOAT(I-1)
              OMEGA=6.2831853*F
              TI=GREEN(2)*TIMUL
              IF(FREF.EQ.0.) THEN
                TR=GREEN(1)
              ELSE
                TR=GREEN(1)-TI*ALOG(F/FREF)/1.5707963
              END IF
              TI=  EXP(-OMEGA*TI)
              AR=TI*COS(OMEGA*TR)
              AI=TI*SIN(OMEGA*TR)
              SEIS1(I)=SEIS1(I)+AR1*AR-AI1*AI
              SEIS2(I)=SEIS2(I)+AR2*AR-AI2*AI
              SEIS3(I)=SEIS3(I)+AR3*AR-AI3*AI
              SEIS4(I)=SEIS4(I)+AI1*AR+AR1*AI
              SEIS5(I)=SEIS5(I)+AI2*AR+AR2*AI
              SEIS6(I)=SEIS6(I)+AI3*AR+AR3*AI
   69       CONTINUE
          END IF
        END IF
      GO TO 30
C
   80 CONTINUE
      IF(NCOM.NE.MCOM) THEN
        CALL WGSE3(LU2)
      ELSE
        IF(LGSE)THEN
          CALL WGSE3(LU2)
        ELSE
          WRITE(LU2,'(A)') '/'
        END IF
      END IF
      CLOSE(LU2)
      CLOSE(LU1)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'gse.for'
C     gse.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C