c DRAPE - computes the surface that lies at least 'clear' vertical units
c         above the terrain and has a maximum slope of 'rate' vertical
c         units per horizontal unit.  The slope is discontinuous over ridge
c         crests.
c
c         written by J. Phillips 10/30/95
c
      character*80 ifile
      parameter (NMAX = 2048)
      character*1 drive
      character*56 id
      character*8 pgm
      character*80 prompt
      real zin(NMAX),zout(NMAX)
c
      data ddval/1.e30/
c
      call pfinit('drape')
      call askin
c      print'(a\)',' input terrain grid:'
c      read(5,100) ifile
c  100 format(a)
c      open(10,file=ifile,status='old',form='unformatted')
c      read(10)id,pgm,nc,nr,nz,xo,dx,yo,dy
      ifile=' '
400   prompt='Input terrain grid'
      call askc(prompt,ifile,ierr)
      if(ierr.eq.-2) stop
      call gopen(10,ifile,'old','read',ierr)
      if(ierr.ne.0) stop 'error opening file'
      call gheader('r',10,id,nc,nr,xo,dx,yo,dy,ierr)
      if(ierr.ne.0) stop 'error reading header'
      if (nc.gt.NMAX) stop 'TOO MANY COLUMNS'
      if(dx.ne.dy) stop 'DX MUST EQUAL DY'
c      print'(a\)',' output surface grid:'
c      read(5,100) ifile
c      open(11,file=ifile,status='unknown',form='unformatted')
c      print'(a\)',' title:'
c      read(5,100) id
c      pgm='drape   '
c      write(11)id,pgm,nc,nr,nz,xo,dx,yo,dy
      ifile=' '
401   prompt='Output surface grid'
      call askc(prompt,ifile,ierr)
      if(ierr.eq.-2) then
        close(10)
        go to 400
      endif
      call gopen(11,ifile,'new','write',ierr)
      if(ierr.ne.0) stop 'error opening file'
402   prompt='Title'
      call askc(prompt,id,ierr)
      if(ierr.eq.-2) then
        close(11)
        go to 401
      endif
      call gheader('w',11,id,nc,nr,xo,dx,yo,dy,ierr)
      if(ierr.ne.0) stop 'error writing header'
c      print*,'enter desired minimum terrain clearance in the'
c      print'(a\)','            same vertical units as the terrain:'
c      read*,clear
c      print*,'enter the maximum rate of climb (in vertical'
c      print'(a\)','            terrain units/horizontal terrain units):'
c      read*,rate
c      print'(a\)',' Enter the (RAM) drive letter for temporary files: '
c      read(5,100) drive
      clear=0.0
403   print*,'Enter desired minimum terrain clearance in the'
      prompt='   same vertical units as the terrain'
      call askf4(prompt,clear,ierr)
      if(ierr.eq.-2) go to 402
      rate=1.
404   print*,'Enter the maximum rate of climb (in vertical'
      prompt='   terrain units/horizontal terrain units)'
      call askf4(prompt,rate,ierr)
      if(ierr.eq.-2) go to 403
      drive='c'
405   prompt='Enter the (RAM) drive letter for temporary files'
      call askc(prompt,drive,ierr)
      if(ierr.eq.-2) go to 404
c
c  create a direct access file on the RAM drive
c
      open(12,access='direct',status='unknown',form='unformatted',
     1 file=drive//':\drapein.tmp',recl=nc*4)
      do 10 j=1,nr
      call srdbin(10,zin,nc)
      call swrda(12,j,zin,nc)
   10 continue
c
c  begin main loops
c
      icount=0
      do 33 j=1,nr
      py=yo+(j-1)*dy
      write(6,101)j,nr
      call srdda(12,j,zin,nc)
      do 32 i=1,nc
      px=xo+(i-1)*dx
      pz=zin(i)
      if(pz.gt.ddval) go to 32
      if(j.eq.1) go to 22
      none=0
c
c start from the current row and go downwards until we find
c a row that needs no modifications
c
      do 21 jj=j-1,1,-1
      qy=yo+(jj-1)*dy
      if(none.eq.1) go to 22
      none=1
      call srdda(12,jj,zout,nc)
      do 20 ii=1,nc
      qz=zout(ii)
      if(qz.ge.pz) go to 20
      if(qz.gt.ddval) go to 20
      qx=xo+(ii-1)*dx
      dist=sqrt((px-qx)**2 + (py-qy)**2)
      if(qz.lt.pz-rate*dist) then
        zout(ii)=pz-rate*dist
        icount=icount+1
        none=0
      endif
   20 continue
      call swrda(12,jj,zout,nc)
   21 continue
   22 continue
      none=0
c
c start from the current row and go upwards until we find
c a row that needs no modifications
c
      do 24 jj=j,nr
      qy=yo+(jj-1)*dy
      if(none.eq.1) go to 32
      none=1
      call srdda(12,jj,zout,nc)
      do 23 ii=1,nc
      qz=zout(ii)
      if(qz.ge.pz) go to 23
      if(qz.gt.ddval) go to 23
      qx=xo+(ii-1)*dx
      dist=sqrt((px-qx)**2 + (py-qy)**2)
      if(qz.lt.pz-rate*dist) then
        zout(ii)=pz-rate*dist
        icount=icount+1
        none=0
      endif
   23 continue
      call swrda(12,jj,zout,nc)
   24 continue
   32 continue
   33 continue
c
c  add the clearance and write it out
c
      do 30 j=1,nr
      call srdda(12,j,zout,nc)
      do 25 i=1,nc
      if(zout(i).gt.ddval) go to 25
      zout(i)=zout(i)+clear
   25 continue
      call swrbin(11,zout,nc)
   30 continue
      print*,icount,' modifications'
  101 format(1h+,'working on row ',i5,' of ',i5)
      return
      end
c***********************************************************************
      subroutine srdbin(no,dat[huge],n)
c       subroutine 'srdbin' reads standard grids (binary).
c***********************************************************************
      dimension dat(n)
      read(no)dum,dat
      return
      end
c***********************************************************************
      subroutine swrbin(no,dat[huge],n)
c       subroutine 'swrbin' writes standard grids (binary).
c***********************************************************************
      dimension dat(n)
      data dum/0.0/
      write(no)dum,dat
      return
      end
c***********************************************************************
       subroutine srdda(no,ipos,dat[huge],n)
c       subroutine 'srrda' reads keyed sequential files.
c***********************************************************************
      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.
c***********************************************************************
      dimension dat(n)
      write(no,rec=ipos)dat
      return
      end

