      subroutine scr_p3(i1,i2,screen,npt,spoint,nser,ymin,ymax,
     +     flag,title)
c...  plotting of 3 component seismogram separated by 3
c...  horizontal lines
      dimension screen(1022,3),ymins(3),ymaxs(3),dummy(1022,3)          
      common/scayle/zmin,zmax
      character*(*) title
      integer*4 spoint      
      diffm=ymax-ymin
      if(diffm.gt.0.) goto 6
      diffm=0.
      do 5 i=1,3
      call scr_lim(screen(1,i),npt,ymins(i),ymaxs(i),flag)
    5 diffm=amax1(diffm,ymaxs(i)-ymins(i))    
    6 diffh=0.5*diffm
      zmin=-1.5*diffm
      zmax=1.5*diffm
      do 15 j=1,npt
      dummy(j,1)=amax1(amin1(screen(j,1)+diffm,zmax),diffh)
      dummy(j,2)=amax1(amin1(screen(j,2),diffh),-diffh)
   15 dummy(j,3)=amax1(amin1(screen(j,3)-diffm,-diffh),zmin)
      call scr_plt(i1,i2,dummy(1,1),npt,spoint,nser,zmin,zmax
     +  ,flag,title)
      call scr_plt(1,i2,dummy(1,2),npt,spoint,nser,zmin,zmax
     +  ,flag,title)
      call scr_plt(1,i2,dummy(1,3),npt,spoint,nser,zmin,zmax
     +  ,flag,title)           
      f1=1
      f2=npt
      call scr_hline(diffm,nser)
      return
      end
c
      subroutine scr_hline(diffm,nser)
      x1=1.
      x2=nser
      y1=-0.5*diffm
      y2=y1
      call scr_fastline(x1,y1,x2,y2,0)
      y1=0.5*diffm
      y2=y1
      call scr_fastline(x1,y1,x2,y2,0)
      call ldoit
      return
      end
c
c
      subroutine getran(ilo,ihi,edc,iopt)
      character*(*) edc 
   50 call scr_xhair(ilo,xval,yval,edc,iopt)
      if(edc.eq.'C') return
      if(iopt.ne.0) return
      call scr_xhair(ihi,xval,yval,edc,iopt)
      if(edc.eq.'C') return
      if(iopt.ne.0) return
      if(ihi-ilo.ge.0) return
      call scr_pline('invalid range - try again')
      goto 50                                      
      end
c
c
      integer*4 function indata(io,ifl,spt,blk)
c  moves a block of data into core.this routine tries to read blk
c  points starting at spt but does so in such a way as to allow
c  for decimation.if blk.le.maxtek (1022) then no decimation is
c  necessary and npt (in env$$$) is equal to nptsrd. if blk.gt.
c  maxtek then a decimation factor is calculated by
c        idfact=nptsrd/maxtek+1
c  the program loops through reading idfact*2 point at a time,
c  computing the max and min of these points and storing these
c  in the appropriate slots in the array screen.if indata reaches
c  the end of a file without reading any data it returns zero.
      parameter(nndim=5000)
      integer*4 spoint,epoint,spt,blk,pread
      common/head$/nsta,nchn,iy,id,ih,im,ss,dt,nscan,title
      common/env$$$/spoint,epoint,nser,npt,ymax,ymin,screen(1022,3)
     +   ,il,ix
      common/scrtch/scrbuf(nndim)
      character*80 title
      if(blk.gt.1022) goto 5
      npt=pread(io,ifl,screen(1,1),spt,spt+blk-1)
      npt1=pread(io,ifl+1,screen(1,2),spt,spt+blk-1)
      npt2=pread(io,ifl+2,screen(1,3),spt,spt+blk-1)
      npt=min0(npt,npt1,npt2)
      indata=npt             
      if(indata.le.1) return
      goto 60
