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 Program GREENSS can also read frequency-dependent coupling-ray-theory
C Green functions corresponding to S waves in weakly anisotropic models.
C
C Version GREENRM of program GREENSS may include the response of the
C transmission through a stack of fine horizontal layers at each
C receiver location.
C The frequency-dependent transmission coefficiens are calculated by
C subroutines of package RMATRIX written by C. J. Thomson (Copyright (c)
C 1996 C. J. Thomson).  Please, refer to package
C RMATRIX for copyright conditions
C and details on compilation.
C Version GREENRM is derived from the basic version of program GREENSS
C by changing the asterisks '*' in the first column of the INCLUDE
C statements at the end of this file.
C
C Version: 5.20
C Date: 1998, November 11
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) 'SEP',/
C     'SEP'...Name of the input formatted SEP file containing the
C             specification of the frequency domain for the calculation
C             of response functions.
C             This file is not read if 'SIGGSE'=' '.
C             Description of file SEP
C     Default: 'SEP'='ss.h'.
C
C                                                   
C Input file 'SEP' in the SEP format:
C     The file has the form of a SEP parameter file.
C     For the description of the SEP format refer to file
C     'sep.for'.
C Names of input files:
C     GREEN='string'... Name of the input formatted file with the Green
C             tensor.  The file is generated by program 'green.for'.
C             Description of file GREEN
C             Default: GREEN='green.out' (mostly sufficient)
C     SOURCE='string'... Name of the input formatted file containing the
C             complex-valued seismic force or moment.
C             Description of file SOURCE
C             Default: SOURCE='source.dat'
C     RMATRIX=string... Name of optional formatted input file with the
C             parameters of fine layered structures at the receiver
C             sites.  If non-blank, the frequency-dependent transmission
C             coefficiens, calculated by subroutines of package RMATRIX,
C             are considered at each receiver.
C             Available just for version GREENRM of program GREENSS, see
C             the comment at the beginning of this file.
C             Default: RMATRIX=' '
C     SIGGSE='string'... Name of the optional input file in the GSE
C             format containing the source time function and its Hilbert
C             transform (for seismic force), or their derivatives (for
C             seismic moment).
C             If blank (default), the frequency-domain response
C             functions are generated.
C             If non-blank, the time-domain synthetic seismograms are
C             generated instead of frequency-domain response functions.
C             It is recommended not to specify SIGGSE (use default).
C             Description of GSE
C             Default: SIGGSE=' ' (recommended)
C Names of output files:
C     RF='string'... Name of the output file containing the
C             frequency-domain Response Function.
C             Generated just if SIGGSE=' ' (default).
C             Description of file RF
C             Default: RF='rf.out' (mostly sufficient)
C     SS='string'... Name of the optional output file in the GSE format
C             containing the ray-theory seismograms calculated in the
C             time domain.  Generated just if SIGGSE.NE.' '.
C             Description of file SS
C             Default: SS='ss.gse'
C Data describing the frequency domain common with program 'ss.for':
C     DT=real... Time step to digitize the source time function and
C             seismograms.
C             Default: DT=1.
C     NFFT=integer... Number of the time samples for the fast Fourier
C             transform.  Will be used to convert the time step to the
C             frequency step.
C             NFFT must be an integer power of 2.
C             Default: NFFT=1
C     FMIN=real... The lowest frequency of the cosine band-pass filter.
C             Default: FMIN=0.
C     FMAX=real... The highest frequency of the cosine band-pass filter.
C             Default: FMAX=0.5/DT (Nyquist frequency)
C Data describing the frequency domain specific to this program:
C     DF...   Frequency step.
C             Default: DF=1/(NPTS*DT) (recommended)
C     OF...   Response functions are calculated for NF frequencies
C             OF,OF+DF,OF+2*DF,...,OF+(NF-1)*DF.
C             Default: OF=DF*NINT(FMIN/DF) (i.e. FMIN rounded to the
C               nearest multiple of the frequency step, recommended)
C     NF...   Number of frequencies.
C             Default: NF=NINT((FMAX-OF)/DF)+1 (recommended)
C     Hint: Do not specify DF, OF and NF.  Use the data for program
C     'ss.for' instead.
C Data specific to this program:
C     TINT=non-negative real... Maximum time interval.  The contribution
C             of the ray to the seismogram is taken into account only if
C             the travel time does not exceed TMIN+TINT, where TMIN is
C             the minimum travel time over the preceding rays arriving
C             at the receiver.  Useful to remove alias in the time
C             domain.
C             TINT=0: Infinite time interval.
C             Default: TINT=0. (mostly sufficient for reasonable NPTS)
C     TIMUL=non-negative real... Multiplication factor for the imaginary
C             part TI of the travel time.  It may be used to globally
C             decrease or increase attenuation in the whole model.
C             Default: TIMUL=1. (mostly sufficient)
C     FREF=non-negative real... Reference frequency for the Futterman's
C             (1962) quasi-causal attenuation.  The travel times TR and
C             TI 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             Default: FREF=(FMIN+FMAX)/2.
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 'SIGGSE'.
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 'SIGGSE'.
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 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     GREEN(33:*)... Analogy of GREEN(15:32) for other frequencies.
C             Present just for frequency-dependent elementary Green
C             functions.
C     /...    A slash after at the end of line.  It may also stand for
C             zeros at the end of line.
C (3) / (a slash).
C
C                                                 
C Optional input formatted file RMATRIX:
C     Parameters of fine layered structures at the receiver sites.
C Once or for each receiver the following data:
C (1) 'REC',/
C     'REC'...Blank string if a single fine layered structure applies to
C             all receivers.  Otherwise the name of the receiver.
C             The order and names of the receivers must be the same as
C             in file 'GREEN'.  The order is given by the receiver file.
C             If the receiver file is specified for the GREEN program,
C             it also determines the receiver names.
C     Default: 'REC'=' '.
C (2) NLAY,IFREE
C     NLAY... Number of layers, including the bottom and top halfspaces.
C     IFREE...IFREE=1: The thickness THICKN(1) of the first (top) layer
C               is used for the calculation of the travel time
C               correction at the geophone situated in the top halfspace
C               at the elevation THICKN(1) above the uppermost
C               interface.
C             IFREE=0: Has the same effect as THICKN(1)=0.
C (3) ISOFLG,THICKN,RHO
C     ISOFLG..ISOFLG=0: Anisotropic layer.
C             ISOFLG=1: Isotropic elastic layer.
C             ISOFLG=2: Fluid layer.  Does not work in the current
C               version.  ISOFLG=1 with a very small S wave velocity is
C               suggested instead.
C     THICKN..Thickness of the layer.  THICKN(1) and THICKN(NLAY) are
C             used to position the layer with respect to the coordinates
C             of the receiver.  The thicknesses may also be negative.
C             The lowermost interface (top of the bottom halfspace) is
C             situated at the elevation THICKN(NLAY) above the receiver
C             coordinates.  The vertical component of the slowness
C             vector is used to linearly extrapolate travel time from
C             the receiver coordinates to the bottom of the layer stack.
C             The geophone is considered to be situated in the elevation
C             THICKN(1)+THICKN(2)+...+THICKN(NLAY) above the receiver.
C             Example:
C               For a receiver at the Earth surface, layer 1 has
C               THICK(1)=0 (geophone exactly at the surface),  RHO(1)
C               equal to the density of the air (about 0.0012kg/l) or
C               smaller, VP(1) equal to the sound speed in the air
C               (about 0.34km/s), and VS(1) very small.  Layers 2 to
C               NLAY-1 correspond to the properties of the fine layered
C               structure.  Layer NLAY usually takes negative thickness
C               THICKN(NLAY)=-THICKN(1)-THICKN(2)-...-THICKN(NLAY)
C               in order to situate geophone exactly at the receiver
C               coordinates.
C     RHO...  Density of the layer.
C (4) For ISOFLG=0 (anisotropic layer) (4.1), otherwise (4.2):
C (4.1) (A(I,J,K,L),I=1,3,J=1,3,K=1,3,L=1,3)
C     A(I,J,K,L)... Matrix of 81 elastic parameters C(I,J,K,L) divided
C             by density RHO.
C (4.2) VP,VS
C     VP...   P wave velocity.
C     VS...   S wave velocity.  May be very small for fluids but cannot
C             be zero in this version.
C
C                                                      
C Output formatted file RF:
C     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 Optional output formatted file SS:
C     Ray-theory time-domain synthetic seismograms (without
C     attenuation) in the GSE format.  They 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
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 MGREEN,MSEIS
      PARAMETER (MGREEN=MRAM/4*3)
      PARAMETER (MSEIS=(MRAM-MGREEN)/6)
      REAL SEIS1(MSEIS),SEIS2(MSEIS),SEIS3(MSEIS)
      REAL SEIS4(MSEIS),SEIS5(MSEIS),SEIS6(MSEIS)
      REAL GREEN(MGREEN)
      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))
      EQUIVALENCE (GREEN,RAM(6*MSEIS+1))
