C
C Program SP (Seismogram Plotting) to plot seismograms previously stored
C in the GSE data exchange format.
C
C Version: 5.40
C Date: 2000, January 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     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of the input files:
C     SS='string'... String with the name of the input data file in the
C             GSE data exchange format, containing the seismograms to be
C             plotted.
C             Description of file GSE
C             Default: SS='ss.gse'
C     SS1='string', SS2='string', ..., SS9='string'... Strings with the
C             additional, optional names of the input data files in the
C             GSE data exchange format, containing the seismograms to be
C             plotted simultaneously with seismograms given by parameter
C             SS into the same figures.  The order of plottting is SS,
C             SS1, SS2, ..., SS9, considering just nonblank filenames.
C             The seismograms are plotted in colours given by parameters
C             KOLOR, KOLOR1, KOLOR2, ..., KOLOR9, respectively.  Refer
C             to file calcops.rgb for the
C             definitions of the colours.  The frame and description
C             remain in colour number 1, which is probably black.
C             Defaults: SS1=' ', SS2=' ', SS3=' ', ..., SS9=' '
C     SRC='string'... String with the name of the input data file
C             with the coordinates of the hypocentre.  If left blank
C             (default), the hypocentral coordinates are deemed zero.
C             The hypocentral coordinates are used only for eventual
C             travel-time reduction on a given velocity, or for
C             amplitude power scaling.  Otherwise, this filename need
C             not be specified.  This file also need not be specified if
C             the seismograms are generated by programs 'greenss.for' or
C             'ss.for', because those programs enter the hypocentral
C             coordinates directly to the comments of the data section
C             in GSE data exchange file.  The hypocentral coordinates
C             from the GSE file take precedence of the ones from the
C             'SRC' file.
C             Description of file SRC
C             Default: 'SRC'=' '
C     REC='string'... String with the name of the input data file
C             with the list of the receivers.
C             If REC=' ' seismograms at all receivers will be plotted,
C             otherwise only seismograms at the receivers listed in the
C             file will be plotted.
C             This file also enables to determine the number NREC of all
C             receivers, in which the seimograms may be plotted, for
C             scaling purposes.  If REC=' ', NREC=999 for scaling.
C             In most cases, file REC will be the same as for the
C             calculation of seismograms.
C             Description of file REC
C             Default: 'REC'=' '
C     SPHILI='string'... String with the name of the input data file
C             with the list of the travel times to be highlighted.
C             If SPHILI=' ', nothing is highlighted.
C             Description of file SPHILI
C             Default: 'SPHILI'=' '
C Names of output files:
C     SP1='string'... String with the name of the output PostScript file
C             with the plot of the first component of the synthetic
C             seismograms.
C             Default: SP1='ss1.gse'
C     SP2='string'... String with the name of the output PostScript file
C             with the plot of the second component of the synthetic
C             seismograms.
C             Default: SP2='ss2.gse'
C     SP3='string'... String with the name of the output PostScript file
C             with the plot of the third component of the synthetic
C             seismograms.
C             Default: SP3='ss3.gse'
C Data to control plotting:
C Colours:
C     KOLOR=positive integer, KOLOR1=positive integer,
C     KOLOR2=positive integer, ..., KOLOR9=positive integer... Colours
C             to plot seismograms of files SS,SS1, SS2, ..., SS9,
C             respectively.  Colour indices correspond to dummy
C             parameter INP of CalComp subroutine
C             NEWPEN.
C             The colours corresponding to the indices may be defined
C             or changed in file calcops.rgb.
C             Default: KOLOR=1, KOLOR1=2, KOLOR2=3, ..., KOLOR9=10
C Data for the time axis (vertical):
C     SPTMIN=real... Minimum time (or reduced time), i.e.
C             the time corresponding to the bottom of the
C             seismogram plot.
C             Default: SPTMIN=0.
C     SPTMAX=real... Maximum time (or reduced time), i.e.
C             the time corresponding to the top of the
C             seismogram plot.
C             Default: SPTMAX=1.
C     SPTLEN=real... Length of the vertical time axis in cm.
C             Default: SPTLEN=10.
C     SPTDIV=real... ABS(SPTDIV) is the number of intervals along the
C             time axis, starting at the bottom.  Must be SPTDIV.NE.0.
C             SPTDIV.GT.0.: the marks at the endpoints of intervals will
C               be supplemented with corresponding times (or reduced
C               times).
C             SPTDIV.LT.0.: the time axis will have no description.
C             Default: SPTDIV=1.
C     SPTSUB=real... Number of subintervals to be marked in each
C             interval.  Must be SPTSUB.GT.0.
C             Default: SPTSUB=1.
C     SPVRED=real... Reduction velocity.  If non-zero, the time at each
C             receiver is reduced by the value of RR/SPVRED, where RR is
C             the hypocentral distance.
C             No time reduction is applied if SPVRED=0.
C             Default: SPVRED=0.
C Data for the distance axis (horizontal):
C     KODESP=integer... Specifies the distribution and description of
C             seismograms in the plot.  See the description of SPXMIN,
C             SPXMAX,SPYMIN,SPYMAX below.
C             Default: KODESP=0
C     SPXMIN=real, SPXMAX=real, SPYMIN=real, SPYMAX=real:
C             For KODESP=0: Horizontal axis is divided into NREC+1
C               intervals, where NREC is the number of seismograms to be
C               plotted, see (1).  The I-th seismogram is plotted at the
C               endpoint of the I-th interval.  Horizontal axis has no
C               description, SPXMIN,SPXMAX,SPYMIN,SPYMAX have no
C               meaning.
C             For KODESP=1: Horizontal axis represents the profile with
C               endpoints (X1,X2)=(SPXMIN,SPYMIN) and
C               (X1,X2)=(SPXMAX,SPYMAX),
C               situated in a horizontal plane.  The seismograms are
C               located at the orthogonal projections of the receivers
C               onto the profile.  If SPXDIV.GT.0., the horizontal axis
C               is supplemented with the values of the X1 coordinates.
C             For KODESP=2: Horizontal axis represents the profile with
C               endpoints (X1,X2)=(SPXMIN,SPYMIN) and
C               (X1,X2)=(SPXMAX,SPYMAX),
C               situated in a horizontal plane.  The seismograms are
C               located at the orthogonal projections of the receivers
C               onto the profile.  If SPXDIV.GT.0., the horizontal axis
C               is supplemented with the values of both X1 and X2
C               coordinates.
C             For KODESP=3: Horizontal axis represents the vertical
C               profile with endpoints X3=SPXMIN and X3=SPXMAX.
C               The seismograms are located at the horizontal
C               projections of the receivers onto the profile.
C               If SPXDIV.GT.0., the horizontal axis is supplemented
C               with the values of X3 coordinate.
C               SPYMIN and SPYMAX have no meaning.
C             For KODESP=4: Horizontal axis represents the hypocentral
C               distance with endpoints RR=SPXMIN and RR=SPXMAX.
C               The seismograms are located according to the hypocentral
C               distances the receivers.
C               If SPXDIV.GT.0., the horizontal axis is supplemented
C               with the values of the hypocentral distance.
C               SPYMIN and SPYMAX have no meaning.
C             Default: SPXMIN=0., SPXMAX=1., SPYMIN=0., SPYMAX=0.
C     SPXLEN=real... Length of the horizontal axis in cm.
C             Default: SPXLEN=FLOAT(NREC+1), where NREC is the number of
C                                           seismograms to be plotted.
C     SPXDIV=real... ABS(SPXDIV) is the number of intervals along the
C               horizontal axis, starting at the left.
C               Must be SPXDIV.NE.0.
C             SPXDIV.GT.0.: the marks at the endpoints of intervals will
C               be supplemented with the corresponding values.
C             SPXDIV.LT.0.: horizontal axis will have no description.
C             Default: SPXDIV=1.
C     SPXSUB=real... Number of subintervals to be marked in each
C             interval.  Must be SPXSUB.GT.0.
C             Default: SPXSUB=1.
C Amplitude scaling:
C     NORMSP=integer... Type of amplitude scaling:
C             NORMSP=0: Maximum amplitude at each trace normalized to
C               the given value.
C               Simple and universal option.  No other amplitude scaling
C               parameter usually needs to be specified for this option,
C               although the input parameters SPAMP and SPAMPi enable
C               for possible amplitude scaling.
C               Amplitude scale AA is
C                 AA=SPAMPi*XD/AMAX
C               where AMAX is the maximum amplitude in each seismogram
C               and XD is the average distance between seismograms, i.e.
C               XD=SPXLEN/(NREC+1), where NREC is the number of plotted
C               seismograms.  The receiver file given by input parameter
C               REC is thus required to determine NREC.
C             Otherwise: All seismograms have the same scaling or power
C               scaling.
C               Amplitude scale AA is
C                 AA=SPAMPi*(RR/SPDIST)**SPOWER
C               where SPDIST is the hypocentral distance of the receiver
C               under consideration.
C             Default: NORMSP=0.
C     SPAMP=real... Amplitude scale for all 3 components.
C             Default: SPAMP=1.
C     SPAMP1=real, SPAMP2=real, ..., SPAMP9=real... Amplitude scales
C             SPAMP individually set for optional input GSE files
C             given by parameters SS1, SS2, ..., SS9, respectively.
C             Defaults: SPAMP1=SPAMP, SPAMP2=SPAMP, ..., SPAMP9=SPAMP
C     SPAMPX1=real, SPAMPX2=real, SPAMPX3=real... Amplitude scale
C             multiplicative factors for individual components.
C             SPAMP will be multiplied by SPAMPX1 for component 1,
C             by SPAMPX2 for component 2, by SPAMPX3 for component 3.
C             Default: SPAMPX1=1, SPAMPX2=1, SPAMPX3=1
C     SPOWER=real... Exponent of the amplitude power scaling.
C             Need not be specified if power scaling is not required.
C             Has no meaning if NORMSP=0.
C             Default: SPOWER=0.
C     SPOWER1=real, SPOWER2=real, ..., SPOWER9=real... Exponents of the
C             amplitude power scaling for optional input GSE files
C             given by parameters SS1, SS2, ..., SS9, respectively.
C             Defaults: SPOWER1=SPOWER, SPOWER2=SPOWER, ...,
C                       SPOWER9=SPOWER
C     SPDIST=real... Reference hypocentral distance for the amplitude
C             power scaling.
C             Need not be specified if power scaling is not required.
C             Has no meaning if NORMSP=0 or SPOWER=0.
C             Default: SPDIST=1.
C     SPEXP=real, SPEXPT=real... Exponential scaling of seismograms with
C             respect to time.  The amplitude scale is multiplied by the
C             factor of EXP(SPEXP*(t-SPEXPT)).
C             Example 1:  Exponential scaling with SPEXP=pi*F/Q, where
C               F is the dominant frequency and Q is the quality factor,
C               may roughly compensate for attenuation when plotting
C               late arrivals.
C             Example 2:  Exponential scaling with SPEXP equal to half
C               the sum of the positive Lyapunov exponents may
C               compensate for large geometrical spreading when plotting
C               late arrivals.
C             Default: SPEXP=0., SPEXPT=0.
C     SPEXP1=real, SPEXP2=real, ..., SPEXP9=real... Exponential scaling
C             set individually for optional input GSE files given by
C             parameters SS1, SS2, ..., SS9, respectively.
C             Defaults: SPEXP1=SPEXP, SPEXP2=SPEXP, ..., SPEXP9=SPEXP
C Other data to control plotting:
C     SPCHRH=real... Character height in cm.
C             Default: SPCHRH=0.4
C     SPHIWI=real... Width of the highlighted area when plotting travel
C             times.
C             Default: SPHIWI=XD,
C               where XD is the average distance between seismograms,
C               i.e. XD=SPXLEN/(NREC+1), where NREC is the number of
C               plotted seismograms, see the amplitude scaling.
C
C                                                     
C Input file GSE with the seismograms to plot:
C     File in the GSE data exchange format, see the description in file
C     'gse.for'.
C     The 'sp.for' program is looking in the comment lines of the
C     waveform identification section for the hypocentral coordinates
C     identified by strings 'XS1 ', 'XS2 ' and 'XS3 '.  If they are
C     present, they overwrite the values from the optional file SRC.
C     Description of format GSE
C
C                                                     
C Input file SRC with the hypocentral coordinates:
C (1) /
C     None to 20 character strings terminated by a slash.
C (2) 'SRCNAME',XS1,XS2,XS3,/
C     'SRCNAME'... Name of the source.  Not used.
C     XS1,XS2,XS3... Coordinates of the hypocentre.
C     Default: XS1=0., XS2=0., XS3=0.
C
C                                                     
C Input file REC containing receiver coordinates:
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) Several times (2.1):
C (2.1) 'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR),/
C     'NAMER(IR)'... String reserved for the name of the receiver.
C             No meaning here, anything in apostrophes may be submitted.
C     X1REC(IR),X2REC(IR),X3REC(IR)... Coordinates of the receiver.
C             The coordinates need not be present in the file.  It may
C             thus be comfortable to omit them if preparing the list for
C             the selection of particular receivers for seismogram
C             plotting.
C     Default: X1REC(IR)=0, X2REC(IR)=0, X3REC(IR)=0.
C (3) / (a slash).
C
C                                                  
C Input file SPHILI containing travel times to highlight:
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) Several times (2.1):
C (2.1) 'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR),TT,TTERR,K1,K2,/
C     'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR)... Same meaning as in
C             the receiver file above.
C     TT...   Travel time.  If given, it will be plotted as horizontal
C             line of width SPHIWI and colour K1.  It is plotted after
C             the frame and the corresponding error bar, but before
C             seismograms.
C     TTERR...Travel time error.  If given, it will be plotted as
C             a solid rectangle of width SPHIWI and colour K2. It is
C             plotted after the frame but before the corresponding
C             travel time and before seismograms.
C     K1...   Colour to plot the travel time.  The travel time is not
C             plotted if K1 is not positive.
C     K2...   Colour to plot the error bar.  The error bar is not
C             plotted if K2 is negative.  If K2=0 (white), the contour
C             of the rectangle if plotted in colour 1 (usually black).
C     Default: X1REC(IR)=0, X2REC(IR)=0, X3REC(IR)=0.,K1=1,K2=0
C (3) / (a slash).
C
C.......................................................................
C                                                
C This Fortran77 file consists of the following external procedures:
C     SP...   Main program to read and plot the seismograms.
C             SP
C     FRAME...Subroutine called by the main to plot the rectangular
C             frame around the seismograms and supplement it with simple
C             descriptions.
C             FRAME
C
C Other external procedures required:
C     RGSE1,RGSE2... Subroutines of the Fortran 77 file 'gse.for'
C             (package MODEL), designed to read seismograms in the GSE
C             data exchange format.
C     PLOTS,PLOT,SYMBOL,NUMBER... CALCOMP plotting subroutines.  For
C             example, Fortran 77 routines of file 'calcops.for'
C             (package MODEL) may be used to generate seismogram plots
C             in the PostScript files.
C
C=======================================================================
C
C     
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working array:
      INTEGER MPTS
      PARAMETER (MPTS=MRAM)
      REAL SS(MPTS)
      EQUIVALENCE (SS,RAM)