c*** too many points - decimate ***
    5 nptsrd=min0(nscan,spt+blk-1)-spt+1
      idfact=nptsrd/1022+1
      idec=2*idfact
      mfac=5000/idec
      mfac=mfac*idec   
      nsave=nptsrd
      do 50 jser=1,3    
      mfl=ifl-1+jser
      nptsrd=nsave
      npt=0
      indata=0
      i1=spt
   10 i2=min0(mfac,nptsrd)+i1-1
      nread=pread(io,mfl,scrbuf,i1,i2)
      indata=indata+nread
      if(nread.le.0) goto 20
      do 15 n=1,nread,idec
      ldec=min0(idec,nread-n+1)
      call decim(scrbuf(n),ldec,screen(npt+1,jser),screen(npt+2,jser))
   15 npt=npt+2
      nptsrd=nptsrd-nread
      if(nptsrd.le.0) goto 20
      i1=i1+nread
      goto 10
   20 if(indata.le.1) return
   50 continue
   60 spoint=spt
      nser=indata
      epoint=spt+nser-1
      return
      end
c
c
      subroutine decim(buf,n,rfirst,rlast)
c  find max and min of a buffer and put in order of occurence
      dimension buf(1)
      imin=1
      imax=1
      do 5 i=1,n
      if(buf(i).gt.buf(imax)) imax=i
    5 if(buf(i).lt.buf(imin)) imin=i
      if(imin.gt.imax) goto 10
      rfirst=buf(imin)
      rlast=buf(imax)
      return
   10 rfirst=buf(imax)
      rlast=buf(imin)
      return
      end
c
c
      subroutine pshenv(lenv)
c  save environment
      common/env$$$/buf(3074)
      common/envsv$/envsav(3074,2)
      lenv=min0(2,lenv+1)
      do 5 i=1,3074
    5 envsav(i,lenv)=buf(i)
      return
      end
c
c
      subroutine popenv(lenv)
c  restore environment
      common/env$$$/buf(3074)
      common/envsv$/envsav(3074,2)
      if(lenv.eq.0) return
      do 5 i=1,3074
    5 buf(i)=envsav(i,lenv)
      lenv=lenv-1
      return
      end
c
c
      logical function nxtf(io,ifl)
c implement the next command. advance to next required file and read header.
      integer*4 spoint,epoint
      common/env$$$/spoint,epoint,nser,npt,ymax,ymin,screen(1022,3)
     +  ,il,ix
      common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,
     + spare(20)
      common/head$/jsta,jchn,jy,jd,jh,jm,ssj,dtj,jscan,title          
      character*80 title,nchn*4
      dimension hbuf(30)
      equivalence (nscan,hbuf)                             
      nxtf=.false.
    1 call gfs_search(io,hbuf,ifl,0)
      if(ifl.lt.0) return
      if(nchn(3:3).eq.'N'.or.nchn(3:3).eq.'R') goto 1
      if(nchn(3:3).eq.'E'.or.nchn(3:3).eq.'T') goto 1
c*** set up internal header 
    5 jsta=nsta
      jchn=nchn
      jy=iy
      jd=id
      jm=im
      jh=ih
      ssj=ss
      dtj=dt                                
      jscan=nscan
      call gfs_rwdir(io,hbuf,ifl+1,'r',ierr)
      jscan=min0(nscan,jscan) 
      if(nchn(3:3).eq.'E'.or.nchn(3:3).eq.'T') goto 1
      if(nchn(3:3).eq.'N'.or.nchn(3:3).eq.'R') goto 6
      ifl=ifl+1
      goto 5
    6 call gfs_rwdir(io,hbuf,ifl+2,'r',ierr)
      jscan=min0(nscan,jscan)               
      if(nchn(3:3).eq.'N'.or.nchn(3:3).eq.'R') goto 1
      if(nchn(3:3).eq.'E'.or.nchn(3:3).eq.'T') goto 7
      ifl=ifl+2
      goto 5
    7 nscan=jscan
      write(title,900) ifl,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,nscan
  900   format(i4,') ',a4,2x,a4,2x,a4,i6,i4,':',i2,':',i2,':',
     &       f6.3,f8.3,i8,'$')
      spoint=1
      epoint=nscan
      nser=nscan
      npt=0
      nxtf=.true.
      return
      end

      logical function nxtbf(io,ifl)