C     Frequency-dependent Green function requires 14+18*NF storage
C     locations.
C
C-----------------------------------------------------------------------
C
C     Input and output files:
      INTEGER LU1,LU2,LU3
      PARAMETER (LU1=1,LU2=2,LU3=3)
      CHARACTER*80 FILE1,FILE2,FILSEP,FILSIG,FILOUT,FILRM,LINE
      CHARACTER*80 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     Strings to identify elementary Green functions:
      CHARACTER*12 TXTOLD,TXTREC,TXTSRC,TXTRM
C     Green function:
      INTEGER NGREEN
C     Input and output arguments for RMATRIX:
      INTEGER NLAY,IFREE
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
      REAL AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3
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,J
      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     Main input file (SEP format):
      FILSEP='ss.h'
      WRITE(*,'(A)') '+GREENSS: Enter input filename: '
      READ(*,*) FILSEP
      WRITE(*,'(A)') '+GREENSS: Working...            '
      CALL RSEP1(LU1,FILSEP)
C
C     Reading seismic force or seismic moment:
      CALL RSEP3T('SOURCE',FILE2 ,'source.dat')
      OPEN(LU1,FILE=FILE2,STATUS='OLD')
      SFR1=0.
      SFI1=0.
      SFR2=0.
      SFI2=0.
      SFR3=0.
      SFI3=0.
      READ(LU1,*) 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(LU1,*) SMR11,SMI11,SMR12,SMI12,SMR13,SMI13,
     *            SMR21,SMI21,SMR22,SMI22,SMR23,SMI23,
     *            SMR31,SMI31,SMR32,SMI32,SMR33,SMI33
      CLOSE(LU1)
