C
C Program SP (Seismogram Plotting) to plot seismograms previously stored
C in the GSE data exchange format.
C
C Version: 5.10
C Date: 1997, September 26
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C 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 Description of data files:
C
C The data are read in by the list directed input (free format).
C In the description of data files, each numbered paragraph indicates
C the beginning of a new input operation (new READ statement).
C If the symbolic name of the input variable is enclosed in apostrophes,
C the corresponding value in input data is of the type CHARACTER, i.e.
C it should be a character string enclosed in apostrophes.  If the first
C letter of the symbolic name is I-N, the corresponding value is of the
C type INTEGER.  Otherwise, the input parameter is of the type REAL and
C may or may not contain a decimal point.
C
C Input data read from the * external unit:
C     The interactive * external unit may also be redirected to the file
C     containing the relevant data.
C (1) 'SPDAT','GSE',KOMP,'SRC',/
C     'SPDAT'... String with the name of the input data file containing
C             the data describing how to plot the seismograms (e.g.,
C             scale, minimum and maximum values of the coordinates,
C             etc.).
C             Description of file SPDAT
C     'GSE'... String with the name of the input data file in the GSE
C             data exchange format, containing the seismograms to be
C             plotted stored .  In the present version, all seismograms
C             must correspond to a single component.  Other components
C             must be stored in separate files.
C             Description of file GSE
C     KOMP... Component of the seismograms to be plotted.
C     'SRC'.. String with the name of the input data file containing the
C             coordinates of the hypocentre.  If left blank (default),
C             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: 'SPDAT'='sp.dat', 'GSE'='rf.out', KOMP=3, 'SRC'=' '.
C
C                                                   
C Input file SPDAT with the data to control plotting:
C (1) I1,J1,K1,I2,J2,K2,I3,J3,K3,I4,J4,K4,I5,J5,K5,I6,J6,K6,/
C     Selection of seismograms to be plotted.
C     First, the I1th to J1th seismograms of input file 'GSE', with
C     increment K1, are plotted.  Then I2th to J2th seismograms with
C     increment K2, and so on, until zero Ki is encountered.
C     The most frequent form of this input line will be
C     1,N,1,/
C     where N is the number of seismograms located in input file 'GSE'.
C (2) TMIN,TMAX,TAXIS,TM,TN,CODE,VRED,/
C     Time axis (vertical).
C     TMIN,TMAX... Minimum and maximum time (or reduced time), i.e.
C             the times corresponding to the bottom and top of the
C             seismogram plot.
C     TAXIS...Length of the vertical time axis in cm.
C     TM...   ABS(TM) is the number of intervals along the time axis,
C               starting at the bottom.  Must be TM.NE.0.
C             TM.GT.0.: the marks at the endpoints of intervals will
C               be supplemented with corresponding times (or reduced
C               times).
C             TM.LT.0.: the time axis will have no description.
C     TN...   Number of subintervals to be marked in each interval.
C             Must be TN.GT.0.
C     CODE... Specifies the distribution and description of seismograms
C             in the plot.  See the description of XMIN,XMAX,YMIN,YMAX
C             below.
C     VRED... Reduction velocity.  If non-zero, the time at each
C             receiver is reduced by the value of RR/VRED, where RR is
C             the hypocentral distance.
C             No time reduction is applied if VRED=0.
C     Default: TMIN=0., TMAX=1., TAXIS=10., TM=1., TN=1., CODE=0.,
C             VRED=0.
C (3) XMIN,XMAX,XAXIS,XM,XN,YMIN,YMAX,HEIGHT,/
C     Distance axis (horizontal).
C     XMIN,XMAX,YMIN,YMAX:
C             For CODE=0.: Horizontal axis is divided into MR+1
C               intervals, where MR 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, XMIN,XMAX,YMIN,YMAX have no meaning.
C             For CODE=1.: Horizontal axis represents the profile with
C               endpoints (X1,X2)=(XMIN,YMIN) and (X1,X2)=(XMAX,YMAX),
C               situated in a horizontal plane.  The seismograms are
C               located at the orthogonal projections of the receivers
C               onto the profile.  If XM.GT.0., the horizontal axis is
C               supplemented with the values of the X1 coordinates.
C             For CODE=2.: Horizontal axis represents the profile with
C               endpoints (X1,X2)=(XMIN,YMIN) and (X1,X2)=(XMAX,YMAX),
C               situated in a horizontal plane.  The seismograms are
C               located at the orthogonal projections of the receivers
C               onto the profile.  If XM.GT.0., the horizontal axis is
C               supplemented with the values of both X1 and X2
C               coordinates.
C             For CODE=3.: Horizontal axis represents the vertical
C               profile with endpoints X3=XMIN and X3=XMAX.
C               The seismograms are located at the horizontal
C               projections of the receivers onto the profile.
C               If XM.GT.0., the horizontal axis is supplemented with
C               the values of X3 coordinates.
C               YMIN and YMAX have no meaning.
C             For CODE=4.: Horizontal axis represents the hypocentral
C               distance with endpoints RR=XMIN and RR=XMAX.
C               The seismograms are located according to the hypocentral
C               distances the receivers.
C               If XM.GT.0., the horizontal axis is supplemented with
C               the values of the hypocentral distance.
C               YMIN and YMAX have no meaning.
C     XAXIS...Length of the horizontal axis in cm.
C     XM...   ABS(XM) is the number of intervals along the horizontal
C               axis, starting at the left.  Must be XM.NE.0.
C             XM.GT.0.: the marks at the endpoints of intervals will
C               be supplemented with the corresponding values.
C             XM.LT.0.: the horizontal axis will have no description.
C     XN...   Number of subintervals to be marked in each interval.
C             Must be XN.GT.0.
C     HEIGHT..Height of the characters in cm.
C     Default: XMIN=0., XMAX=1., XAXIS=FLOAT(MR+1), XM=1., XN=1.,
C             YMIN=0., YMAX=0., HEIGHT=0.4, where MR is the number of
C             seismograms to be plotted, see (1).
C (4) AMP,B1,EPIC,EPS,/
C     AMP...  AMP=0: Maximum amplitude at each trace normalized to the
C               fixed value.
C             Otherwise: All seismograms have the same scaling or power
C               scaling.
C     B1,EPIC,EPS... Amplitude scaling parameters.  EPIC and EPS are
C             not used if AMP=0.
C             Amplitude scale AA is for AMP=0:
C               AA=B1*XD/AMAX
C               where XD is the average distance between seismograms,
C               and AMAX is the maximum amplitude in each seismogram.
C             Otherwise:
C               AA=B1*(RR/EPIC)**EPS
C               where EPIC is the hypocentral distance of each receiver.
C     Default: AMP=0., B1 =1., EPIC=1., EPS=0.
C Example of data SPDAT
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'.  In the present version, all seismograms must
C     correspond to a single component.  Other components must be stored
C     in separate files.
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 '.
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 Example of data SRC
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     Filenames:
      CHARACTER*80 FSPDAT,FILERF,FILSRC