c implement the next command. advance to next required file and read header.
      integer*4 spoint,epoint,qyr,qjday,qhr,qmn
      common/env$$$/spoint,epoint,nser,npt,ymax,ymin,screen(1022,3)
     +  ,il,ix
      common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,
     + qcolat,qcolon,qdep,qyr,qjday,qhr,qmn,qsc,spare(12)
      common/head$/jsta,jchn,jy,jd,jh,jm,ssj,dtj,jscan,title          
      common/dhead$/d0,th,ph,ky,kd,kh,km,ssk,dstuns
      character*80 title,nchn*4
      dimension hbuf(30)
      equivalence (nscan,hbuf)                             
      nxtbf=.false.
    1 call gfs_search(io,hbuf,ifl,0)
      if(ifl.lt.0) return
      if(nchn(3:3).eq.'N'.or.nchn(3:3).eq.'R') goto 1
      if(nchn(3:3).eq.'E'.or.nchn(3:3).eq.'T') goto 1
c*** set up internal header 
    5 jsta=nsta
      jchn=nchn
      jy=iy
      jd=id
      jm=im
      jh=ih
      ssj=ss
      dtj=dt                                
      jscan=nscan
c... the following bit is only for 
c this piece of code assumes source latitudes to be in the
      call gfs_rwdir(io,hbuf,ifl+1,'r',ierr)
      jscan=min0(nscan,jscan) 
      if(nchn(3:3).eq.'E'.or.nchn(3:3).eq.'T') goto 1
      if(nchn(3:3).eq.'N'.or.nchn(3:3).eq.'R') goto 6
      ifl=ifl+1
      goto 5
    6 call gfs_rwdir(io,hbuf,ifl+2,'r',ierr)
      jscan=min0(nscan,jscan)               
      if(nchn(3:3).eq.'N'.or.nchn(3:3).eq.'R') goto 1
      if(nchn(3:3).eq.'E'.or.nchn(3:3).eq.'T') goto 7
      ifl=ifl+2
      goto 5
    7 nscan=jscan
      write(title,900) ifl,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,nscan
  900   format(i4,') ',a4,2x,a4,2x,a4,i6,i4,':',i2,':',i2,':',
     &       f6.3,f8.3,i8,'$')
      spoint=1
      epoint=nscan
      nser=nscan
      npt=0
c data gfs header!!
      d0=qdep
      th=qcolat
      qcolat=90.-qcolat
      ph=qcolon
      ky=qyr
      kd=qjday
      kh=qhr
      km=qmn
      ssk=qsc
      nxtbf=.true.
      return
      end
c
c
      integer*4 function pread(io,ifl,buf,lo,hi)
c  primitive read from input file,this version for gfs
      integer*4 lo,hi
      dimension buf(1)             
      nwd=hi-lo+1
      call gfs_rwdata(io,buf,ifl,lo,nwd,'r',ierr)
      if(ierr.ne.0) nwd=0
      pread=nwd
      return
      end
c
c
      subroutine updat(n,jy,jd,jh,jm,ssj)
      real*8 t,sec,c1,c2,c3,c4,dt
      common/head$/nsta,nchn,iy,id,ih,im,ss,si,nscan,title
      character*80 title
      data c1,c2,c3,c4/60.d0,24.d0,3600.d0,86400.d0/
      sec=ss
      dt=si                                  
      t=((id*c2+ih)*c1+im)*c1+sec+n*dt-dt
      jd=t/c4
      t=t-jd*c4
      jh=t/c3
      t=t-jh*c3
      jm=t/c1
      ssj=t-jm*c1
      jy=iy
    1 idpy=365
      if(jy-4*(jy/4).eq.0) idpy=366
      if(jd.le.idpy) return
      jy=jy+1
      jd=jd-idpy
      goto 1
      end
c
      subroutine axis(xmin,xmax,ymin,ymax,dsam,ityp)
c...  draw tics on x-axis
      rx=alog10(xmax-xmin)
      ix=nint(rx-.5)
      if(rx-ix.gt..5)then
        dx=10**float(ix)
        ibar=2
      else
        dx=10**float(ix)/2.
        ibar=5
      endif
      yrang=ymax-ymin
      yd=ymin+yrang/30.
      yd1=ymin+yrang/60.
      xval0=(int(xmin/dx)+1)*dx