C
      CALL RSEP3T('SIGGSE',FILSIG,' ')
      IF(FILSIG.EQ.' ') THEN
C
C       Reading the data for the frequency domain:
        CALL RSEP3I('NFFT ',NPTS  ,1)
        CALL RSEP3R('DT   ',TD    ,1.)
        CALL RSEP3R('FMIN ',FMIN  ,0.)
        CALL RSEP3R('FMAX ',FMAX  ,0.5/TD)
        CALL RSEP3R('TINT ',TINT  ,0.)
        CALL RSEP3R('TIMUL',TIMUL ,1.)
        CALL RSEP3R('FREF ',FREF  ,(FMIN+FMAX)/2.)
        CALL RSEP3R('DF   ',FSTEP ,1./(FLOAT(NPTS)*TD))
        CALL RSEP3R('OF   ',FMIN  ,FSTEP*NINT(FMIN/FSTEP))
        CALL RSEP3I('NF   ',NF    ,NINT((FMAX-FMIN)/FSTEP)+1)
        IF(14+18*NF.GT.MGREEN) THEN
C         GREENSS-01
          CALL ERROR('GREENSS-01: Small array for Green functions')
C         Frequency-dependent Green function requires 14+18*NF storage
C         locations.  Array GREEN located within array RAM declared in
C         include file ram.inc has
C         MGREEN=MRAM/4*3 locations.
        END IF
        CALL RSEP3T('RMATRIX',FILRM,' ')
C       Parameters for the response functions: FMIN,FSTEP,NF,TIMUL,FREF.
        NCOM=MCOM
        CALL RSEP3T('RF',FILOUT,'rf.out')
C
      ELSE
C
C       Reading the source time function and its Hilbert transform:
        OPEN(LU1,FILE=FILSIG,STATUS='OLD')
        CALL RGSE2(LU1,STAT,CHAN,
     *                 I,XS1,XS2,XS3,TSTAR0,TSTEP0,NSEIS4,MSEIS,SEIS4)
        CALL RGSE2(LU1,STAT,CHAN,
     *                 I,XS1,XS2,XS3,TSTARH,TSTEP ,NSEIS5,MSEIS,SEIS5)
        CLOSE(LU1)
        IF(TSTEP0.NE.TSTEP) THEN