C
      CHARACTER*1  STAT,CHAN
      CHARACTER*80 TEXT,COMLIN
      INTEGER KR(18),NR(50)
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
      DATA KX /0,1,1,1,19/
      DATA KXTEXT/' ','X1','X1','X3','HYPOCENTRAL DISTANCE'/
C
C.......................................................................
C
C     Opening data files:
      FSPDAT='sp.dat'
      FILERF='rf.out'
      KOMP=3
      FILSRC=' '
      WRITE(*,'(A)') ' Enter ''sp.dat'', ''rf.out'', KOMP, ''SRC'': '
      READ(*,*) FSPDAT,FILERF,KOMP,FILSRC
      VCOM(1)=0.
      VCOM(2)=0.
      VCOM(3)=0.
      IF(FILSRC.NE.' ') THEN
        OPEN(5,FILE=FILSRC,STATUS='OLD')
        READ(5,*) (TEXT,I=1,20)
        READ(5,*) STAT,VCOM(1),VCOM(2),VCOM(3)
        CLOSE(5)
      END IF
      OPEN(5,FILE=FSPDAT,STATUS='OLD')
      OPEN(4,FILE=FILERF,STATUS='OLD')
C
C     Selection of receivers:
      DO 1 I=1,18
        KR(I)=0
    1 CONTINUE
      READ (5,*) (KR(I),I=1,18)
      MR= 0
      DO 3 K=1,16,3
        L = KR(K)
        M = KR(K+1)
        N = KR(K+2)
        IF(N.EQ.0) THEN
          GO TO 10
        END IF
        DO 2 J=L,M,N
          MR= MR+1
          NR(MR)= J
    2   CONTINUE
    3 CONTINUE
   10 CONTINUE
