C 
C Program SP (Seismogram Plotting) to plot seismograms previously stored
C in the GSE data exchange format.
C
C Version: 5.90
C Date: 2005, June 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             If the source names are specified in the comment lines
C               of the GSE file, the hypocentral coordinates of the
C               sources are taken, in the order of preference, from file
C               PTS (if the filename is specified), file SRC (if the
C               filename is specified), corresponding GSE file SS* (if
C               the coordinates are present), default to zeros.
C             If the source name is not specified in the comment lines
C               of the GSE file, the hypocentral coordinates of the
C               sources are taken, in the order of preference, from the
C               the first point in file SRC (if the filename is
C               specified), GSE file SS (if the coordinates are
C               present), default to zeros.
C             The coordinates of the receivers are taken, in the order
C               of preference, from file PTS (if the filename is
C               specified), file REC (if the filename is specified),
C               GSE file SS.
C             The source and receiver names cannot be longer than
C               6 characters.
C             Description of GSE file SS
C             Default: SS='ss.gse'
C     Hereinafter, # represents the value of integer constant MFILSS.
C             E.g., if MFILSS=12, SS# stands for SS12, KOLOR# stands for
C             KOLOR12, KOMP1# stands for KOMP112, etc.
C     SS1='string', SS2='string', ..., SS#='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.
C             The order of plottting is SS, SS1, SS2, ..., SS#,
C             considering just nonblank filenames.
C             The seismograms are plotted in colours given by parameters
C             KOLOR, KOLOR1, KOLOR2, ..., KOLOR#, 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=' ', ..., SS#=' '
C     SRC='string'... String with the name of the input data file
C             with the name(s) and coordinates of the hypocentre(s).
C             If specified, this file serves (a) to select the sources
C             for plotting seismograms, (b) to update the hypocentre
C             coordinates if PTS=' ' or if the source name is not
C             specified in the comment lines of the GSE file and file
C             PTS thus cannot be used.
C             If SRC is not blank and the source names are specified
C               comment lines of the GSE file, only seismograms
C               corresponding to the sources listed in file SRC will be
C               plotted, otherwise no selection according to the source
C               names is performed.
C             If the source names are specified in the comment lines
C               of the GSE file, the hypocentral coordinates of the
C               sources are taken, in the order of preference, from file
C               PTS (if the filename is specified), file SRC (if the
C               filename is specified), corresponding GSE file SS* (if
C               the coordinates are present), default to zeros.
C             If the source name is not specified in the comment lines
C               of the GSE file, the hypocentral coordinates are taken,
C               in the order of preference, from the the first point in
C               file SRC (if the filename is specified), corresponding
C               GSE file SS* (if the coordinates are present), default
C               to zeros.
C             The source names cannot be longer than 6 characters.
C             The hypocentral coordinates are used only for optional
C             travel-time reduction on a given velocity, or for
C             amplitude power scaling.  File SRC thus usually need not
C             be specified if the seismograms are generated by programs
C             'greenss.for' or 'ss.for', because those programs enter
C             the hypocentral coordinates directly to the comments of
C             the data section in GSE data exchange file.
C             Description of file SRC
C             Default: SRC=' '
C     REC='string'... String with the name of the input data file
C             containing the list of receivers.
C             If specified, this file serves (a) to select the receivers
C             for plotting seismograms, (b) to update the receiver
C             coordinates if PTS=' '.
C             If REC=' ' seismograms at all receivers will be plotted,
C               otherwise only seismograms at the receivers listed in
C               file REC will be plotted.
C             The coordinates of the receivers are taken, in the order
C               of preference, from file PTS (if the filename is
C               specified), file REC (if the filename is specified),
C               corresponding GSE file SS*.
C             File REC 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             Only the first 6 characters of receiver names are
C               significant.
C             The receiver names cannot be longer than 6 characters.
C             In most cases, file REC will be the same as for the
C             calculation of synthetic seismograms.
C             Parameter REC has to be specified and cannot be blank if
C             KODESP=0, see below.  If KODESP=0, the horizontal position
C             of a plotted seismogram corresponds to the position of the
C             corresponding receiver in file REC.
C             Description of file REC
C             Default: REC=' '
C     FTT='string'... String with the name of the input data file
C             with the list of travel times to be highlighted.
C             The travel times are specified in dependence on the
C             source and receiver names.  The coordinates of the source
C             and receivers are taken from file PTS.  If PTS=' ', the
C             coordinates of the source are taken from file SRC and the
C             coordinates of receivers are taken from file REC.
C             If PTS=' ' and SRC=' ', the source coordinates are assumed
C             zero.  If FTT.NE.' ' and PTS=' ', REC cannot be blank.
C             If FTT=' ', no travel times are highlighted.
C             Description of file FTT
C             Default: FTT=' '
C     PTS='string'... String with the name of the input data file
C             with the coordinates corresponding to the source and
C             receiver names.  Since this file is given just to specify
C             the coordinates, the coordinates from this file have the
C             highest priority.  This feature is especially useful when
C             changing the coordinate system.  The seismogram files
C             need not be changed and may preserve the old coordinates.
C             The same advantage applies to files SRC, REC and SPHILI.
C             The point names cannot be longer than 6 characters.
C             Description of file PTS
C             Default: PTS=' '
C     SPHILI='string'... String with the name of the input data file
C             with the list of times to be highlighted.
C             The times are specified in dependence on the names and
C             coordinates of receivers.
C             If REC is not blank and the receiver names are specified
C               in file SPHILI, only times corresponding to the
C               receivers listed in file REC will be plotted, otherwise
C               no selection according to the receiver names is
C               performed.
C             The coordinates of the receivers are taken, in the order
C               of preference, from file PTS (if filename PTS is
C               specified and the point names in file SPHILI are not
C               blank), from file SPHILI.
C             If SPHILI=' ', no times are highlighted.
C             Description of file SPHILI
C             Default: SPHILI=' '
C Names of output files:
C     SP1='string'... String with the name of the first output
C             PostScript file, usually contaning the plot of the first
C             component of the seismograms.
C             If blank, the file is not created.
C             Default: SP1='ss1.gse'
C     SP2='string'... String with the name of the second output
C             PostScript file, usually contaning the plot of the second
C             component of the seismograms.
C             If blank, the file is not created.
C             Default: SP2='ss2.gse'
C     SP3='string'... String with the name of the third output
C             PostScript file, usually contaning the plot of the third
C             component of the seismograms.
C             If blank, the file is not created.
C             Default: SP3='ss3.gse'
C Component selection (mostly not needed):
C     KOMP1=integer... Component of the seismograms of file given by
C             parameter SS, plotted into the output file given by
C             parameter SP1.
C             Default: KOMP1=1
C     KOMP11=integer, KOMP12=integer, ..., KOMP1#=integer...
C             Components of the seismograms of files given by
C             parameters SS1, SS2, SS3, ..., SS#, respectively,
C             plotted into the output file given by parameter SP1.
C             Defaults: KOMP11=KOMP1, KOMP12=KOMP1, ..., KOMP1#=KOMP1
C     KOMP2=integer... Component of the seismograms of file given by
C             parameter SS, plotted into the output file given by
C             parameter SP2.  Analogous to KOMP1.
C             Default: KOMP2=2
C     KOMP21=integer, KOMP22=integer, ..., KOMP2#=integer...
C             Analogous to KOMP11 to KOMP1#, but for file SP2.
C             Defaults equal the value of KOMP2.
C     KOMP3=integer... Component of the seismograms of file given by
C             parameter SS, plotted into the output file given by
C             parameter SP3.  Analogous to KOMP1.
C             Default: KOMP3=3
C     KOMP31=integer, KOMP32=integer, ..., KOMP3#=integer...
C             Analogous to KOMP11 to KOMP1#, but for file SP2.
C             Defaults equal the value of KOMP3.
C Data to control plotting:
C Colours:
C     KOLOR=positive integer, KOLOR1=positive integer,
C     KOLOR2=positive integer, ..., KOLOR#=positive integer... Colours
C             to plot seismograms of files SS,SS1, SS2, ..., SS#,
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, ..., KOLOR#=#+1
C     KOLORTT=positive integer... Colour to plot the travel times of
C             optional file FTT.  It is also used as the default colour
C             for optional file SPHILI.
C             Travel times are not plotted if KOLORTT is not positive.
C             Default: KOLORTT=1
C     KOLORTD=integer... Colour to plot the error bar of optional file
C             FTT.  It is also used as the default colour for optional
C             file SPHILI.
C             The error bar is not plotted if KOLORTD is negative.
C             If KOLORTD=0 (white), the contour of the rectangle is
C             plotted in colour KOLORTT.
C             Default: KOLORTD=0
C Data for the time axis (vertical):
C     SPTMIN=real... Time (or reduced time) corresponding to the bottom
C             of the seismogram plot.
C             Default: SPTMIN=0.
C     SPTMAX=real... Time (or reduced time) corresponding to the top of
C             the seismogram plot.
C             SPTMAX may be chosen smaller than SPTMIN to point the time
C             axis downwards.
C             Default: SPTMAX=1.
C     SPTLEN=positive 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=positive real... Number of subintervals to be marked in
C             each 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 represents the index of the
C               receiver, corresponding to the receiver position in
C               the file given by parameter REC.
C               SPXMIN and SPXMAX are the receiver indices corresponding
C               to the lef-hand and right-hand sides of the frame.
C               Example: For default values SPXMIN=0 and SPXMAX=NREC+1,
C                 where NREC is the number of receivers in file REC, the
C                 horizontal axis is divided into NREC+1 intervals.
C                 The seismogram at the Ith receiver is then plotted at
C                 the endpoint of the Ith interval.
C               SPYMIN and SPYMAX have no meaning.
C               If SPXDIV.GT.0., the horizontal axis is denoted by the
C               names of the receivers corresponding to the plotted
C               seismograms, otherwise it has no description.
C               Parameter REC has to be specified and cannot be blank in
C               this case.
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=NREC+1, SPYMIN=0., SPYMAX=0.,
C               where NREC is the number of receivers in file REC.
C     SPXMIN1=real, SPXMIN2=real, ..., SPXMIN#=real...
C             Analogous to SPXMIN, but for files SS1 to SS#.
C             Do not influence the description of the horizontal axis
C             nor highlighting of times.
C             Usually need not be specified.
C             Defaults equal the value of SPXMIN.
C     SPXMAX1=real, SPXMAX2=real, ..., SPXMAX#=real...
C             Analogous to SPXMAX, but for files SS1 to SS#.
C             Do not influence the description of the horizontal axis
C             nor highlighting of times.
C             Usually need not be specified.
C             Defaults equal the value of SPXMAX.
C     SPYMIN1=real, SPYMIN2=real, ..., SPYMIN#=real...
C             Analogous to SPYMIN, but for files SS1 to SS#.
C             Do not influence the description of the horizontal axis
C             nor highlighting of times.
C             Usually need not be specified.
C             Defaults equal the value of SPYMIN.
C     SPYMAX1=real, SPYMAX2=real, ..., SPYMAX#=real...
C             Analogous to SPYMAX, but for files SS1 to SS#.
C             Do not influence the description of the horizontal axis
C             nor highlighting of times.
C             Usually need not be specified.
C             Defaults equal the value of SPYMAX.
C     SPXLEN=positive real... Length of the horizontal axis in cm.
C             Default: SPXLEN=FLOAT(NREC+1), where NREC is the number of
C                      receivers in file REC.
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=positive real... Number of subintervals to be marked in
C             each interval.  Must be SPXSUB.GT.0.
C             Default: SPXSUB=1.
C Amplitude scaling:
C     NORMSP=integer... Type of amplitude scaling:
C             NORMSP.LT.0: Like NORMSP=0, but the maximum amplitude is
C               calculated only over the plotted part of seismogram.
C             NORMSP.EQ.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
C               receivers in file REC.  The receiver file given by input
C               parameter REC is thus required to determine NREC.
C             NORMSP.GT.0: All seismograms have the same scaling or
C               power 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, ..., SPAMP#=real... Amplitude scales
C             SPAMP individually set for optional input GSE files
C             given by parameters SS1, SS2, ..., SS#, respectively.
C             Defaults: SPAMP1=SPAMP, SPAMP2=SPAMP, ..., SPAMP#=SPAMP
C     SPAMPX1=real, SPAMPX2=real, SPAMPX3=real... Amplitude scale
C             multiplicative factors for individual output files SP1,
C             SP2 and SP3, usually corresponding to different
C             components.
C             SPAMP will be multiplied by SPAMPX1 for file SP1 (usually
C             component 1), by SPAMPX2 for file SP1 (usually component
C             2), by SPAMPX3 for file SP1 (usually 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.LE.0.
C             Default: SPOWER=0.
C     SPOWER1=real, SPOWER2=real, ..., SPOWER#=real... Exponents of the
C             amplitude power scaling for optional input GSE files
C             given by parameters SS1, SS2, ..., SS#, respectively.
C             Defaults: SPOWER1=SPOWER, SPOWER2=SPOWER, ...,
C                       SPOWER#=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.LE.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, ..., SPEXP#=real... Exponential scaling
C             set individually for optional input GSE files given by
C             parameters SS1, SS2, ..., SS#, respectively.
C             Defaults: SPEXP1=SPEXP, SPEXP2=SPEXP, ..., SPEXP#=SPEXP
C Other data to control plotting:
C     SPTEXT1='string'... Text to be written above the left-hand top
C             corner of the frame.
C             Default: SPTEXT1=' '
C     SPTEXT2='string'... Text to be written above the right-hand top
C             corner of the frame.
C             Default: SPTEXT2=' '
C     SPTEXT3='string'... Text to be written below the frame.
C             Default: SPTEXT3=' '
C     SPTEXT4='string'... Text to be written below text SPTEXT3.
C             Default: SPTEXT4=' '
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               receivers in file REC, see the amplitude scaling.
C     CALCOPS='string'... String with the PostScript instructions, see
C             file calcops.for.
C
C                                                     
C Input GSE files SS* 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 source name identified by
C     string 'NAMESRC', and for the hypocentral coordinates identified
C     by strings 'X1SRC', 'X2SRC' and 'X3SRC'.
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',X1SRC,X2SRC,X3SRC,/
C     'SRCNAME'... Name of the source.  Used only if the comment lines
C             of the GSE file contain the source name.
C             The source names cannot be longer than 6 characters.
C     X1SRC,X2SRC,X3SRC... Coordinates of the hypocentre.
C     Default: X1SRC=0., X2SRC=0., X3SRC=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 containing the name of the receiver.
C             The receiver names cannot be longer than 6 characters.
C             The limitation of the receiver names to 6 characters is
C             imposed by the GSE standard.
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 is plotted in colour K1.
C     Default: X1REC(IR)=0, X2REC(IR)=0, X3REC(IR)=0., K1=KOLORTD,
C             K2=KOLORTD
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
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C     Allocation of working arrays:
      INTEGER MPTS,MSS
      PARAMETER (MPTS=3000,MSS=MRAM-3*MPTS)
      CHARACTER*6 PTS(MPTS)
      COMMON/PTSC/PTS
C
C-----------------------------------------------------------------------
C
C     External functions and subroutines:
      EXTERNAL LENGTH
      INTEGER  LENGTH
      EXTERNAL STRIND
C     The length of character function STRIND is declared later on.
C
C     Input and output data filenames and logical unit numbers:
      INTEGER LU,LUPAR,MFILSS,MDIGSS
      PARAMETER (LU=1,LUPAR=2,MFILSS=29,MDIGSS=2)
C     MFILSS..Maximum number of input files with synthetic seismograms.
C     MDIGSS..Number of digits of MFILSS.
      CHARACTER*80 FILSEP,FILPAR,FILPTS,FILSRC,FILREC,FILFTT,FILHIL
      CHARACTER*80 FILOLD(0:MFILSS),FILESS(0:MFILSS),FILEPS(3)
C
C     Storing seismograms in memory
      INTEGER IFILO(0:MFILSS),IFILE(0:MFILSS),ISSRAM(0:MFILSS+1)
C     IFILO(I:I)... Index of the data corresponding to file FILOLD(I:I).
C     IFILE(I:I)... Index of the data corresponding to file FILESS(I:I).
C     ISSRAM(IFILE:IFILE)... End of IFILEth data stored in RAM.
C
C     Parameters and small working arrays:
      REAL UNDEF
      PARAMETER (UNDEF=-999999.)
      INTEGER KOLOR(0:MFILSS),KOMP(0:MFILSS,3)
      INTEGER ISS,ISP
      REAL SPXMIN(0:MFILSS),SPXMAX(0:MFILSS)
      REAL SPYMIN(0:MFILSS),SPYMAX(0:MFILSS)
      REAL SPAMP(0:MFILSS,3),SPOWER(0:MFILSS),SPEXP(0:MFILSS)
      REAL SPAMPX(3),XPTS(4),YPTS(4)
C
      CHARACTER*(6+MDIGSS) STRIND
      CHARACTER*5  STRKOM
      CHARACTER*1  CHAN,TEXT
      CHARACTER*6  NAMSRC,NAMREC
      CHARACTER*80 NAMPTS
C
C     Line of optional SEP parameter file or comments in the GSE file
      CHARACTER*80 LINE
C
C     Data specifying labels in the plot:
      CHARACTER*80 TEXT1,TEXT2,TEXT3,TEXT4
      CHARACTER*20 KXTEXT(5)
C
C     Lists of point coordinates, sources and receivers:
C     INTEGER NPTS,NSRC,NREC
C
      INTEGER MFILS1
      PARAMETER (MFILS1=MFILSS+1)
      DATA FILESS/MFILS1*' '/
      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
      IF (FILSEP.EQ.' ') THEN
C       SP-01
        CALL ERROR('SP-01: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
C     Reading all data from the SEP file into the memory:
      CALL RSEP1(LU,FILSEP)
      WRITE(*,'(A)') '+SP: Working...            '
C
C     Loop over SP executions in optional SP parameter file:
      CALL RSEP3T('SPPAR',FILPAR,' ')
      IF(FILPAR.NE.' ') THEN
        OPEN(LUPAR,FILE=FILPAR,STATUS='OLD')
        NSSRAM=0
        ISSRAM(0)=0
      END IF
  100 CONTINUE
      IF(FILPAR.NE.' ') THEN
C       Loop over lines the SP parameter file
  110   CONTINUE
          READ(LUPAR,'(A)',END=999) LINE
          CALL RSEP2(LINE)
          I=INDEX(LINE,'#')
          IF(I.EQ.0) THEN
            I=LENGTH(LINE)
          END IF
          I=INDEX(LINE(1:I),'sp:')
          IF(I.GT.1) THEN
            IF(LINE(I-1:I-1).NE.' ') THEN
              I=0
            END IF
          END IF
        IF(I.EQ.0) GO TO 110
C       SP execution is prescribed at the current positions in the SP
C       parameter file.
C       Copying filenames corresponding to the data stored in RAM:
        DO 111 I2=0,MFILSS
          FILOLD(I2)=FILESS(I2)
          IFILO (I2)=IFILE (I2)
  111   CONTINUE
      END IF
C
C     Input and output filenames:
      CALL RSEP3T('SS' ,FILESS(0),'ss.gse')
      DO 112 I1=1,MFILSS
        CALL RSEP3T(STRIND('SS',I1),FILESS(I1),' ')
  112 CONTINUE
      CALL RSEP3T('SP1',FILEPS(1),'ss1.ps')
      CALL RSEP3T('SP2',FILEPS(2),'ss2.ps')
      CALL RSEP3T('SP3',FILEPS(3),'ss3.ps')
      IF(FILPAR.EQ.' ') THEN
        ISS0=0
      ELSE
        DO 129 I2=0,MFILSS
          IF(FILOLD(I2).NE.' ') THEN
            DO 121 I1=0,MFILSS
              IF(FILESS(I1).EQ.FILOLD(I2)) THEN
                GO TO 128
              END IF
  121       CONTINUE
C             Removing seismograms from the memory
              I=ISSRAM(IFILE(I2))-ISSRAM(IFILE(I2)-1)
              DO 122 I1=ISSRAM(IFILE(I2))+1,ISSRAM(NSSRAM)
                RAM(I1-I)=RAM(I1)
  122         CONTINUE
              DO 123 I1=IFILE(I2),NSSRAM
                ISSRAM(I1-1)=ISSRAM(I1)-I
  123         CONTINUE
              NSSRAM=NSSRAM-1
              DO 124 I1=I2+1,MFILSS
                IF(FILOLD(I1).EQ.FILOLD(I2)) THEN
                  FILOLD(I1)=' '
                END IF
  124         CONTINUE
              FILOLD(I2)=' '
  128       CONTINUE
          END IF
  129   CONTINUE
        DO 139 I2=0,MFILSS
          IF(FILESS(I2).NE.' ') THEN
            DO 131 I1=0,MFILSS
              IF(FILESS(I2).EQ.FILOLD(I1)) THEN
                IFILE(I2)=IFILO(I1)
                GO TO 138
              END IF
  131       CONTINUE
            DO 132 I1=0,I2-1
              IF(FILESS(I2).EQ.FILESS(I1)) THEN
                IFILE(I2)=IFILE(I1)
                GO TO 138
              END IF
  132       CONTINUE
C             Reading seismograms into the memory
              WRITE(*,'(2A)') '+SP: Reading ',FILESS(I2)(1:66)
              OPEN(LU,FILE=FILESS(I2),STATUS='OLD')
              CALL RGSE1(LU,TEXT)
              ISS0=ISSRAM(NSSRAM)+1
  133         CONTINUE
                CALL RGSE2(LU,NAMREC,CHAN,I,X1R,X2R,X3R,T0,TD,
     *                                      NSS,MSS-ISS0-22,RAM(ISS0+1))
                IF(NSS.LE.-1) THEN
C                 End of the GSE file
                  GO TO 137
                END IF
  134           CONTINUE
                  CALL RGSE2C(LINE,*135)
                  CALL RSEP2(LINE)
                GO TO 134
  135           CONTINUE
                CALL RSEP3T('NAMESRC',NAMSRC,' ')
                CALL RSEP3R('X1SRC',X1S,0.)
                CALL RSEP3R('X2SRC',X2S,0.)
                CALL RSEP3R('X3SRC',X3S,0.)
                CALL WRAM2(ISS0,NAMSRC,X1S,X2S,X3S,
     *                          NAMREC,X1R,X2R,X3R,I,T0,TD,NSS,IRAM,RAM)
              GO TO 133
  137         CONTINUE
              NSSRAM=NSSRAM+1
              ISSRAM(NSSRAM)=ISS0
              IFILE(I2)=NSSRAM
              CLOSE(LU)
  138       CONTINUE
          END IF
  139   CONTINUE
      END IF
C
C     Reading lists of point coordinates, sources and receivers:
C     Point coordinates
      NPTS=0
      CALL RSEP3T('PTS',FILPTS,' ')
      IF(FILPTS.NE.' ') THEN
        OPEN(LU,FILE=FILPTS,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        DO 11 I=1,MPTS
          I0=MSS+3*I
          NAMPTS='$$$$$$$'
          RAM(I0-2)=0.
          RAM(I0-1)=0.
          RAM(I0  )=0.
          READ(LU,*,END=12) NAMPTS,RAM(I0-2),RAM(I0-1),RAM(I0)
          IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN
            GO TO 12
          END IF
          IF(NAMPTS(7:).NE.' ') THEN
C           SP-13
            CALL ERROR('SP-13: Point name exceeds 6 characters')
C           Names of points in file given by input parameter PTS cannot
C           be longer than 6 characters.  This limitation is imposed by
C           the GSE standard.
          END IF
          PTS(I)=NAMPTS(1:6)
   11   CONTINUE
C         SP-02
          CALL ERROR('SP-02: Array dimension MPTS small for points')
   12   CONTINUE
        NPTS=I-1
        CLOSE(LU)
      END IF
C     Sources
      NSRC=0
      CALL RSEP3T('SRC',FILSRC,' ')
      IF(FILSRC.NE.' ') THEN
        OPEN(LU,FILE=FILSRC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        DO 13 I=NPTS+1,MPTS
          I0=MSS+3*I
          NAMPTS='$$$$$$$'
          RAM(I0-2)=0.
          RAM(I0-1)=0.
          RAM(I0  )=0.
          READ(LU,*,END=14) NAMPTS,RAM(I0-2),RAM(I0-1),RAM(I0)
          IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN
            GO TO 14
          END IF
          IF(NAMPTS(7:).NE.' ') THEN
C           SP-14
            CALL ERROR('SP-14: Source name exceeds 6 characters')
C           Names of sources in file given by input parameter SRC cannot
C           be longer than 6 characters.  This is the limitation imposed
C           by the GSE standard on the receiver names and applied here
C           also to the source names.
          END IF
          PTS(I)=NAMPTS(1:6)
   13   CONTINUE
C         SP-03
          CALL ERROR('SP-03: Array dimension MPTS small for sources')
   14   CONTINUE
        NSRC=I-NPTS-1
        CLOSE(LU)
      END IF
C     Receivers
      NREC=0
      RECNUM=999.
      CALL RSEP3T('REC',FILREC,' ')
      IF(FILREC.NE.' ') THEN
        OPEN(LU,FILE=FILREC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        DO 15 I=NPTS+NSRC+1,MPTS
          I0=MSS+3*I
          NAMPTS='$$$$$$$'
          RAM(I0-2)=0.
          RAM(I0-1)=0.
          RAM(I0  )=0.
          READ(LU,*,END=16) NAMPTS,RAM(I0-2),RAM(I0-1),RAM(I0)
          IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN
            GO TO 16
          END IF
          IF(NAMPTS(7:).NE.' ') THEN
C           SP-15
            CALL ERROR('SP-15: Receiver name exceeds 6 characters')
C           Names of receivers in file given by input parameter REC
C           cannot be longer than 6 characters.  This limitation is
C           imposed by the GSE standard.
          END IF
          PTS(I)=NAMPTS(1:6)
   15   CONTINUE
C         SP-04
          CALL ERROR('SP-04: Array dimension MPTS small for points')
   16   CONTINUE
        NREC=I-NPTS-NSRC-1
        RECNUM=FLOAT(NREC)
        CLOSE(LU)
      END IF
C
C     Components:
      DO 152 I2=1,3
        STRKOM=STRIND('KOMP',I2)
        CALL RSEP3I(STRKOM,KOMP(0,I2),I2)
        DO 151 I1=1,MFILSS
          CALL RSEP3I(STRIND(STRKOM,I1),KOMP(I1,I2),KOMP(0,I2))
  151   CONTINUE
  152 CONTINUE
C
C     Colours:
      CALL RSEP3I('KOLOR ',KOLOR(0),1)
      DO 153 I1=1,MFILSS
        CALL RSEP3I(STRIND('KOLOR',I1),KOLOR(I1),I1+1)
  153 CONTINUE
      CALL RSEP3I('KOLORTT',KOLORT,1)
      CALL RSEP3I('KOLORTD',KOLORD,0)
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('SPXLEN',SPXLEN,RECNUM+1.)
      CALL RSEP3R('SPXDIV',SPXDIV, 1.)
      CALL RSEP3R('SPXSUB',SPXSUB, 1.)
      CALL RSEP3R('SPXMIN ',SPXMIN(0), 0.)
      DO 154 I1=1,MFILSS
        CALL RSEP3R(STRIND('SPXMIN',I1),SPXMIN(I1),SPXMIN(0))
  154 CONTINUE
      CALL RSEP3R('SPXMAX ',SPXMAX(0),RECNUM+1.)
      DO 155 I1=1,MFILSS
        CALL RSEP3R(STRIND('SPXMAX',I1),SPXMAX(I1),SPXMAX(0))
  155 CONTINUE
      CALL RSEP3R('SPYMIN ',SPYMIN(0), 0.)
      DO 156 I1=1,MFILSS
        CALL RSEP3R(STRIND('SPYMIN',I1),SPYMIN(I1),SPYMIN(0))
  156 CONTINUE
      CALL RSEP3R('SPYMAX ',SPYMAX(0), 0.)
      DO 157 I1=1,MFILSS
        CALL RSEP3R(STRIND('SPYMAX',I1),SPYMAX(I1),SPYMAX(0))
  157 CONTINUE
C     Characters:
      CALL RSEP3T('SPTEXT1',TEXT1,' ')
      CALL RSEP3T('SPTEXT2',TEXT2,' ')
      CALL RSEP3T('SPTEXT3',TEXT3,' ')
      CALL RSEP3T('SPTEXT4',TEXT4,' ')
      CALL RSEP3R('SPCHRH',SPCHRH, 0.4)
C
C     Amplitude scaling:
      CALL RSEP3I('NORMSP',NORMSP,0)
      CALL RSEP3R('SPAMP ',SPAMP(0,1),1.)
      DO 161 I1=1,MFILSS
        CALL RSEP3R(STRIND('SPAMP',I1),SPAMP(I1,1),SPAMP(0,1))
  161 CONTINUE
      CALL RSEP3R('SPAMPX1',SPAMPX(1),1.)
      CALL RSEP3R('SPAMPX2',SPAMPX(2),1.)
      CALL RSEP3R('SPAMPX3',SPAMPX(3),1.)
      DO 168 I2=3,1,-1
        DO 167 I1=0,MFILSS
          SPAMP(I1,I2)=SPAMP(I1,1)*SPAMPX(I2)
  167   CONTINUE
  168 CONTINUE
      CALL RSEP3R('SPDIST' ,SPDIST,1.)
      CALL RSEP3R('SPOWER' ,SPOWER(0),0.)
      DO 164 I1=1,MFILSS
        CALL RSEP3R(STRIND('SPOWER',I1),SPOWER(I1),SPOWER(0))
  164 CONTINUE
      CALL RSEP3R('SPEXP ',SPEXP(0),0.)
      DO 165 I1=1,MFILSS
        CALL RSEP3R(STRIND('SPEXP',I1),SPEXP(I1),SPEXP(0))
  165 CONTINUE
      CALL RSEP3R('SPEXPT',SPEXPT,0.)
C
      SSTMIN=AMIN1(SPTMIN,SPTMAX)
      SSTMAX=AMAX1(SPTMIN,SPTMAX)
      SCY= SPTLEN/(SPTMAX-SPTMIN)
      XA = SPXMAX(0)-SPXMIN(0)
      YA = SPYMAX(0)-SPYMIN(0)
C
C     Higlighting given areas (e.g., travel times with error bars):
      CALL RSEP3T('FTT'   ,FILFTT,' ')
      CALL RSEP3T('SPHILI',FILHIL,' ')
      CALL RSEP3R('SPHIWI',SPHIWI,SPXLEN/(RECNUM+1.))
C
C.......................................................................
C
C     Loop over 3 seismogram files:
      DO 99 ISP=1,3
        IF(FILEPS(ISP).NE.' ') THEN
          WRITE(*,'(2A)') '+SP: Plotting ',FILEPS(ISP)(1:65)
C
C         Initialization of CALCOMP:
          CALL PLOTN(FILEPS(ISP),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
          CALL FRAME(SPXLEN,SPTLEN,ABS(SPXDIV),SPXSUB,
     *                             ABS(SPTDIV),SPTSUB,0,IX,IT,SPCHRH,
     *               SPXMIN(0),SPXMAX(0),KXTEXT(KODESP+1),
     *               SPYMIN(0),SPYMAX(0),'X2',
     *               SPTMIN,SPTMAX,'TIME',
     *               TEXT1,0,0.,0,TEXT2,0,0.,0,TEXT3,TEXT4)
C
C         Higlighting travel times of file FTT:
          IF(FILFTT.NE.' ') THEN
            OPEN(LU,FILE=FILFTT,STATUS='OLD')
            READ(LU,*) (TEXT,I=1,20)
C           Loop for areas to highlight:
   20       CONTINUE
              NAMSRC='$$$$$$'
              T0=0.
              TD=0.
              READ(LU,*,END=29) NAMSRC,NAMREC,T0,TD
              IF(NAMSRC.EQ.'$$$$$$') THEN
C               End of travel-time file
                GO TO 29
              END IF
C             Selecting the receiver:
              IF(FILREC.NE.' ') THEN
C               Loop for receivers
                DO 21 I=NPTS+NSRC+1,NPTS+NSRC+NREC
                  IF(PTS(I).EQ.NAMREC) THEN
                    IREC=I-NPTS-NSRC
                    GO TO 22
                  END IF
   21           CONTINUE
                GO TO 20
              END IF
   22         CONTINUE
C             Selecting the source:
              IF(FILSRC.NE.' ') THEN
C               Loop for sources
                DO 23 I=NPTS+1,NPTS+NSRC
                  IF(PTS(I).EQ.NAMSRC) GO TO 24
   23           CONTINUE
                GO TO 20
              END IF
   24         CONTINUE
C             Finding the receiver coordinates:
              DO 25 I=1,NPTS
                IF(PTS(I).EQ.NAMREC) THEN
                  I0=MSS+3*I
                  X1R=RAM(I0-2)
                  X2R=RAM(I0-1)
                  X3R=RAM(I0)
                  GO TO 26
                END IF
   25         CONTINUE
C               SP-05
                LINE='SP-05: Receiver '//NAMREC(1:LENGTH(NAMREC))
     *                                 //' not found in file PTS'
                CALL ERROR(LINE)
C               If file FTT with travel times is given, file PTS has
C               to contain all receiver names of file REC (if REC is
C               specified) or all receiver names of file FTT (if
C               REC=' ').
   26         CONTINUE
C             Finding the source coordinates:
              DO 27 I=1,NPTS
                IF(PTS(I).EQ.NAMSRC) THEN
                  I0=MSS+3*I
                  X1S=RAM(I0-2)
                  X2S=RAM(I0-1)
                  X3S=RAM(I0)
                  GO TO 28
                END IF
   27         CONTINUE
C               SP-06
                LINE='SP-06: Source '//NAMSRC(1:LENGTH(NAMSRC))
     *                               //' not found in file PTS'
                CALL ERROR(LINE)
C               If file FTT with travel times is given, file PTS has
C               to contain all receiver names of file REC (if REC is
C               specified) or all receiver names of file FTT (if
C               REC=' ').
   28         CONTINUE
C             Reduction of the travel time
              IF(SPVRED.NE.0.) THEN
                RR=SQRT((X1R-X1S)**2+(X2R-X2S)**2+(X3R-X3S)**2)
                T0=T0-RR/SPVRED
              END IF
              IF(KODESP.EQ.0) THEN
                IF(FILREC.EQ.' ') THEN
C                 SP-07
                  CALL ERROR('SP-07: No receiver file specified')
C                 For KODESP=0, filename REC must be specified in the
C                 input data.
                END IF
                X=SPXLEN*(FLOAT(IREC)-SPXMIN(0))
     *                  /(SPXMAX(0)-SPXMIN(0))
              ELSE IF(KODESP.EQ.1.OR.KODESP.EQ.2) THEN
                X=SPXLEN*((X1R-SPXMIN(0))*XA+(X2R-SPYMIN(0))*YA)
     *                  /(XA*XA+YA*YA)
              ELSE IF(KODESP.EQ.3) THEN
                X=SPXLEN*(X3R-SPXMIN(0))/XA
              ELSE
                X=SPXLEN*(RR-SPXMIN(0))/XA
              END IF
C             Plotting the highlight:
              IF(KOLORD.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(KOLORD.GT.0) THEN
                  CALL NEWPEN(KOLORD)
                  CALL FILL(XPTS,YPTS,4)
                ELSE
                  CALL NEWPEN(KOLORT)
                  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)
                  CALL PLOT(XPTS(1),YPTS(1),2)
                END IF
              END IF
              IF(KOLORT.GT.0) THEN
                CALL NEWPEN(KOLORT)
                CALL PLOT(X-0.5*SPHIWI,SCY*(T0-SPTMIN),3)
                CALL PLOT(X+0.5*SPHIWI,SCY*(T0-SPTMIN),2)
              END IF
            GO TO 20
   29       CONTINUE
C           Closing highlighting file
            CLOSE(LU)
          END IF
C
C         Higlighting times given in file SPHILI:
          IF(FILHIL.NE.' ') THEN
            OPEN(LU,FILE=FILHIL,STATUS='OLD')
            READ(LU,*) (TEXT,I=1,20)
C           Loop for areas to highlight:
   30       CONTINUE
              NAMREC='$$$$$$'
              T0=UNDEF
              TD=UNDEF
              K1=KOLORT
              K2=KOLORD
              READ(LU,*,END=39) NAMREC,X1R,X2R,X3R,T0,TD,K1,K2
              IF(NAMREC.EQ.'$$$$$$') THEN
                GO TO 39
              END IF
C             Selecting the receiver:
              IF(FILREC.NE.' ') THEN
C               Loop for receivers
                DO 31 I=NPTS+NSRC+1,NPTS+NSRC+NREC
                  IF(PTS(I).EQ.NAMREC) THEN
                    IREC=I-NPTS-NSRC
                    GO TO 32
                  END IF
   31           CONTINUE
                GO TO 30
              END IF
   32         CONTINUE
C             Updating the coordinates:
              IF(FILPTS.NE.' '.AND.NAMREC.NE.' ') THEN
C               Receiver coordinates
                DO 33 I=1,NPTS
                  IF(PTS(I).EQ.NAMREC) THEN
                    I0=MSS+3*I
                    X1R=RAM(I0-2)
                    X2R=RAM(I0-1)
                    X3R=RAM(I0)
                    GO TO 34
                  END IF
   33           CONTINUE
C                 SP-08
                  LINE='SP-08: Receiver '//NAMREC(1:LENGTH(NAMREC))
     *                                   //' not found in file PTS'
                  CALL ERROR(LINE)
C                 If file PTS with the coordinates of points is given
C                 and file SPHILI contains receiver names, file PTS has
C                 to contain all receiver names of file REC (if REC is
C                 specified) or all receiver names of file SPHILI (if
C                 REC=' ').
   34           CONTINUE
              END IF
              IF(SPVRED.NE.0.) THEN
                IF(FILSRC.EQ.' ') THEN
                  X1S=0.
                  X2S=0.
                  X3S=0.
                ELSE
                  I0=MSS+3*NPTS
                  X1S=RAM(I0+1)
                  X2S=RAM(I0+2)
                  X3S=RAM(I0+3)
                END IF
                RR=SQRT((X1R-X1S)**2+(X2R-X2S)**2+(X3R-X3S)**2)
                T0=T0-RR/SPVRED
              END IF
              IF(KODESP.EQ.0) THEN
                IF(FILREC.EQ.' ') THEN
C                 SP-09
                  CALL ERROR('SP-09: No receiver file specified')
C                 For KODESP=0, filename REC must be specified in the
C                 input data.
                END IF
                X=SPXLEN*(FLOAT(IREC)-SPXMIN(0))
     *                  /(SPXMAX(0)-SPXMIN(0))
              ELSE IF(KODESP.EQ.1.OR.KODESP.EQ.2) THEN
                X=SPXLEN*((X1R-SPXMIN(0))*XA+(X2R-SPYMIN(0))*YA)
     *                  /(XA*XA+YA*YA)
              ELSE IF(KODESP.EQ.3) THEN
                X=SPXLEN*(X3R-SPXMIN(0))/XA
              ELSE
                X=SPXLEN*(RR-SPXMIN(0))/XA
              END IF
C             Plotting the highlight:
              IF(T0.EQ.UNDEF) THEN
                K1=-1
                K2=-1
              END IF
              IF(TD.EQ.UNDEF) THEN
                K2=-1
              END IF
              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(K1)
                  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 30
   39       CONTINUE
C           Closing highlighting file
            CLOSE(LU)
          END IF
C
C         Plotting seismograms:
C         Loop for GSE files
          DO 90 ISS=0,MFILSS
            IF(FILESS(ISS).NE.' ') THEN
              CALL NEWPEN(KOLOR(ISS))
              XA=SPXMAX(ISS)-SPXMIN(ISS)
              YA=SPYMAX(ISS)-SPYMIN(ISS)
C             Opening input GSE file with the seismograms:
              IF(FILPAR.EQ.' ') THEN
                OPEN(LU,FILE=FILESS(ISS),STATUS='OLD')
                CALL RGSE1(LU,TEXT)
              ELSE
                ISS0=ISSRAM(IFILE(ISS)-1)+1
              END IF
C             Loop for seismograms
              IREC=0
   40         CONTINUE
   41           CONTINUE
C               Selecting the component:
   42           CONTINUE
                  ISS1=ISS0
                  IF(FILPAR.EQ.' ') THEN
                    CALL RGSE2
     *                  (LU,NAMREC,CHAN,I,X1R,X2R,X3R,T0,TD,NSS,MSS,RAM)
                  ELSE
                    CALL RRAM2(ISS0,NAMSRC,X1S,X2S,X3S,
     *                          NAMREC,X1R,X2R,X3R,I,T0,TD,NSS,IRAM,RAM)
                  END IF
                  IF(NSS.LE.-1) THEN
C                   End of the GSE file
                    GO TO 80
                  END IF
                IF(I.NE.KOMP(ISS,ISP)) GO TO 42
C               Selecting the receiver:
                IF(FILREC.EQ.' ') THEN
                  IREC=IREC+1
                ELSE
C                 Loop for receivers
                  DO 51 I=NPTS+NSRC+1,NPTS+NSRC+NREC
                    IF(PTS(I).EQ.NAMREC) THEN
                      I0=MSS+3*I
                      X1R=RAM(I0-2)
                      X2R=RAM(I0-1)
                      X3R=RAM(I0)
                      IREC=I-NPTS-NSRC
                      GO TO 52
                    END IF
   51             CONTINUE
                  GO TO 41
                END IF
   52           CONTINUE
C               Reading the source information:
                IF(FILPAR.EQ.' ') THEN
   53             CONTINUE
                    CALL RGSE2C(LINE,*54)
                    CALL RSEP2(LINE)
                  GO TO 53
   54             CONTINUE
                  CALL RSEP3T('NAMESRC',NAMSRC,' ')
                  CALL RSEP3R('X1SRC',X1S,0.)
                  CALL RSEP3R('X2SRC',X2S,0.)
                  CALL RSEP3R('X3SRC',X3S,0.)
                END IF
C               Selecting the source:
                IF(FILSRC.NE.' ') THEN
                  IF(NAMSRC.EQ.' ') THEN
                    I0=MSS+3*NPTS
                    X1S=RAM(I0+1)
                    X2S=RAM(I0+2)
                    X3S=RAM(I0+3)
                  ELSE
C                   Loop for sources
                    DO 55 I=NPTS+1,NPTS+NSRC
                      IF(PTS(I).EQ.NAMSRC) THEN
                        I0=MSS+3*I
                        X1S=RAM(I0-2)
                        X2S=RAM(I0-1)
                        X3S=RAM(I0)
                        GO TO 56
                      END IF
   55               CONTINUE
                    GO TO 41
                  END IF
                END IF
   56           CONTINUE
C
C               Updating the coordinates:
                IF(FILPTS.NE.' ') THEN
C                 Receiver coordinates
                  DO 61 I=1,NPTS
                    IF(PTS(I).EQ.NAMREC) THEN
                      I0=MSS+3*I
                      X1R=RAM(I0-2)
                      X2R=RAM(I0-1)
                      X3R=RAM(I0)
                      GO TO 62
                    END IF
   61             CONTINUE
C                   SP-10
                    LINE='SP-10: Receiver '//NAMREC(1:LENGTH(NAMREC))
     *                                     //' not found in file PTS'
                    CALL ERROR(LINE)
C                   If file PTS with the coordinates of points is given,
C                   it has to contain all receiver names of file REC
C                   (if REC is specified) or all receiver names of the
C                   GSE file (if REC=' ').
   62             CONTINUE
C                 Source coordinates
                  IF(NAMSRC.NE.' ') THEN
                    DO 63 I=1,NPTS
                      IF(PTS(I).EQ.NAMSRC) THEN
                        I0=MSS+3*I
                        X1S=RAM(I0-2)
                        X2S=RAM(I0-1)
                        X3S=RAM(I0)
                        GO TO 64
                      END IF
   63               CONTINUE
C                   SP-11
                    CALL ERROR('SP-11: Source not found in file PTS')
C                   If file PTS with the coordinates of points is given
C                   and the GSE file contains source names, file PTS has
C                   to contain all source names of file SRC (if SRC is
C                   specified) or all source names of the GSE file (if
C                   SRC=' ').
   64               CONTINUE
                  END IF
                END IF
C
C               Plotting the seismogram:
                RR=SQRT((X1R-X1S)**2+(X2R-X2S)**2+(X3R-X3S)**2)
                IF(SPEXP(ISS).NE.0.) THEN
                  DO 76 I=1,NSS
                    RAM(ISS1+I)=RAM(ISS1+I)
     *                        *EXP(SPEXP(ISS)*(T0+TD*FLOAT(I-1)-SPEXPT))
   76             CONTINUE
                END IF
                IF(SPVRED.NE.0.) THEN
                  T0=T0-RR/SPVRED
                END IF
C
                IF(NORMSP.LE.0) THEN
                  IF(NORMSP.LT.0) THEN
                    I1=MAX0(INT((SSTMIN-T0)/TD+2.),1)
                    I2=MIN0(INT((SSTMAX-T0)/TD+1.),NSS)
                  ELSE
                    I1=1
                    I2=NSS
                  END IF
                  AA=0.
                  DO 77 I=I1,I2
                    AA=AMAX1(ABS(RAM(ISS1+I)),AA)
   77             CONTINUE
                  IF(AA.NE.0.) THEN
                    AA=SPAMP(ISS,ISP)*SPXLEN/(RECNUM+1.)/AA
                  END IF
                ELSE
                  IF(SPOWER(ISS).EQ.0.) THEN
                    AA=SPAMP(ISS,ISP)
                  ELSE
                    AA=SPAMP(ISS,ISP)*((RR/SPDIST)**SPOWER(ISS))
                  END IF
                END IF
C
                IF(KODESP.LE.0) THEN
                  IF(FILREC.EQ.' ') THEN
C                   SP-12
                    CALL ERROR('SP-12: No receiver file specified')
C                   For KODESP=0, filename REC must be specified in the
C                   input data.
                  END IF
                  X=SPXLEN*(FLOAT(IREC)-SPXMIN(ISS))
     *                    /(SPXMAX(ISS)-SPXMIN(ISS))
                ELSE IF(KODESP.EQ.1.OR.KODESP.EQ.2) THEN
                  X=SPXLEN*((X1R-SPXMIN(ISS))*XA+(X2R-SPYMIN(ISS))*YA)
     *                    /(XA*XA+YA*YA)
                ELSE IF(KODESP.EQ.3) THEN
                  X=SPXLEN*(X3R-SPXMIN(ISS))/XA
                ELSE
                  X=SPXLEN*(RR-SPXMIN(ISS))/XA
                END IF
                Y=SCY*(SSTMIN-SPTMIN)
C
C               Denoting the horizontal axis by the receiver name:
C               IF(ISS.EQ.0) THEN
                  IF(KODESP.LE.0.AND.SPXDIV.GT.0.) THEN
                    CALL SYMBOL(X-0.5*(LENGTH(NAMREC)-0.43)*SPCHRH,
     *                     Y-1.8*SPCHRH,SPCHRH,NAMREC,0.,LENGTH(NAMREC))
                  END IF
C               END IF
C
                CALL PLOT(X,Y,3)
                DO 78 I=1,NSS
                  T = T0+TD*FLOAT(I-1)
                  IF (SSTMIN.LE.T.AND.T.LE.SSTMAX) THEN
                    IF (T.LT.SSTMIN+TD) THEN
C                     Bottom intersection point of the seismogram
                      A=(T-SSTMIN)/TD
                      IF(I.GT.1) THEN
                        A=(1.-A)*RAM(ISS1+I)+A*RAM(ISS1+I-1)
                      ELSE
                        A=(1.-A)*RAM(ISS1+I)
                      END IF
                      A=X-AA*A
                      Y=SCY*(SSTMIN-SPTMIN)
                      CALL PLOT(A,Y,2)
                    ELSE IF (I.EQ.1) THEN
C                     Straight line before the seismogram
                      Y = SCY*(T-TD-SPTMIN)
                      CALL PLOT(X,Y,2)
                    END IF
                    A = X-AA*RAM(ISS1+I)
                    Y = SCY*(T-SPTMIN)
                    CALL PLOT(A,Y,2)
                    IF (T.GT.SSTMAX-TD) THEN
C                     Top intersection point of the seismogram
                      A=(SSTMAX-T)/TD
                      IF(I.LT.NSS) THEN
                        A=(1.-A)*RAM(ISS1+I)+A*RAM(ISS1+I+1)
                      ELSE
                        A=(1.-A)*RAM(ISS1+I)
                      END IF
                      A=X-AA*A
                      Y=SCY*(SSTMAX-SPTMIN)
                      CALL PLOT(A,Y,2)
                    ELSE IF (I.EQ.NSS) THEN
C                     Straight line after the seismogram
                      Y = SCY*(T+TD-SPTMIN)
                      CALL PLOT(X,Y,2)
                    END IF
                  END IF
   78           CONTINUE
                Y=SCY*(SSTMAX-SPTMIN)
                CALL PLOT(X,Y,2)
C
              GO TO 40
   80         CONTINUE
C             End of loop for receivers
C             Closing GSE file
              IF(FILPAR.EQ.' ') THEN
                CLOSE(LU)
              END IF
            END IF
   90     CONTINUE
C         End of loop for GSE files
          CALL PLOT(0.,0.,999)
        END IF
   99 CONTINUE
C
C     End of loop over SP executions in optional SP parameter file:
      IF(FILPAR.NE.' ') THEN
        GO TO 100
      END IF
  999 CONTINUE
      IF(FILPAR.NE.' ') THEN
        CLOSE(LUPAR)
      END IF
C
      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,KX1,XL2,XR2,KX2,YL,YR,KY,
     *                 K1,M1,A1,N1,K2,M2,A2,N2,K3,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,M1,N1,M2,N2
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     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     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     KY...   Label of the vertical axis.  String to be written to the
C             left of the vertical axis.
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     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     K3...   String to be written below the frame, starting 5.3*HEIGHT
C             below the left bottom corner.
C     K4...   String to be written below the frame, starting 7.0*HEIGHT
C             below the left bottom corner.
C
C-----------------------------------------------------------------------
C
      EXTERNAL LENGTH
      INTEGER  LENGTH
C
      INTEGER LX1,LX2,LY,L1,L2,L3,L4
C
C     Lengths of input strings
      LX1=LENGTH(KX1)
      LX2=LENGTH(KX2)
      LY =LENGTH(KY)
      L1 =LENGTH(K1)
      L2 =LENGTH(K2)
      L3 =LENGTH(K3)
      L4 =LENGTH(K4)
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
      SUBROUTINE WRAM2(ISS0,NAMSRC,X1S,X2S,X3S,
     *                      NAMREC,X1R,X2R,X3R,KOMP,T0,TD,NSS,IRAM,RAM)
      CHARACTER*(*) NAMSRC,NAMREC
      INTEGER IRAM(*)
      REAL RAM(*)
C
      CHARACTER*6 NAME
      IRAM(ISS0)=22+NSS
      I=ISS0+IRAM(ISS0)
      NAME=NAMSRC
      IRAM(I-21)=ICHAR(NAME(1:1))
      IRAM(I-20)=ICHAR(NAME(2:2))
      IRAM(I-19)=ICHAR(NAME(3:3))
      IRAM(I-18)=ICHAR(NAME(4:4))
      IRAM(I-17)=ICHAR(NAME(5:5))
      IRAM(I-16)=ICHAR(NAME(6:6))
      NAME=NAMREC
      IRAM(I-15)=ICHAR(NAME(1:1))
      IRAM(I-14)=ICHAR(NAME(2:2))
      IRAM(I-13)=ICHAR(NAME(3:3))
      IRAM(I-12)=ICHAR(NAME(4:4))
      IRAM(I-11)=ICHAR(NAME(5:5))
      IRAM(I-10)=ICHAR(NAME(6:6))
      RAM(I-9)=X1S
      RAM(I-8)=X2S
      RAM(I-7)=X3S
      RAM(I-6)=X1R
      RAM(I-5)=X2R
      RAM(I-4)=X3R
      IRAM(I-3)=KOMP
      RAM(I-2)=T0
      RAM(I-1)=TD
      ISS0=I
      IRAM(ISS0)=0
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE RRAM2(ISS0,NAMSRC,X1S,X2S,X3S,
     *                      NAMREC,X1R,X2R,X3R,KOMP,T0,TD,NSS,IRAM,RAM)
      CHARACTER*(*) NAMSRC,NAMREC
      INTEGER IRAM(*)
      REAL RAM(*)
C
      CHARACTER*6 NAME
      NSS=IRAM(ISS0)-22
      I=ISS0+IRAM(ISS0)
      IF(NSS.GT.-1) THEN
        NAME(1:1)=CHAR(IRAM(I-21))
        NAME(2:2)=CHAR(IRAM(I-20))
        NAME(3:3)=CHAR(IRAM(I-19))
        NAME(4:4)=CHAR(IRAM(I-18))
        NAME(5:5)=CHAR(IRAM(I-17))
        NAME(6:6)=CHAR(IRAM(I-16))
        NAMSRC=NAME
        NAME(1:1)=CHAR(IRAM(I-15))
        NAME(2:2)=CHAR(IRAM(I-14))
        NAME(3:3)=CHAR(IRAM(I-13))
        NAME(4:4)=CHAR(IRAM(I-12))
        NAME(5:5)=CHAR(IRAM(I-11))
        NAME(6:6)=CHAR(IRAM(I-10))
        NAMREC=NAME
        X1S=RAM(I-9)
        X2S=RAM(I-8)
        X3S=RAM(I-7)
        X1R=RAM(I-6)
        X2R=RAM(I-5)
        X3R=RAM(I-4)
        KOMP=IRAM(I-3)
        T0=RAM(I-2)
        TD=RAM(I-1)
      END IF
      ISS0=I
      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