C         GREENSS-02
          CALL ERROR('GREENSS-02: Different time steps.')
        END IF
        FILRM=' '
        IF(FILRM.NE.' ') THEN
C         GREENSS-04
          CALL ERROR('GREENSS-04: Layer stack in time domain')
C         The response of the transmission through a stack of fine
C         horizontal layers at the receiver locations cannot be
C         calculated in the time domain.
        END IF
        NCOM=3
        CALL RSEP3T('SS',FILOUT,'ss.gse')
C
      END IF
C
C.......................................................................
C
C     Opening input file with the Green function:
      CALL RSEP3T('GREEN' ,FILE1 ,'green.out' )
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
C
C     Opening the output file:
      OPEN(LU3,FILE=FILOUT)
      IF(NCOM.NE.MCOM) THEN
C       Time domain
        CALL WGSE1(LU3,' ')
      ELSE
C       Frequency domain
        IF(LGSE) THEN
          CALL WGSE1(LU3,' ')
        ELSE
          LWRITE=.TRUE.
        END IF
      END IF
C
C     Removed from the input data in version 5.20:
C     CALL RSEP3I('KOMP',KOMP,0)
      KOMP=0
C     KOMP=integer... The default value is recommended.
C             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: KOMP=0 (recommended)
C
C     Loop over rays:
      TXTREC='$'
      TXTRM='$'
   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
C       Reading one elementary Green function:
        IF(NCOM.NE.MCOM) THEN
C         Time domain
          NGREEN=33
        ELSE
C         Frequency domain
          NGREEN=14+18*NF
        END IF
        DO 32 I=1,NGREEN
          GREEN(I)=0.
   32   CONTINUE
        GREEN(33)=UNDEF
        READ(LU1,*) TXTREC,TXTSRC,(GREEN(I),I=1,NGREEN)
        IF(GREEN(33).EQ.UNDEF) THEN
C         Frequency-independent Green function:
          NGREEN=32
        ELSE
C         Frequency-dependent Green function:
          IF(NCOM.NE.MCOM) THEN
C           Time domain calculation
C           GREENSS-05
            CALL ERROR('GREENSS-05: Frequency dependent Green function')
C           Frequency-dependent Green function has been identified in
C           the input file during the calculation of the seismograms in
C           the time domain.
          END IF
          NGREEN=14+18*NF
        END IF
        IF(TXTOLD.NE.'$'.AND.TXTREC.NE.TXTOLD) THEN
C         Writing the previous seismogram:
C         --------------------------------
          IF(NCOM.NE.MCOM) THEN
C           Time domain
            IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN
              CALL WGSE2(LU3,TXTOLD,' ',1,XR1,XR2,XR3,
     *                             TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS1)
            END IF
            IF(KOMP.LE.0.OR.KOMP.EQ.2) THEN
              CALL WGSE2(LU3,TXTOLD,' ',2,XR1,XR2,XR3,
     *                             TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS2)
            END IF
            IF(KOMP.LE.0.OR.KOMP.EQ.3) THEN
              CALL WGSE2(LU3,TXTOLD,' ',3,XR1,XR2,XR3,
     *                             TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS3)
            END IF
            IF(NSEIS.GT.MSEIS) THEN
C             GREENSS-51
              WRITE(TEXT,'(A,I12,2A)') 'GREENSS-51:',NSEIS,
     *                           ' non-zero samples at receiver ',TXTOLD
              CALL WARN(TEXT)
            END IF
          ELSE