c...  xval1=(int(xmax/dx))*dx
      idx=int(xmax/dx)-int(xmin/dx)
      do 1 i=1,idx
         xval=xval0+(i-1)*dx
         xval1=(xval-ityp*xmin)*dsam
         call scr_fastline(xval1,ymin,xval1,yd,0)
         do 1 j=1,ibar-1
            xvald=xval+j*dx/ibar
            xval1=(xvald-ityp*xmin)*dsam
            call scr_fastline(xval1,ymin,xval1,yd1,0)
 1    continue
      xval0=int(xmin/dx)*dx          
      do 2 j=1,ibar-1
           xvald=xval0+j*dx/ibar
           call scr_fastline(xvald*dsam,ymin,xvald*dsam,yd1,0)
 2    continue
      call ldoit
      return
      end

      subroutine limit(rmin,rmax,a1,a2,a3,npts)
c  find max and min of all data
      dimension a1(1),a2(1),a3(1)
      rmin=1.e+30
      rmax=-1.e+30
      do 5 i=1,npts                          
      rmax=amax1(rmax,a1(i))
      rmin=amin1(rmin,a1(i))
      rmax=amax1(rmax,a2(i))
      rmin=amin1(rmin,a2(i))
      rmax=amax1(rmax,a3(i))
      rmin=amin1(rmin,a3(i))
    5 continue
      return
      end

      subroutine limit1(rmin,rmax,xmin,xmax,x,a1,npts)
c... same as limit, but rmin, rmax are given and only
c... one array is considered
      dimension a1(1),x(1)
      do 5 i=1,npts
        if(x(i).lt.xmin)goto 5
        if(x(i).gt.xmax)goto 6
        rmax=amax1(rmax,a1(i))
        rmin=amin1(rmin,a1(i))
 5    continue
 6    continue
      return
      end

      subroutine limit2(rmin,rmax,xmin,xmax,x,a1,a2,npts)
c... same as limit, but rmin, rmax are given and only
c... one array is considered
      dimension a1(1),a2(1),x(1)
      do 5 i=1,npts
        valh=a1(i)+abs(a2(i))
        vall=a1(i)-abs(a2(i))
        if(x(i).lt.xmin)goto 5
        if(x(i).gt.xmax)goto 6
        rmax=amax1(rmax,valh)
        rmin=amin1(rmin,vall)
 5    continue
 6    continue
      return
      end

      subroutine avbnd(x,y,ery,xmin,xmax,av,er,n)
c...  calculate mean value and error of weighted data
      dimension y(1),x(1),ery(1)             
      SUM1=0.
      SUM2=0.
      do 10 i=1,n
      if(x(i).lt.xmin) goto 10
      if(x(i).gt.xmax) goto 20
      var=ery(i)*ery(i)
      SUM1=SUM1+y(i)/VAR
      SUM2=SUM2+1./VAR
   10 continue
   20 av=0.
      er=0.
      if(sum2.le.0.) return
      AV=SUM1/SUM2   
      ER=1./SQRT(SUM2)
      return
      end

      subroutine avevar(x,y,xmin,xmax,ave,var,n)
c...  calculate mean value and standard deviation
c...  without prior knowledge of number of points
      dimension x(1),y(1)
      ave=0.
      var=0.
      np=0.
      do 1 i=1,n
      if(x(i).lt.xmin) goto 1
      if(x(i).gt.xmax) goto 2
      np=np+1
      ave=ave+y(i)
      var=var+y(i)*y(i)
 1    continue
 2    ave=ave/np
      var=(var-ave*ave*np)/float(np-1)
      if(var.gt.1.e-09)var=sqrt(var)
      return
      end

      subroutine smooth(a1,smth,nout)
      dimension a1(1),smth(1)
      smth(1)=a1(1)
      smth(2)=(a1(1)+a1(2)+a1(3))/3.
      do 1 i=3,nout-2
 1    smth(i)=(a1(i-2)+a1(i-1)+a1(i)+a1(i+1)
     &        +a1(i+2))/5.
      smth(nout-1)=(a1(nout)+a1(nout-1)+a1(nout-2))/3.
      smth(nout)=a1(nout)
      return
      end

      double precision function geocen(arg)