C
C-----------------------------------------------------------------------
C
C     Input and output data filenames and logical unit numbers:
      INTEGER LU
      PARAMETER (LU=1)
      CHARACTER*80 FILSEP,FILESS(0:9),FILSRC,FILREC,FILHIL,FILEPS(3)
C
      REAL UNDEF
      PARAMETER (UNDEF=-999999.)
      INTEGER KOLOR(0:9)
      REAL SPAMP(0:9),SPAMPX(3),SPOWER(0:9),SPEXP(0:9),XPTS(4),YPTS(4)
C
      CHARACTER*1  CHAN,TEXT
      CHARACTER*80 COMLIN
C
C     Data specifying labels in the plot:
      CHARACTER*20 KXTEXT(5)
      INTEGER KX(5)
C
C     Source coordinates transferred through the GSE file:
      INTEGER NCOM
      PARAMETER (NCOM=3)
      REAL VCOM(NCOM),DEF
C
C     List of receivers for plotting:
      INTEGER MREC,NREC
      PARAMETER (MREC=1000)
      CHARACTER*6 STAT,REC(MREC)
C
      DATA KX /0,1,1,1,19/
      DATA KXTEXT/' ','X1','X1','X3','HYPOCENTRAL DISTANCE'/
C
C.......................................................................
C
C     Reading name of SEP file with input data:
      FILSEP=' '
      WRITE(*,'(A)') '+SP: Enter input filename '
      READ(*,*) FILSEP
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU,FILSEP)
      ELSE
