c*****EXTRACTS PROFILE COORDINATES AND VALUES FROM A CONTOUR PLOT
c*****ARRAYS MUST BE DIMENSION 640 IN CALLING PROGRAM
      subroutine profil(igrid,ifile,z,x1,y1,x2,y2,pz,wx,np,ierr)
      integer*4 x1,y1,x2,y2
      integer*2 xx(640),yy(640),np
      character*(*) ifile
      character*64 title
      real*4 z(1),pz(1),wx(1)
      logical opn
      data iprof/7/,dval/1.0e37/
      ierr = 0
      call raster(x1,y1,x2,y2,xx,yy,np)
      read(igrid) title,ncol,nrow,nz,xo,dx,yo,dy
      xo2 = xo + (x1-1) * dx
      xo3 = xo + (x2-1) * dx
      dx2 = (xo3 - xo2) / np
      yo2 = yo + (y1-1) * dy
      yo3 = yo + (y2-1) * dy
      dy2 = (yo3 - yo2) / np
      diffx = x2 - x1
      diffy = y2 - y1
      if(diffy.ge.0..and.diffx.ge.0.) then
        is = 1
        ie = np
        incr = 1
        k = 0
        kinc = 1
      else if(diffy.ge.0.and.diffx.lt.0.) then
        is = 1
        ie = np
        incr = 1
        kinc = - 1
        k = np - kinc
        dx2 = -dx2
        xo2 = xo3
        dy2 = -dy2
        yo2 = yo3
      else if(diffy.lt.0.and.diffx.ge.0.) then
        is = np
        ie = 1
        incr = -1
        kinc = -1
        k = np - kinc
      else if(diffy.lt.0.and.diffx.lt.0.) then
        is = np
        ie = 1
        incr = -1
        k = 0
        kinc = 1
        dx2 = -dx2
        xo2 = xo3
        dy2 = -dy2
        yo2 = yo3
      endif
      do 50 j = 1,nrow
      call rread(igrid,z,ncol)
        do 30 i = is,ie,incr
          if(yy(i).eq.j) then
            k = k + kinc
            pz(k) = z(xx(i))
            if(pz(k).gt.dval) then
              print *,'profile interrupted because large value found'
              ierr = 1
              return
            endif
            wx(k) = float(k)
          else
            go to 40
          endif
30      continue
40    is = i
50    continue
c       print *,'second dx2,dy2,xo2,yo2',dx2,dy2,xo2,yo2
c       print *,'pz,',(pz(i),i=1,5)
c       print *,'wx,',(wx(i),i=1,5)
c       pause
      rewind igrid

      inquire(file='profile.dat',opened=opn)
      if(.not.opn) then
        open(iprof,file='profile.dat',status='unknown')
      endif
      write(iprof,'(2a)') '    GRID FILE NAME: ',ifile
      write(iprof,'(2a)') '    NO. POINTS        X(1)',
     & '         DELTA X          Y(1)         DELTA Y'
      write(iprof,'(i10,5x,4e15.7)') np,xo2,dx2,yo2,dy2 
      nump = np/5
      nrem = np - nump * 5
      is = 1
      ie = 5
      do 60 j = 1,nump
      write(iprof,'(5e15.7)') (pz(i),i=is,ie)
      is = is + 5
      ie = ie + 5
60    continue
      if(nrem.ne.0) then
        write(iprof,'(5e15.7)') (pz(i),i=is,is+nrem-1)
      endif

      return
      end
c*********************************************************************      
      subroutine raster(x1,y1,x2,y2,xx,yy,np)
c*****BRESENHAM'S LINE ALGORITHM
c*****CODED BY R.GODSON-9/84
      implicit integer*2 (a-z)
      integer*4 x1,y1,x2,y2
      dimension xx(1),yy(1)
      if(x2.gt.x1) then
        xinc=1
        dx=x2-x1
      else
        xinc=-1
        dx=x1-x2
      endif
      if(y2.gt.y1) then
        yinc=1
        dy=y2-y1
      else
        yinc=-1
        dy=y1-y2
      endif
      xx(1)=x1
      yy(1)=y1
      if(dx.gt.dy) then
c*******slope between 0 and 1
        incr1=2*dy
        d=incr1-dx
        incr2=2*(dy-dx)
        np=dx+1
        y=y1
c*******loop to set pixels
        do 10 i=2,np
        xx(i)=xx(i-1)+xinc
        if(d.lt.0) then
          d=d+incr1
        else
          d=d+incr2
          y=y+yinc
        endif
        yy(i)=y
   10   continue
      else
c*******slope greater than 1
        incr1=2*dx
        d=incr1-dy
        incr2=2*(dx-dy)
        np=dy+1
        x=x1
c*******loop to set pixels
        do 20 i=2,np
        yy(i)=yy(i-1)+yinc
        if(d.lt.0) then
          d=d+incr1
        else
          d=d+incr2
          x=x+xinc
        endif
        xx(i)=x
   20   continue
      endif
      return
      end
c*********************************************************************      
      subroutine rread(infile,z,ncol)
      dimension z(ncol)
      read(infile) dum,z
      return
      end
