c sf2prof - grid to saki profile
c
      character ifile*80,unit*1,id*56,pgm*8,ans*1
      common /grid/nc,nr,xo,dx,yo,dy,g(2048)
      data ddval/1.e30/
c
      print'(a\)',' anomaly grid file:'
      read(5,100) ifile
  100 format(a)
      open(10,file=ifile,status='old',form='unformatted')
      read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
      xmax=xo+(nc-1)*dx
      ymax=yo+(nr-1)*dy
      write(6,200) id,pgm,nc,dx,nr,dy,xo,xmax,yo,ymax
  200 format(1h0,a56,1x,a8,/,i5,' cols at ',g12.6,' km intervals',/,
     1i5,' rows at ',g12.6,' km intervals',/,' x-range = ',g12.6,', ',
     2g12.6,/,' y-range = ',g12.6,', ',g12.6)
      if(nc.gt.2048) stop 'too many columns'
      ifile=' '
      print'(a\)',' topography grid file (or return):'
      read(5,100) ifile
      if(ifile(1:1).eq.' ') then
        itopo=0
        print'(a\)',' altitude of data in km (negative upward):'
        read(5,*) alt
        go to 8
      else
        itopo=1
      endif
      open(11,file=ifile,status='old',form='unformatted')
      read(11) id,pgm,nc1,nr1,nz,xo1,dx1,yo1,dy1
      if(nc1.ne.nc.or.nr1.ne.nr.or.xo1.ne.xo.or.yo1.ne.yo.or.dx1.ne.dx
     1   .or.dy1.ne.dy) stop 'grids not registered'
    7 print'(a\)',' units of topography (ft,m,km)?:'
      read(5,100) unit
      fact=0.
      if(unit.eq.'f'.or.unit.eq.'F') fact=-.0003048
      if(unit.eq.'m'.or.unit.eq.'M') fact=-.001
      if(unit.eq.'k'.or.unit.eq.'K') fact=-.1
      if(fact.eq.0.) go to 7
      print'(a\)',' terrain clearance in same units as topography:'
      read(5,*) tc
    8 open(12,status='scratch',form='unformatted',access='direct',
     1     recl=nc*4)
      if(itopo.eq.1) then
        open(13,status='scratch',form='unformatted',access='direct',
     1     recl=nc*4)
      endif
      do 10 j=1,nr
      call rowio(g,yr,nz,nc,10,1)
      call swrda(12,j,g,nc)
      if(itopo.eq.0) go to 10
      call rowio(g,yr,nz,nc,11,1)
      do 9 i=1,nc
      if(g(i).le.ddval) g(i)=(g(i)+tc)*fact
    9 continue
      call swrda(13,j,g,nc)
   10 continue
      close(10)
      if(itopo.eq.1) close(11)
   11 print'(a\)',' output profile file:'
      read(5,100) ifile
      open(14,file=ifile,status='unknown',form='formatted')
      print'(a\)',' enter profile sample interval:'
      read(5,*) del
      tlen=0.
      start=0.
      d0=0.
   12 print'(a\)',' enter first x,y pair:'
      read(5,*) x1,y1
      if(x1.lt.xo.or.x1.gt.xmax.or.y1.lt.yo.or.y1.gt.ymax) then
        print*,'coordinate lies outside the grid'
        go to 12
      endif
   13 print'(a\)',' enter next x,y pair:'
      read(5,*) x2,y2
      if(x2.lt.xo.or.x2.gt.xmax.or.y2.lt.yo.or.y2.gt.ymax) then
        print*,'coordinate lies outside the grid'
        go to 13
      endif
      slen=sqrt((x2-x1)**2+(y2-y1)**2)
      tlen=tlen+slen
      print*,'segment length = ',slen,' total length =',tlen
      npts=(slen-start)/del+1
      c=(x2-x1)/slen
      s=(y2-y1)/slen
      xstart=x1+start*c
      ystart=y1+start*s
      do 20 i=1,npts
      d=d0+(i-1)*del
      x=xstart+(i-1)*c*del
      y=ystart+(i-1)*s*del
      call interp(12,x,y,value)
      if(itopo.eq.0) then
        write(14,300) d,alt,value,x,y
      else
        call interp(13,x,y,t)
        write(14,300) d,t,value,x,y
      endif
   20 continue
  300 format(5g16.6)
      start=slen-d
      d0=d+del
      x1=x2
      y1=y2
      print'(a\)',' more x,y pairs?:'
      read(5,100) ans
      if(ans.eq.'y'.or.ans.eq.'Y') go to 13
      print'(a\)',' additional profiles?:'
      read(5,100) ans
      if(ans.eq.'y'.or.ans.eq.'Y') go to 11
      stop
      end

c***********************************************************************
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
c***********************************************************************
       subroutine srdda(no,ipos,dat[huge],n)
c       subroutine 'srrda' reads keyed sequential files.
      dimension dat(n)
      read(no,rec=ipos)dat
      return
      end
c***********************************************************************
      subroutine swrda(no,ipos,dat[huge],n)
c       subroutine 'swrda' writes keyed sequential files.
      dimension dat(n)
      write(no,rec=ipos)dat
      return
      end
c***********************************************************************
      subroutine interp(iunit,xpt,ypt,value)
      common /grid/nc,nr,xo,dx,yo,dy,g(2048)
      data dval/1.e38/,ddval/1.e30/
c
      x=(xpt-xo)/dx+1.0
      y=(ypt-yo)/dy+1.0
      if(x.lt.1.0.or.x.gt.real(nc)
     1   .or.y.lt.1.0.or.y.gt.real(nr)) then
        value=dval
        return
      endif
      i0=ifix(x)
      j0=ifix(y)
      i1=i0+1
      j1=j0+1
      r=x-real(i0)
      s=y-real(j0)
      call srdda(iunit,j0,g,nc)
      z00=g(i0)
      z10=g(i1)
      call srdda(iunit,j1,g,nc)
      z01=g(i0)
      z11=g(i1)
      if(z00.gt.ddval.or.z10.gt.ddval.or.
     1   z01.gt.ddval.or.z11.gt.ddval) then
         value=dval
      else
         value=(z00*(1-r) + z10*r)*(1-s)
     1        +(z01*(1-r) + z11*r)*s
      endif
      return
      end