C       SP-04
        CALL ERROR('SP-04: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
      WRITE(*,'(A)') '+SP: Working...           '
C
C     Input and output filenames:
      CALL RSEP3T('SS' ,FILESS(0),'ss.gse')
      CALL RSEP3T('SS1',FILESS(1),' ')
      CALL RSEP3T('SS2',FILESS(2),' ')
      CALL RSEP3T('SS3',FILESS(3),' ')
      CALL RSEP3T('SS4',FILESS(4),' ')
      CALL RSEP3T('SS5',FILESS(5),' ')
      CALL RSEP3T('SS6',FILESS(6),' ')
      CALL RSEP3T('SS7',FILESS(7),' ')
      CALL RSEP3T('SS8',FILESS(8),' ')
      CALL RSEP3T('SS9',FILESS(9),' ')
      CALL RSEP3T('SP1',FILEPS(1),'ss1.ps')
      CALL RSEP3T('SP2',FILEPS(2),'ss2.ps')
      CALL RSEP3T('SP3',FILEPS(3),'ss3.ps')
C
C     Hypocentral coordinates:
      VCOM(1)=0.
      VCOM(2)=0.
      VCOM(3)=0.
      CALL RSEP3T('SRC',FILSRC,' ')
      IF(FILSRC.NE.' ') THEN
        OPEN(LU,FILE=FILSRC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        READ(LU,*) STAT,VCOM(1),VCOM(2),VCOM(3)
        CLOSE(LU)
      END IF
C
C     Number and list of receivers:
      NREC=999
      CALL RSEP3T('REC',FILREC,' ')
      IF(FILREC.NE.' ') THEN
        OPEN(LU,FILE=FILREC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        DO 1 I=1,MREC
          REC(I)='$$$$$$'
          READ(LU,*,END=2) REC(I)
          IF(REC(I).EQ.'$$$$$$') THEN
            GO TO 2
          END IF
    1   CONTINUE
C         SP-01
          CALL ERROR('SP-01: Array dimension MREC small')
    2   CONTINUE
        NREC=I-1
        CLOSE(LU)
      END IF
C
C     Colours:
      CALL RSEP3I('KOLOR ',KOLOR(0),1)
      CALL RSEP3I('KOLOR1',KOLOR(1),2)
      CALL RSEP3I('KOLOR2',KOLOR(2),3)
      CALL RSEP3I('KOLOR3',KOLOR(3),4)
      CALL RSEP3I('KOLOR4',KOLOR(4),5)
      CALL RSEP3I('KOLOR5',KOLOR(5),6)
      CALL RSEP3I('KOLOR6',KOLOR(6),7)
      CALL RSEP3I('KOLOR7',KOLOR(7),8)
      CALL RSEP3I('KOLOR8',KOLOR(8),9)
      CALL RSEP3I('KOLOR9',KOLOR(9),10)
C
C     Initial values for plotting frame:
C     Time axis:
      CALL RSEP3R('SPTMIN',SPTMIN, 0.)
      CALL RSEP3R('SPTMAX',SPTMAX, 1.)
      CALL RSEP3R('SPTLEN',SPTLEN,10.)
      CALL RSEP3R('SPTDIV',SPTDIV, 1.)
      CALL RSEP3R('SPTSUB',SPTSUB, 1.)
      CALL RSEP3R('SPVRED',SPVRED, 0.)
C     Distance axis:
      CALL RSEP3I('KODESP',KODESP, 0 )
      CALL RSEP3R('SPXMIN',SPXMIN, 0.)
      CALL RSEP3R('SPXMAX',SPXMAX, 1.)
      CALL RSEP3R('SPYMIN',SPYMIN, 0.)
      CALL RSEP3R('SPYMAX',SPYMAX, 0.)
      CALL RSEP3R('SPXLEN',SPXLEN,FLOAT(NREC+1))
      CALL RSEP3R('SPXDIV',SPXDIV, 1.)
      CALL RSEP3R('SPXSUB',SPXSUB, 1.)
C     Characters:
      CALL RSEP3R('SPCHRH',SPCHRH, 0.4)
C
      SCY= SPTLEN/(SPTMAX-SPTMIN)
      XD = SPXLEN/FLOAT(NREC+1)
      XA = SPXMAX-SPXMIN
      YA = SPYMAX-SPYMIN
C
C     Amplitude scaling:
      CALL RSEP3I('NORMSP',NORMSP,0)
      CALL RSEP3R('SPAMP ',SPAMP(0),1.)
      CALL RSEP3R('SPAMP1',SPAMP(1),SPAMP(0))
      CALL RSEP3R('SPAMP2',SPAMP(2),SPAMP(0))
      CALL RSEP3R('SPAMP3',SPAMP(3),SPAMP(0))
      CALL RSEP3R('SPAMP4',SPAMP(4),SPAMP(0))
      CALL RSEP3R('SPAMP5',SPAMP(5),SPAMP(0))
      CALL RSEP3R('SPAMP6',SPAMP(6),SPAMP(0))
      CALL RSEP3R('SPAMP7',SPAMP(7),SPAMP(0))
      CALL RSEP3R('SPAMP8',SPAMP(8),SPAMP(0))
      CALL RSEP3R('SPAMP9',SPAMP(9),SPAMP(0))
      CALL RSEP3R('SPAMPX1',SPAMPX(1),1.)
      CALL RSEP3R('SPAMPX2',SPAMPX(2),1.)
      CALL RSEP3R('SPAMPX3',SPAMPX(3),1.)
      CALL RSEP3R('SPDIST' ,SPDIST,1.)
      CALL RSEP3R('SPOWER' ,SPOWER(0),0.)
      CALL RSEP3R('SPOWER1',SPOWER(1),SPOWER(0))
      CALL RSEP3R('SPOWER2',SPOWER(2),SPOWER(0))
      CALL RSEP3R('SPOWER3',SPOWER(3),SPOWER(0))
      CALL RSEP3R('SPOWER4',SPOWER(4),SPOWER(0))
      CALL RSEP3R('SPOWER5',SPOWER(5),SPOWER(0))
      CALL RSEP3R('SPOWER6',SPOWER(6),SPOWER(0))
      CALL RSEP3R('SPOWER7',SPOWER(7),SPOWER(0))
      CALL RSEP3R('SPOWER8',SPOWER(8),SPOWER(0))
      CALL RSEP3R('SPOWER9',SPOWER(9),SPOWER(0))
      CALL RSEP3R('SPEXP ',SPEXP(0),0.)
      CALL RSEP3R('SPEXP1',SPEXP(1),SPEXP(0))
      CALL RSEP3R('SPEXP2',SPEXP(2),SPEXP(0))
      CALL RSEP3R('SPEXP3',SPEXP(3),SPEXP(0))
      CALL RSEP3R('SPEXP4',SPEXP(4),SPEXP(0))
      CALL RSEP3R('SPEXP5',SPEXP(5),SPEXP(0))
      CALL RSEP3R('SPEXP6',SPEXP(6),SPEXP(0))
      CALL RSEP3R('SPEXP7',SPEXP(7),SPEXP(0))
      CALL RSEP3R('SPEXP8',SPEXP(8),SPEXP(0))
      CALL RSEP3R('SPEXP9',SPEXP(9),SPEXP(0))
      CALL RSEP3R('SPEXPT',SPEXPT,0.)
C
C     Higlighting given areas (e.g., travel times with error bars):
      CALL RSEP3T('SPHILI',FILHIL,' ')
      CALL RSEP3R('SPHIWI',SPHIWI,XD)
C
C.......................................................................
C
C     Loop over 3 seismogram components:
      DO 99 KOMP=1,3
        IF(FILEPS(KOMP).NE.' ') THEN
C
C         Initialization of CALCOMP:
          CALL PLOTN(FILEPS(KOMP),0)
          CALL PLOTS(0,0,0)
C
C         Plotting frame:
          CALL NEWPEN(1)
          IX = KODESP
          IT = -1
          IF(KODESP.GE.3)  IX=1
          IF(SPXDIV.LT.0.) IX=0
          IF(SPTDIV.LT.0.) IT=0
          KX1= KX(KODESP+1)+1
          CALL FRAME(SPXLEN,SPTLEN,ABS(SPXDIV),SPXSUB,
     *                             ABS(SPTDIV),SPTSUB,0,IX,IT,SPCHRH,
     *               SPXMIN,SPXMAX,KX1,KXTEXT(KODESP+1),
     *               SPYMIN,SPYMAX,2,'X2',
     *               SPTMIN,SPTMAX,4,'TIME',
     *               1,' ',0,0.,0,1,' ',0,0.,0,1,' ',1,' ')
C
C         Higlighting given areas (e.g., travel times with error bars):
          IF(FILHIL.NE.' ') THEN
            OPEN(LU,FILE=FILHIL,STATUS='OLD')
            READ(LU,*) (TEXT,I=1,20)
C           Loop for areas to highlight:
   10       CONTINUE
              STAT='$$$$$$'
              T0=UNDEF
              TD=UNDEF
              K1=1
              K2=0
              READ(LU,*,END=19) STAT,XR,YR,ZR,T0,TD,K1,K2
              IF(STAT.EQ.'$$$$$$') THEN
                GO TO 19
              END IF
              IF(T0.EQ.UNDEF) THEN
                K1=-1
                K2=-1
              END IF
              IF(TD.EQ.UNDEF) THEN
                K2=-1
              END IF
C             Determining the receiver:
              IF(FILREC.NE.' ') THEN
C               Loop for receivers
                DO 11 IREC=1,NREC
                  IF(REC(IREC).EQ.STAT) GO TO 12
   11           CONTINUE
                GO TO 10
              END IF
   12         CONTINUE
              IF(SPVRED.NE.0.) THEN
                RR=SQRT((XR-VCOM(1))**2+(YR-VCOM(2))**2+(ZR-VCOM(3))**2)
                T0=T0-RR/SPVRED
              END IF
              IF(KODESP.EQ.0) THEN
                IF(FILREC.EQ.' ') THEN
C                 SP-02
                  CALL ERROR('SP-02: No receiver file specified')
C                 For KODESP=0, filename REC must be specified in the
C                 input data.
                END IF
                X=XD*FLOAT(IREC)
              ELSE IF(KODESP.EQ.1.OR.KODESP.EQ.2) THEN
                X=SPXLEN*((XR-SPXMIN)*XA+(YR-SPYMIN)*YA)/(XA*XA+YA*YA)
              ELSE IF(KODESP.EQ.3) THEN
                X=SPXLEN*(ZR-SPXMIN)/XA
              ELSE
                X=SPXLEN*(RR-SPXMIN)/XA
              END IF
C             Plotting the highlight:
              IF(K2.GE.0) THEN
                XPTS(1)=X-0.5*SPHIWI
                XPTS(2)=X+0.5*SPHIWI
                XPTS(3)=X+0.5*SPHIWI
                XPTS(4)=X-0.5*SPHIWI
                YPTS(1)=SCY*(T0-TD-SPTMIN)
                YPTS(2)=SCY*(T0-TD-SPTMIN)
                YPTS(3)=SCY*(T0+TD-SPTMIN)
                YPTS(4)=SCY*(T0+TD-SPTMIN)
                IF(K2.GT.0) THEN
                  CALL NEWPEN(K2)
                  CALL FILL(XPTS,YPTS,4)
                ELSE
                  CALL NEWPEN(1)
                  CALL PLOT(XPTS(1),YPTS(1),3)
                  CALL PLOT(XPTS(2),YPTS(2),2)
                  CALL PLOT(XPTS(3),YPTS(3),2)
                  CALL PLOT(XPTS(4),YPTS(4),2)
                END IF
              END IF
              IF(K1.GT.0) THEN
                CALL NEWPEN(K1)
                CALL PLOT(X-0.5*SPHIWI,SCY*(T0-SPTMIN),3)
                CALL PLOT(X+0.5*SPHIWI,SCY*(T0-SPTMIN),2)
              END IF
            GO TO 10
   19       CONTINUE
C           Closing highlighting file
            CLOSE(LU)
          END IF
C
C         Plotting seismograms:
C         Loop for GSE files
          DO 49 ISS=0,9
            IF(FILESS(ISS).NE.' ') THEN
              CALL NEWPEN(KOLOR(ISS))
C
C             Opening input GSE file with the seismograms:
              OPEN(LU,FILE=FILESS(ISS),STATUS='OLD')
              CALL RGSE1(LU,TEXT)
              IREC=0
C
C             Loop for seismograms
   20         CONTINUE
   21           CONTINUE
C
C               Selecting the component:
   22           CONTINUE
                  CALL RGSE2
     *                     (LU,STAT,CHAN,I,XR,YR,ZR,T0,TD,NPTS,MPTS,SS)
                  IF(NPTS.LE.-1) THEN
C                   End of the GSE file
                    GO TO 40
                  END IF
                IF(I.NE.KOMP) GO TO 22
C
C               Determining the receiver:
                IF(FILREC.EQ.' ') THEN
                  IREC=IREC+1
                ELSE
C                 Loop for receivers
                  DO 31 IREC=1,NREC
                    IF(REC(IREC).EQ.STAT) GO TO 32
   31             CONTINUE
                  GO TO 21
                END IF
   32           CONTINUE
C
C               Plotting the seismogram:
   33           CONTINUE
                  CALL RGSE2C(COMLIN,*34)
                  CALL RSEP2(COMLIN)
                GO TO 33
   34           CONTINUE
                DEF=VCOM(1)
                CALL RSEP3R('XS1',VCOM(1),DEF)
                DEF=VCOM(2)
                CALL RSEP3R('XS2',VCOM(2),DEF)
                DEF=VCOM(3)
                CALL RSEP3R('XS3',VCOM(3),DEF)
                RR=SQRT((XR-VCOM(1))**2+(YR-VCOM(2))**2+(ZR-VCOM(3))**2)
C
                IF(SPEXP(ISS).NE.0.) THEN
                  DO 35 I=1,NPTS
                   SS(I)=SS(I)*EXP(SPEXP(ISS)*(T0+TD*FLOAT(I-1)-SPEXPT))
   35             CONTINUE
                END IF
                IF(SPVRED.NE.0.) THEN
                  T0=T0-RR/SPVRED
                END IF
C
                IF(NORMSP.EQ.0) THEN
                  AA=ABS(SS(1))
                  DO 36 I=2,NPTS
                    AA=AMAX1(ABS(SS(I)),AA)
   36             CONTINUE
                  AA=SPAMP(ISS)*SPAMPX(KOMP)*XD/AA
                ELSE
                  IF(SPOWER(ISS).EQ.0.) THEN
                    AA=SPAMP(ISS)*SPAMPX(KOMP)
                  ELSE
                    AA=SPAMP(ISS)*SPAMPX(KOMP)
     *                                       *((RR/SPDIST)**SPOWER(ISS))
                  END IF
                END IF

                IF(KODESP.EQ.0) THEN
                  IF(FILREC.EQ.' ') THEN
C                   SP-03
                    CALL ERROR('SP-03: No receiver file specified')
C                   For KODESP=0, filename REC must be specified in the
C                   input data.
                  END IF
                  X=XD*FLOAT(IREC)
                ELSE IF(KODESP.EQ.1.OR.KODESP.EQ.2) THEN
                  X=SPXLEN*((XR-SPXMIN)*XA+(YR-SPYMIN)*YA)/(XA*XA+YA*YA)
                ELSE IF(KODESP.EQ.3) THEN
                  X=SPXLEN*(ZR-SPXMIN)/XA
                ELSE
                  X=SPXLEN*(RR-SPXMIN)/XA
                END IF
                Y = 0.
                CALL PLOT(X,Y,3)
C
                DO 38 I=1,NPTS
                  T = T0+TD*FLOAT(I-1)
                  IF (SPTMIN.LE.T.AND.T.LE.SPTMAX) THEN
                    A = X-AA*SS(I)
                    Y = SCY*(T-SPTMIN)
                    IF (T.LT.SPTMIN+TD) THEN
C                     Straight line before the seismogram
                      CALL PLOT(X,Y,2)
                    END IF
                    CALL PLOT(A,Y,2)
                  END IF
   38           CONTINUE
C
C               Straight line after the seismogram
                CALL PLOT(X,Y,2)
                Y = SPTLEN
                CALL PLOT(X,Y,2)
              GO TO 20
   40         CONTINUE
C             End of loop for receivers
C
C             Closing GSE file
              CLOSE(LU)
            END IF
   49     CONTINUE
C         End of loop for GSE files
C
          CALL PLOT(0.,0.,999)
        END IF
   99 CONTINUE
      WRITE(*,'(A)') '+SP: Done.                '
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE FRAME(XX,YY,XM,XN,YM,YN,MM,IX,IY,HEIGHT,
     *                 XL1,XR1,LX1,KX1,XL2,XR2,LX2,KX2,YL,YR,LY,KY,
     *                 L1,K1,M1,A1,N1,L2,K2,M2,A2,N2,L3,K3,L4,K4)
C
      CHARACTER*(*) KX1,KX2,KY,K1,K2,K3,K4
      REAL XX,YY,XM,XN,YM,YN,HEIGHT,XL1,XR1,XL2,XR2,YL,YR,A1,A2
      INTEGER MM,IX,IY,LX1,LX2,LY,L1,M1,N1,L2,M2,N2,L3,L4
C
C Input:
C     XX...   Length of the horizontal axis in cm.
C     YY...   Length of the vertical axis in cm.
C     XM...   The number of intervals along the horizontal axis,
C             starting at the left.  The intervals are marked with
C             long ticks and are supplemented with the numerical
C             values if IX is positive.  XM must be positive.
C     XN...   Number of subintervals to be marked in each interval
C             with short ticks.  XN must be positive.
C     YM...   The number of intervals along the vertical axis,
C             starting at the bottom.  The intervals are marked with
C             long ticks and are supplemented with the numerical
C             values if IY is positive.  YM must be positive.
C     YN...   Number of subintervals to be marked in each interval
C             with short ticks.  YN must be positive.
C     MM...   If negative, the top line of the frame is not plotted.
C     IX...   IX=0: No labeling of the horizontal axis.
C             IX=1: Labeling of the horizontal axis with the first
C                   variable only.
C             IX=2: Labeling of the horizontal axis with both variables.
C     IY...   IY=0: No labeling of the vertical axis.
C             IY=1: Labeling of the vertical axis.
C     HEIGHT... Height of the characters in cm.
C     XL1,XR1... Values of the first variable along the horizontal axis,
C             corresponding to left-hand and right-hand corners.
C     LX1...  Length of string KX1.
C     KX1...  First label of the horizontal axis.  String to be written
C             below the horizontal axis.
C     XL2,XR2... Values of the second variable along the horizontal
C             axis, corresponding to left-hand and right-hand corners.
C     LX2...  Length of string KX2.
C     KX2...  Second label of the horizontal axis.  String to be written
C             below the horizontal axis.
C     YL,YR...Values of the variable along the vertical axis,
C             corresponding to left-hand and right-hand corners.
C     LY...   Length of string KY.
C     KY...   Label of the vertical axis.  String to be written to the
C             left of the vertical axis.
C     L1...   Length of string K1.
C     K1,A1...String and the number to be written above the frame,
C             starting 0.5*HEIGHT above the left top corner.
C     M1...   Number A1 is written only if M2 is positive.
C     N1...   Number of decimal places of A1.
C     L2...   Length of string K1.
C     K2,A2...String and the number to be written above the frame,
C             ending 0.5*HEIGHT above the right top corner.
C     M2...   Width of A2 in characters.  A2 is written only if M2 is
C             positive.
C     N2...   Number of decimal places of A2.
C     L2...   Length of string K2.
C     L3...   Length of string K3.
C     K3...   String to be written below the frame, starting 5.3*HEIGHT
C             below the left bottom corner.
C     L4...   Length of string K4.
C     K4...   String to be written below the frame, starting 7.0*HEIGHT
C             below the left bottom corner.
C
C-----------------------------------------------------------------------
C
      HEIGH0=.215*HEIGHT
      HEIGH2=.500*HEIGHT
C
C     Initial values for plotting frame:
      I0= 0
      JX= IABS(IX)-1
      MX= INT(XM*XN+0.001)
      NX= INT(XN+0.5)
      MY= INT(YM*YN+0.001)
      NY= INT(YN+0.5)
      XD= XX/XM/XN
      YD= YY/YM/YN
      XH1= (XR1-XL1)/XM
      XH2= (XR2-XL2)/XM
      YH = (YR -YL )/YM
C
C     Plotting border 29.7*21.0 cm (landscape):
C     CALL PLOT( 0. , 0. ,3)
C     CALL PLOT(29.7, 0. ,2)
C     CALL PLOT(29.7,21.0,2)
C     CALL PLOT( 0. ,21.0,2)
C     CALL PLOT( 0. , 0. ,2)
C
C     Plotting frame XX*YY cm centred on the A4 sheet:
C     Landscape
*     X = (29.7-XX)/2.
*     Y = (21.0-YY)/2.
C     Portrait
      X = (21.0-XX)/2.
      Y = (29.7-YY)/2.
C     Leaving 2 cm belts for description of axes
      IF(IX.GE.0) THEN
        X=X+2.5*HEIGHT
      END IF
      IF (IY.NE.0) THEN
        Y=Y+2.5*HEIGHT
      END IF
C     Shifting the origin and plotting the frame
      CALL PLOT(X,Y,-3)
      CALL PLOT(XX,0.,2)
      CALL PLOT(XX,YY,2)
      IF(MM.GE.0) THEN
        I = 2
      ELSE
        I = 3
      END IF
      CALL PLOT(0.,YY,I)
      CALL PLOT(0.,0.,2)
C
C     Description of axes:
C
C     Bottom horizontal axis:
      X = 0.
      X1= XL1-XH1
      X2= XL2-XH2
      DO 16 I=I0,MX
        IF (MOD(I,NX).NE.0) THEN
          A = 0.1
        ELSE
          A = 0.2
          IF (JX.GE.0) THEN
            X1= X1+XH1
            X2= X2+XH2
            IF(IX.GE.0.OR.MOD(I,MX).NE.0) THEN
C             Determination of the number of decimal places:
ccc           J = INT(.99-ALOG10(AMAX1(ABS(XH1),ABS(XH2),0.001)))
              DO 11 J=0,5
                IF(AMOD(ABS(XH1)+0.000001,0.1**J).LT.0.000002) THEN
                  GO TO 12
                END IF
   11         CONTINUE
   12         CONTINUE
C             J is the preferable number of decimal places.
              JMAX=MAX0(INT(ALOG10(ABS(X1)+0.5*0.1**J))+1,1)
              IF(X1.LT.0.) JMAX=JMAX+1
C             JMAX is the number of digits left to the decimal point.
              IF (JX.GT.0) THEN
                DO 13 J=J,5
                  IF(AMOD(ABS(XH2)+0.000001,0.1**J).LT.0.000002) THEN
                    GO TO 14
                  END IF
   13           CONTINUE
   14           CONTINUE
C               J is the preferable number of decimal places.
                JMAX=MAX0(INT(ALOG10(ABS(X2)+0.5*0.1**J))+1,JMAX)
                IF(X1.GE.0..AND.X2.LT.0.) JMAX=JMAX+1
C               JMAX is the number of digits left to the decimal point.
              END IF
              JMAX=INT(XX/XM/HEIGHT)-JMAX-1
C             JMAX is the maximum number of decimal places.
              J=MIN0(J,JMAX)
              IF(J.LE.0) THEN
                J=-1
              END IF
C             J is the number of decimal places.
C
              B = X-( .5*FLOAT(1+JX)*AINT(ALOG10(ABS(X1)+.5))+0.285
     *                  -FLOAT(JX)+.5*FLOAT(J+1) )*HEIGHT
              IF(X1.LT.0.) THEN
                B=B-HEIGHT
              END IF
              CALL NUMBER(B,-1.8*HEIGHT,HEIGHT,X1,0.,J)
              IF (JX.GT.0) THEN
                B = X-HEIGHT*AINT(ALOG10(ABS(X2)+.5))+.715*HEIGHT
                IF(X2.LT.0.) B=B-HEIGHT
                CALL NUMBER(B,-3.3*HEIGHT,HEIGHT,X2,0.,J)
              END IF
            END IF
          END IF
        END IF
        IF (MOD(I,MX).NE.0) THEN
          CALL PLOT(X,0.,3)
          CALL PLOT(X,A ,2)
        END IF
        X = X+XD
   16 CONTINUE
      IF (JX.EQ.0) THEN
        B = XX-XX/XM/2.-HEIGH2*FLOAT(LX1)+HEIGH0
        CALL SYMBOL(B,-3.3*HEIGHT,HEIGHT,KX1,0.,LX1)
      ELSE IF (JX.GT.0) THEN
        CALL SYMBOL(XX+2.715*HEIGHT,-1.8*HEIGHT,HEIGHT,KX1,0.,LX1)
        CALL SYMBOL(XX+2.715*HEIGHT,-3.3*HEIGHT,HEIGHT,KX2,0.,LX2)
      END IF
C
C     Right-hand vertical axis:
      Y = 0.
      M = MY
      DO 23 I=1,M
        Y = Y+YD
        A = 0.1
        IF (MOD(I,NY).EQ.0) THEN
          A = 0.2
        END IF
        CALL PLOT(XX,Y,3)
        CALL PLOT(XX-A,Y,2)
   23 CONTINUE
C
C     Top horizontal axis:
      IF (MM.GE.0) THEN
        X = XX
        M = MX-1
        DO 33 I=1,M
          X = X-XD
          IF (MOD(I,NX).NE.0) THEN
            A = 0.1
          ELSE
            A = 0.2
          END IF
          CALL PLOT(X,YY,3)
          CALL PLOT(X,YY-A,2)
   33   CONTINUE
      END IF
C
C     Left-hand vertical axis:
      Y = 0.
      Y1= YL-YH
      DO 45 I=I0,MY
        IF (MOD(I,NY).NE.0) THEN
          A = 0.1
        ELSE
          A = 0.2
          IF (IY.NE.0) THEN
            Y1= Y1+YH
C
C           Determination of the number of decimal places:
ccc         J = INT(.99-ALOG10(AMAX1(ABS(YH),0.001)))
            DO 41 J=0,5
              IF(AMOD(ABS(YH)+0.000001,0.1**J).LT.0.000002) THEN
                GO TO 42
              END IF
   41       CONTINUE
   42       CONTINUE
            IF(J.LE.0) THEN
              J=-1
            END IF
ccc         IF(J.LE.0.OR.IY.GT.0) THEN
ccc           J=IY
ccc         END IF
C           J is the number of decimal places.
C
            B = -( AINT(ALOG10(ABS(Y1)+.5))+2.785+FLOAT(J) )*HEIGHT
            IF(Y1.LT.0.) THEN
              B=B-HEIGHT
            END IF
            CALL NUMBER(B,Y-HEIGH2,HEIGHT,Y1,0.,J)
          END IF
        END IF
        IF (I.NE.0) THEN
          CALL PLOT(0.,Y,3)
          CALL PLOT(A,Y ,2)
        END IF
        Y = Y+YD
   45 CONTINUE
      IF (IY.NE.0) THEN
        A = -HEIGHT*FLOAT(LY)-1.285*HEIGHT
        IF (YM-FLOAT(MY/NY).GE.0.25) THEN
          B = YY-HEIGHT
        ELSE
          B = YY-YY/YM/2.-HEIGH2
        END IF
        CALL SYMBOL(A,B,HEIGHT,KY,0.,LY)
      END IF
C
C     Writing texts:
      CALL SYMBOL(HEIGH0,YY+HEIGH2,HEIGHT,K1,0.,L1)
      B = HEIGHT*FLOAT(L1)+1.215*HEIGHT
      IF(M1.GT.0) THEN
        CALL NUMBER(B,YY+HEIGH2,HEIGHT,A1,0.,N1)
      END IF
      B = XX-HEIGHT*FLOAT(L2+M2)-.785*HEIGHT
      CALL SYMBOL( B,YY+HEIGH2,HEIGHT,K2,0.,L2)
      B = XX-HEIGHT*FLOAT(M2)   +.215*HEIGHT
      IF(M2.GT.0) THEN
        CALL NUMBER(B,YY+HEIGH2,HEIGHT,A2,0.,N2)
      END IF
      CALL SYMBOL(HEIGH0,-5.3*HEIGHT,HEIGHT,K3,0.,L3)
      CALL SYMBOL(HEIGH0,-7.0*HEIGHT,HEIGHT,K4,0.,L4)
C
      RETURN
      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 'calcops.for'
C     calcops.for
C
C=======================================================================
C