C           Frequency domain
            IF(LGSE) THEN
              IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN
                CALL WGSE2(LU3,TXTOLD,' ', 1,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS1)
                CALL WGSE2(LU3,TXTOLD,' ',11,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS4)
              END IF
              IF(KOMP.LE.0.OR.KOMP.EQ.2) THEN
                CALL WGSE2(LU3,TXTOLD,' ', 2,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS2)
                CALL WGSE2(LU3,TXTOLD,' ',12,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS5)
              END IF
              IF(KOMP.LE.0.OR.KOMP.EQ.3) THEN
                CALL WGSE2(LU3,TXTOLD,' ', 3,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS3)
                CALL WGSE2(LU3,TXTOLD,' ',13,XR1,XR2,XR3,
     *                                              FMIN,FSTEP,NF,SEIS6)
              END IF
            ELSE
              IF(LWRITE)THEN
                WRITE(LU3,'(A)') '/'
                WRITE(LU3,'(3(F7.3,1X),A)') VCOM(1),VCOM(2),VCOM(3),'/'
                WRITE(LU3,'(2(E11.5,1X),I8,1X,A)') FMIN,FSTEP,NF,'/'
                LWRITE=.FALSE.
              END IF
              AMAX=0.
              DO 33 I=1,NF
                AMAX=AMAX1(ABS(SEIS1(I)),ABS(SEIS2(I)),ABS(SEIS3(I)),
     *                     ABS(SEIS4(I)),ABS(SEIS5(I)),ABS(SEIS6(I)),
     *                                                             AMAX)
   33         CONTINUE
              WRITE(LU3,'(3A,3(F7.3,1X),3(E11.5,1X),A)')
     *            '''',TXTOLD,''' ',XR1,XR2,XR3,VCOM(4),VCOM(5),AMAX,'/'
              IF(VCOM(4).LE.VCOM(5)) THEN
                DO 34 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(LU3,'(12I6)') I01,I02,I03,I04,I05,I06,
     *                                  I07,I08,I09,I10,I11,I12
                  ELSE
                    WRITE(LU3,'(12I6)') I01,I02,I03,I04,I05,I06
                  END IF
   34           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         -------------
          IF(NCOM.NE.MCOM) THEN
C           Time domain
            ISHIFT=INT(GREEN(1)/TSTEP)
            TSTART=TSTAR0+FLOAT(ISHIFT)*TSTEP
            NSEIS=0
          ELSE
C           Frequency domain
            VCOM(4)= 999999.
            VCOM(5)=-999999.
            DO 41 I=1,NF
              SEIS1(I)=0.
              SEIS2(I)=0.
              SEIS3(I)=0.
              SEIS4(I)=0.
              SEIS5(I)=0.
              SEIS6(I)=0.
   41       CONTINUE
          END IF
C
C         Reading fine layered structure at the receiver location
          IF(FILRM.NE.' ') THEN
            IF(TXTRM.EQ.'$') THEN
C             Opening file RMATRIX with fine layered structures at the
C             receiver sites:
              OPEN(11 ,STATUS='SCRATCH')
              OPEN(LU2,FILE=FILRM,STATUS='OLD')
            END IF
            IF(TXTRM.NE.' ') THEN
C             Searching for the receiver in file RMATRIX
   43         CONTINUE
                TXTRM=' '
                READ(LU2,*,ERR=44,END=44) TXTRM
                CALL INPUT(LU2,NLAY,IFREE)
              IF(TXTRM.NE.' '.AND.TXTRM.NE.TXTREC) GO TO 43
              GO TO 45
C             ..........................................................
C             The receiver has not been found
   44         CONTINUE
C             GREENSS-06
              CALL ERROR('GREENSS-06: Receiver is not in file RMATRIX')
C             The receiver has not been found in file
C             RMATRIX
C             ..........................................................
   45         CONTINUE
            END IF
          END IF
        END IF
C
        DO 49 I=15,NGREEN
          GREEN(I)=GREEN(I)/1000000.
   49   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         Time domain
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=MIN0(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=MIN0(MSEIS,-ISHIFT),1,-1
              SEIS1(I)=0.
              SEIS2(I)=0.
              SEIS3(I)=0.
   52       CONTINUE
            ISHIFT=0
          END IF
          DO 53 I=NSEIS+1,MIN0(NSEIS4+ISHIFT+1,MSEIS)
            SEIS1(I)=0.
            SEIS2(I)=0.
            SEIS3(I)=0.
   53     CONTINUE
          DO 54 I=1,MIN0(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)
   54     CONTINUE
          ISHIFT=ISHIFT+1
          NSEIS=MAX0(NSEIS,NSEIS4+ISHIFT)
          DO 55 I=1,MIN0(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)
   55     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
          DO 63 I=NSEIS+1,MIN0(NSEIS5+ISHIFT+1,MSEIS)
            SEIS1(I)=0.
            SEIS2(I)=0.
            SEIS3(I)=0.
   63     CONTINUE
          NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT)
          DO 64 I=1,MIN0(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)
   64     CONTINUE
          ISHIFT=ISHIFT+1
          NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT)
          DO 65 I=1,MIN0(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)
   65     CONTINUE
        ELSE
C         Frequency domain
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.LE.0..OR.TI.LE.0.) THEN
                TR=GREEN(1)
              ELSE
                TR=GREEN(1)-TI*ALOG(F/FREF)/1.5707963
              END IF