c Returns geocentric colat given geographic colat as arg
c  (n.b. fac=(1-f)**2,   argument in radians.
      implicit real*8(a-h,o-z)
      data pi2,fac/1.570796326794895d0,0.993305621334896d0/
      geocen=pi2-datan(fac*dcos(arg)/dmax1(1.d-200,dsin(arg)))
      return
      end

      subroutine dist(adis,azi,tstart,fmin0,fmax,ffiler,ffilet)
c...  evaluation of epicentral distance and beginning
c...  of record rel. to origin time
      parameter(no=7)
      real*8 th0,geocen
      integer*4 nper
      character*80 ffiler,ffilet
      common/trav/trminr(no),trmaxr(no),trminl(no),trmaxl(no)
      common/srce/jy,jd,jh,jm,sec,th,ph,d0
      common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,spare(20)
      data pi/3.14159265358979/
      rad=180./pi
      fmin=fmin0
      if(fmin.lt..0028)fmin=.0028
c     if(ih.eq.23.and.im.gt.30)id=id+1
      call sourcep
c     if(ih.eq.23.and.im.gt.30)id=id-1
      adis=0.
      tstart=0.
      if(iy.lt.0) return
      t0=th/rad
      p0=ph/rad
      th0=t0
      t0=geocen(th0)
      c0=cos(t0)
      s0=sin(t0)                         
      t1=0.
      p1=0.
      call allnsta(nsta,t1,p1)
      if(abs(t1-1.57079).lt.1.e-04.and.abs(p1).lt.1.e-04)then
        adis=0.
      else
        th0=t1
        t1=geocen(th0)
        c1=cos(t1)
        s1=sin(t1)
        dp=p1-p0
        co=c0*c1+s0*s1*cos(dp)
        si=dsqrt(1.d0-co*co)
        saz=-s0*sin(dp)/si
        caz=(c0-co*c1)/(si*s1)
        azi=atan2(saz,caz)*rad
        adis=atan2(si,co)*6371.
      endif
      tstart=((((iy-jy)*366.+(id-jd))*24.+(ih-jh))*60.+(im-jm))*60.+
     +   ss-sec
      if(ffiler.eq.' ')goto 40
c...  get dispersion model
      open(2,file=ffiler)
c...  calculate travel time range of rayleigh waves
      read(2,1)dummy
      read(2,*)nper
      if(nper.eq.0)nper=100000
      read(2,1)dummy
      read(2,1)dummy
      do 50 i=1,no
      trminr(i)=1.e20
      trmaxr(i)=-1.e20
 50   continue
      do 10 i=1,nper
         read(2,*,end=30)fdum,ug
         if(1./fdum.lt.fmin)goto 10
         if(1./fdum.gt.fmax)goto 30
         trtime=adis/ug-tstart
         if(trtime.lt.trminr(1))then
           trminr(1)=trtime
           trminr(2)=(2.*pi*6371.-adis)/ug-tstart
           trminr(3)=trminr(1)+2.*pi*6371./ug
           trminr(4)=(4.*pi*6371.-adis)/ug-tstart
           trminr(5)=trminr(1)+4.*pi*6371./ug
           trminr(6)=(6.*pi*6371.-adis)/ug-tstart
           trminr(7)=trminr(1)+6.*pi*6371./ug
         endif
         if(trtime.gt.trmaxr(1))then
           trmaxr(1)=trtime
           trmaxr(2)=(2.*pi*6371.-adis)/ug-tstart
           trmaxr(3)=trmaxr(1)+2.*pi*6371./ug
           trmaxr(4)=(4.*pi*6371.-adis)/ug-tstart
           trmaxr(5)=trmaxr(1)+4.*pi*6371./ug
           trmaxr(6)=(6.*pi*6371.-adis)/ug-tstart
           trmaxr(7)=trmaxr(1)+6.*pi*6371./ug
         endif
 10   continue
 30   continue
      close(2)
      open(3,file=ffilet)
c...  calculate travel time range of love waves
      read(3,1)dummy
      read(3,*)nper
      if(nper.eq.0)nper=100000
      read(3,1)dummy
      read(3,1)dummy
      do 60 i=1,no
      trminl(i)=1.e20
      trmaxl(i)=-1.e20
 60   continue
      do 20 i=1,nper
         read(3,*,end=40)fdum,ug
         if(1./fdum.lt.fmin)goto 20
         if(1./fdum.gt.fmax)goto 40
         trtime=adis/ug-tstart
         if(trtime.lt.trminl(1))then
           trminl(1)=trtime
           trminl(2)=(2.*pi*6371.-adis)/ug-tstart
           trminl(3)=trminl(1)+2.*pi*6371./ug
           trminl(4)=(4.*pi*6371.-adis)/ug-tstart
           trminl(5)=trminl(1)+4.*pi*6371./ug
           trminl(6)=(6.*pi*6371.-adis)/ug-tstart
           trminl(7)=trminl(1)+6.*pi*6371./ug
         endif
         if(trtime.gt.trmaxl(1))then
           trmaxl(1)=trtime
           trmaxl(2)=(2.*pi*6371.-adis)/ug-tstart
           trmaxl(3)=trmaxl(1)+2.*pi*6371./ug
           trmaxl(4)=(4.*pi*6371.-adis)/ug-tstart
           trmaxl(5)=trmaxl(1)+4.*pi*6371./ug
           trmaxl(6)=(6.*pi*6371.-adis)/ug-tstart
           trmaxl(7)=trmaxl(1)+6.*pi*6371./ug
         endif
 20   continue
 40   continue
      close(3)
 1    format(a)
      return
      end

      subroutine sourcep
      common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,slat,slon,
     &             sdep,jys,jds,jhs,jms,sss,spare(12)

      common/srce/jy,jd,jh,jm,sec,t0,p0,d0
c...  transfer source information             
      jy=jys
      jd=jds
      jh=jhs
      jm=jms
      sec=sss
      t0=90.-slat
      p0=slon
      d0=sdep
      return
      end

      subroutine plt_trav(spoint)                     
c*** plot lines at arrival times
      parameter(no=7)
      common/scayle/zmin,zmax
      common/trav/trminr(no),trmaxr(no),trminl(no),trmaxl(no)

      integer*4 spoint

      do 10 i=1,no
      tmin=trminr(i)-spoint
      tmax=trmaxr(i)-spoint
      call scr_fastline(tmin,zmin,tmin,zmax,0)
      call scr_fastline(tmax,zmin,tmax,zmax,0)
      tmin=trminl(i)-spoint
      tmax=trmaxl(i)-spoint
      call scr_fastline(tmin,zmin,tmin,zmax,5)
      call scr_fastline(tmax,zmin,tmax,zmax,5)
      do 10 j=1,i
      cmin=zmin+(zmax-zmin)*j/100.
      call scr_fastline(tmin-10.,cmin,tmin+10.,cmin,0)
 10   continue
      call ldoit
      return
      end

      subroutine scr_pl2(icl,b,x,npt,ymn,ymx,yval,xmn,xmx,title,ityp)
c
c  generate a screen plot - this version allows no decimation.
c  Plots b(x) for b(1)...b(npt). The first call for a plot must be made
c  with icl=0, afterwards segments of the initial plot may be replotted
c  or overplotted by calling with icl.ne.0. The segment to be replotted
c  is expected to be in b(1) to b(npt) where npt is now the number of
c  points in the segment
c
c  if icl=0 clear screen before plotting
c  if icl>0 overplot segment
c  if icl<0 replot segment
c
c...npt = input -- number of points
c   x   = array of x-data points
c   b   = array of y-data points
c   xmn = input -- first  (chosen) start value of x-data array
c   xmx = input -- last value...
c   ymn = input -- minimum value of b-data array
c   ymx = input -- maximum value...
c   title = name of the plot
      parameter (iymin=6,iymax=603,ixmin=6,ixmax=1017,ix0=5,iy0=5)
      parameter (nch=7,ncht=13,nchmax=1012/nch)
      character*(nchmax) text
      character ylab*12
      character*(*) title
      real b(*),x(*)
      common/scr1$$/iscrx(2000),iscry(2000)
      integer cent(2)
      common/scr$$$/idraw,ifill,flags,xfac,bxp,yfac,byp,nropt,jbeg
      if(npt.le.1) return
c... set up plot scale factors whenever screen is cleared 
      if(icl.eq.0) then
        ypmin=iymin
        ypmax=iymax
        xpmin=ixmin
        xpmax=ixmax
        flags=1.e+33
        ymaxs=ymx
        ymins=ymn
        xmaxs=xmx
        xmins=xmn
        if(ymaxs.le.ymins) call scr_lim(b,npt,ymins,ymaxs,flags)
        if(xmaxs.le.xmins) call scr_lim(x,npt,xmins,xmaxs,flags)
        if(xmaxs.le.xmins) return
        if(ymaxs.le.ymins) then
               ymaxs=ymaxs+0.5
               ymins=ymins-0.5
        end if
        xfac=(xpmax-xpmin)/(xmaxs-xmins)
        yfac=(ypmax-ypmin)/(ymaxs-ymins)
        jbeg=1
        call ldrawfilledbox(0,0,1024,619)
        call lsetcolor(idraw)
        call ldrawhollowbox(ix0,iy0,ixmax-ixmin+2,iymax-iymin+2)
        call lsetcolor(ifill)
      end if
c... determine x and y pixel range to be plotted
      byp=iymin+0.5+yfac*ymaxs
      bxp=ixmin+0.5-xfac*xmins
      ip1=bxp+xfac*x(1)
      ip2=bxp+xfac*x(npt)
      ip1=max(ixmin,ip1)
      ip2=min(ixmax,ip2)
      if(icl .lt. 0) call scr_clrseg(ip1,ip2)
      if(ifill.eq.0) then
         call lsetcolor(1)
      else
         call lsetcolor(0)
      end if
c*****************************************************************      
c...Plot symbols
c    Isymb will be expressed in modulo 10, and the  numbers
c    less then 10 represent a symbol for every datapoint b(i).
c    For more information see subroutine scr_symb
c       
9030  isymb=-iopt
      mode=0
9031  If (isymb.le.10) goto 9032
      mode=mode+1
      isymb=isymb-10
      goto 9031
9032  if (mode.gt.3) mode=0
      do 9300 i=1,npt
             xp=xfac*x(i)+bxp
             if (xp.lt.xpmin.or.xp.gt.xpmax) goto 9300
             yp=-yfac*yval+byp
             if (yp.lt.ypmin.or.yp.gt.ypmax) goto 9300
             cent(1)=xp
c... term "+ityp" in next line shifts limit down to .1 for filled symbols
c... for plotting components (ityp=0 for sing. values)
             ityp1=0
             if(ityp.gt.0)ityp1=1
             ib=b(i)**(ityp1+1)/0.1+7*ityp1
             ib=ib-2
             if(ib.lt.1)ib=1
             if(ib.gt.8)ib=8
             goto (10,20,30,40,50,60,70,90)ib
c... open circle for d1<0.7
 10          isymb=1+ityp
             rad=5.0
             goto 80
 20          isymb=1+ityp
             rad=4.0
             goto 80
 30          isymb=1+ityp
             rad=3.0
             goto 80
 40          isymb=1+ityp
             rad=2.0
             goto 80
c... filled circle for d1>0.7
 50          isymb=2+ityp
             rad=2.0
             goto 80
 60          isymb=2+ityp
             rad=3.0
             goto 80
 70          isymb=2+ityp
             rad=4.0
             goto 80
 90          isymb=2+ityp
             rad=5.0
 80          continue
             cent(2)=yp
             call scr_sym2(cent,isymb,rad)
9300  continue
c*****************************************************************
9050  call lsetcolor(ifill)
      if(icl .ne. 0) goto 9999
c...Plot labels if required
c set character height 
c write the labels for ymin and ymax
      write (ylab,9900) ymaxs
      call scr_pltchr(ix0+2,ncht+iy0+1,ylab)
      write (ylab,9900) ymins  
      call scr_pltchr(ix0+2,iymax-1,ylab)      
9900  format(g12.4)
      nlen=min(len(title),nchmax-24)
      text(1:nchmax)=' '
      nstr=0   
      do 51 i=1,nlen
      if(title(i:i).eq.'$') goto 55
   51 nstr=nstr+1
   55 write(text(1:12),9900) xmins
      write(text(nchmax-11:nchmax),9900) xmaxs
      kc=0
      do 550 i=1,12
      if(text(i:i).ne.' ') goto 551
  550 kc=kc+1  
  551 do 552 i=1,12
      ii=i+kc  
  552 text(i:i)=text(ii:ii)
      kc=0
      do 553 i=nchmax,nchmax-11,-1
      if(text(i:i).ne.' ') goto 554
  553 kc=kc+1
  554 do 555 i=nchmax,nchmax-11,-1
      ii=i-kc
  555 text(i:i)=text(ii:ii)
      mbeg=nchmax/2-nstr/2
      mend=mbeg+nstr-1
      text(mbeg:mend)=title(1:nstr)
      call scr_pltchr(ix0+2,iymax+ncht+3,text)
9999  call ldoit
      return
      end

      subroutine scr_sym2(center,isymb,crad)
c*** internal routine for plotting symbols, the symbol
c is plotted centered at 'center' (a two element integer array)
c  and isymb determines the type of symbol (see below)
c
c    uses general drawing primitive
c    primitive defaults values and # of points are
c    prim. = 0 -- circle              |    npt = 2
c          = 1 -- circular arc        |        = 3
c          = 2 -- rectangle           |        = 2
c          = 3 -- circular disk       |        = 2
c          = 4 -- closed circular arc |        = 3
c          = 5 -- solid  rectangle    |        = 2

      integer center(2),rad1,rad2
      data rad1,rad2/1,2/
      irad=crad
      goto (10,60,50,30,40,20,70),isymb
c*** open circle
   10 call ldrawhollowcircle(center(1),center(2),irad)
      return    
c*** x   
   20 ix=center(1)-irad
      iy=center(2)+irad
      ix1=ix-2*irad
      iy1=iy-2*irad
      call lmoveto(ix,iy)
      call llineto(ix1,iy1)
      call lmoveto(ix1,iy)
      call llineto(ix,iy1)
      call lmoveto(center(1),center(2))
      return
c*** +
   30 iy=center(2)+irad
      ix=center(1)-irad
      call lmoveto(center(1),iy)
      call llineto(center(1),iy-2*irad)
      call lmoveto(ix,center(2))
      call llineto(ix+2*irad,center(2))
      call lmoveto(center(1),center(2))
      return
c*** *
   40 ix=center(1)+irad
      iy=center(2)+irad
      ix1=ix-2*irad
      iy1=iy-2*irad
      call lmoveto(ix,iy)
      call llineto(ix1,iy1)
      call lmoveto(ix1,iy)
      call llineto(ix,iy1)
      call lmoveto(center(1),iy)
      call llineto(center(1),iy1)
      call lmoveto(ix1,center(2))
      call llineto(ix,center(2))
      call lmoveto(center(1),center(2))
      return             
c*** open square
   50 call ldrawhollowbox(center(1)-irad,center(2)-irad,2*irad,2*irad)
      return
c*** filled circle
   60 call ldrawfilledcircle(center(1),center(2),irad)
      return
c*** filled square 
   70 call ldrawfilledbox(center(1)-irad,center(2)-irad,2*irad+1,
     +    2*irad+1)
      return
      end

      subroutine scr_fastline(x1,y1,x2,y2,idash)
c...Plot straight line with color icol and dash length idash pixels
c*** coordinates are user coordinates (not pixels)
c*** because of raster ops -- surprising colors will result
      parameter (iymin=6,iymax=603,ixmin=6,ixmax=1017)
      common/scr$$$/idraw,ifill,flags,xfac,bxp,yfac,byp,nropt,jbeg
      iy1=-y1*yfac+byp
      iy1=max(min(iy1,iymax),iymin)
      iy2=-y2*yfac+byp
      iy2=max(min(iy2,iymax),iymin)
      ix1=x1*xfac+bxp
      ix1=max(min(ix1,ixmax),ixmin)
      ix2=x2*xfac+bxp
      ix2=max(min(ix2,ixmax),ixmin)
      call lsetcolor(10)
      call llinestyle(0,idash,0,0)
      call lsetrasterops(6) 
      call lmoveto(ix1,iy1) 
      call llineto(ix2,iy2) 
      call lsetrasterops(3) 
      call llinestyle(0,0,0,0)
      call lsetcolor(ifill)
      return
      end