C
C
C     Initial values for plotting frame:
      TMIN=0.
      TMAX=1.
      TAXIS=10.
      TM=1.
      TN=1.
      CODE=0.
      VRED=0.
      READ (5,*) TMIN,TMAX,TAXIS,TM,TN,CODE,VRED
      XMIN=0.
      XMAX=1.
      XAXIS=FLOAT(MR+1)
      XM=1.
      XN=1.
      YMIN=0.
      YMAX=0.
      HEIGHT=0.4
      READ (5,*) XMIN,XMAX,XAXIS,XM,XN,YMIN,YMAX,HEIGHT
      KOD= INT(CODE+0.5)
      SCY= TAXIS/(TMAX-TMIN)
      XD = XAXIS/FLOAT(MR+1)
      XA = XMAX-XMIN
      YA = YMAX-YMIN
C
C     Plotting frame:
      CALL RGSE1(4,TEXT)
      IX = KOD
      IT = -1
      IF(KOD.GE.3) IX=1
      IF(XM.LT.0.) IX=0
      IF(TM.LT.0.) IT=0
      KX1= KX(KOD+1)+1
      CALL FRAME(XAXIS,TAXIS,ABS(XM),XN,ABS(TM),TN,0,IX,IT,HEIGHT,
     *           XMIN,XMAX,KX1,KXTEXT(KOD+1),YMIN,YMAX,2,'X2',
     *           TMIN,TMAX,4,'TIME',
     *           1,' ',0,0.,0,1,' ',0,0.,0,80,TEXT,1,' ')
C
C     Amplitude scaling:
      AMP=0.
      B1 =1.
      EPIC=1.
      EPS=0.
      READ (5,*) AMP,B1,EPIC,EPS
      IF(B1.EQ.0.)   B1=1.
      IF(EPIC.LE.0.) EPIC=1.
C
C
C     Plotting seismograms:
      JR= 0
      DO 39 IR=1,MR
   31   CONTINUE
          JR= JR+1
   32     CONTINUE
            CALL RGSE2(4,STAT,CHAN,I,XR,YR,ZR,T,TD,NPTS,MPTS,SS)
          IF(I.NE.KOMP) GO TO 32
        IF(JR.LT.NR(IR)) GO TO 31
   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)
C
        RR= SQRT((XR-VCOM(1))**2+(YR-VCOM(2))**2+(ZR-VCOM(3))**2)
        IF(VRED.NE.0.) THEN
          T=T-RR/VRED
        END IF
        IF(AMP.EQ.0.) THEN
          AA=ABS(SS(1))
          DO 35 I=2,NPTS
            AA=AMAX1(ABS(SS(I)),AA)
   35     CONTINUE
          AA= B1*XD/AA
        ELSE
          IF(EPS.EQ.0.) THEN
            AA=B1
          ELSE
            AA=B1*((RR/EPIC)**EPS)
          END IF
        END IF
        X = XD*FLOAT(IR)
        IF(KOD.EQ.1.OR.KOD.EQ.2)
     *  X = XAXIS*((XR-XMIN)*XA+(YR-YMIN)*YA)/(XA*XA+YA*YA)
        IF(KOD.EQ.3) X= XAXIS*(ZR-XMIN)/XA
        IF(KOD.EQ.4) X= XAXIS*(RR-XMIN)/XA
        Y = 0.
        CALL PLOT(X,Y,3)
C
        DO 36 I=1,NPTS
          IF (TMIN.LE.T.AND.T.LE.TMAX) THEN
            A = X-AA*SS(I)
            Y = SCY*(T-TMIN)
            CALL PLOT(A,Y,2)
          END IF
          T = T+TD
   36   CONTINUE
C
        Y = TAXIS
        CALL PLOT(X,Y,2)
   39 CONTINUE
C
      CALL PLOT(0.,0.,999)
      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     LX1= MOD(MOD(KX1(I0),256)+256,256)+1
C     LX2= MOD(MOD(KX2(I0),256)+256,256)+1
C     LY = MOD(MOD(KY (I0),256)+256,256)+1
C     L1 = MOD(MOD(K1 (I0),256)+256,256)+1
C     L2 = MOD(MOD(K2 (I0),256)+256,256)+1
C     L3 = MOD(MOD(K3 (I0),256)+256,256)+1
C     L4 = MOD(MOD(K4 (I0),256)+256,256)+1
C
C     Initialization of CALCOMP:
      CALL PLOTS(0,0,0)
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
        B = YY-YY/YM/2.-HEIGH2
        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 '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