C
C             Frequency-dependent green function
              IF(NGREEN.GT.32) THEN
                J=18*(I-1)
                AR1=GREEN(J+15)*FR1+GREEN(J+21)*FR2+GREEN(J+27)*FR3
     *             -GREEN(J+16)*FI1-GREEN(J+22)*FI2-GREEN(J+28)*FI3
                AR2=GREEN(J+17)*FR1+GREEN(J+23)*FR2+GREEN(J+29)*FR3
     *             -GREEN(J+18)*FI1-GREEN(J+24)*FI2-GREEN(J+30)*FI3
                AR3=GREEN(J+19)*FR1+GREEN(J+25)*FR2+GREEN(J+31)*FR3
     *             -GREEN(J+20)*FI1-GREEN(J+26)*FI2-GREEN(J+32)*FI3
                AI1=GREEN(J+15)*FI1+GREEN(J+21)*FI2+GREEN(J+27)*FI3
     *             +GREEN(J+16)*FR1+GREEN(J+22)*FR2+GREEN(J+28)*FR3
                AI2=GREEN(J+17)*FI1+GREEN(J+23)*FI2+GREEN(J+29)*FI3
     *             +GREEN(J+18)*FR1+GREEN(J+24)*FR2+GREEN(J+30)*FR3
                AI3=GREEN(J+19)*FI1+GREEN(J+25)*FI2+GREEN(J+31)*FI3
     *             +GREEN(J+20)*FR1+GREEN(J+26)*FR2+GREEN(J+32)*FR3
              END IF
C
C             Adding the response of the fine layered structure at the
C             receiver location to the complex-valued amplitude:
              IF(FILRM.EQ.' ') THEN
                BR1=AR1
                BI1=AI1
                BR2=AR2
                BI2=AI2
                BR3=AR3
                BI3=AI3
              ELSE
                CALL RM(NLAY,IFREE,I,
     *                  OMEGA,TR,GREEN(9),GREEN(10),GREEN(11),
     *                  AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3)
              END IF
C
C             Complex-valued phase factor (AR,AI)
              TI=  EXP(-OMEGA*TI)
              AR=TI*COS(OMEGA*TR)
              AI=TI*SIN(OMEGA*TR)
C
C             Contribution of the ray to the response function:
              SEIS1(I)=SEIS1(I)+BR1*AR-BI1*AI
              SEIS2(I)=SEIS2(I)+BR2*AR-BI2*AI
              SEIS3(I)=SEIS3(I)+BR3*AR-BI3*AI
              SEIS4(I)=SEIS4(I)+BI1*AR+BR1*AI
              SEIS5(I)=SEIS5(I)+BI2*AR+BR2*AI
              SEIS6(I)=SEIS6(I)+BI3*AR+BR3*AI
   69       CONTINUE
          END IF
        END IF
      GO TO 30
C
   80 CONTINUE
      IF(NCOM.NE.MCOM) THEN
C       Time domain
        CALL WGSE3(LU3)
      ELSE
C       Frequency domain
        IF(LGSE)THEN
          CALL WGSE3(LU3)
        ELSE
          WRITE(LU3,'(A)') '/'
        END IF
      END IF
      CLOSE(LU3)
      IF(FILRM.NE.' ') THEN
        CLOSE(LU2)
        CLOSE(11)
      END IF
      CLOSE(LU1)
      WRITE(*,'(A)') '+GREENSS: 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 Response of fine layers:
C Include just one set of the following files:
C (a) Version GREENSS: No response of a stack of fine horizontal layers:
      INCLUDE 'rmnul.for'
C     rmnul.for
C (b) Version GREENRM: Stacks of fine horizontal layers at receivers:
*     INCLUDE 'rm.for'
C     rm.for
*     INCLUDE 'rmatrix1.f'
C     rmatrix1.f
*     INCLUDE 'rmatrix2.f'
C     rmatrix2.f
*     INCLUDE 'isotroc.f'
C     isotroc.f
*     INCLUDE 'aniso6cg.f'
C     aniso6cg.f
*     INCLUDE 'aniso6c1.f'
C     aniso6c1.f
*     INCLUDE 'aniso3c.f'
C     aniso3c.f
*     INCLUDE 'deigv.f'
C     deigv.f
*     INCLUDE 'cg.f'
C     cg.f
C
C=======================================================================
C