c            *********************************************
c            *      Topographic Illumination Program     *
c            *********************************************
c
c     Originally coded for Multics by Linda Crownover,
c     Branch of Petrophysics & Remote Sensing,
c     Denver, CO, 1979.
c     Converted to DEC VAX 11/780 by M.H. Podwysocki,
c     U.S.G.S., BP&RS, Sept., 1981.
c     Converted again by J. Phillips, 11/02/83
c
c      parameter nsize = 3000, dval = 1.e37, ddval = '37777677777'o
      parameter (nsize = 3000, dval = 1.e37, ddval = 1.e+38)
      integer ttyout, ttyin
      character ifile*50, ofile*50
      dimension z(nsize), id(14), pgm(2), a(nsize), b(nsize), c(nsize)
      data ttyin, ttyout /5,6/,istop/0/
      data z/nsize*ddval/
      print *,'ifile:'
      read(5,101) ifile
  101 format(a50)
      open(10,file=ifile,form='unformatted',status='old')
      read(10) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      print *,'ofile:'
      read(5,101) ofile
c
c     Get Parameters
c
   54 write (ttyout,55)
   55 format (/1h, 'Enter Solar Altitude and Sun Azimuth in Degrees:',$)
      read (ttyin,*) sunang, azim
      sunang = sunang * 0.0174533
      azim = azim * 0.0174533
      cx = cos (sunang) * cos (azim) / sin (sunang)
      cy = cos (sunang) * sin (azim) / sin (sunang)
c     
      iy=0
c
c     Check Grid Parameters
c
      if (dx .eq. dy) go to 12
      write (ttyout,10)
   10 format(/1h, 'Grid spacing is not square.')
      go to 999
   12 if (ncol .lt. (nsize -1)) go to 14
      write (ttyout,13)
   13 format (/1h, 'Number of columns in grid exceeds array size.')
      go to 999
   14 grd = dx
c
      open(11,file=ofile,form='unformatted',status='new')
      write(11) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      call rowio(ncol,z,0,10,11,istop)
c
      call rowio(ncol,a,-1,10,10,istop)
      call rowio(ncol,b,-1,10,10,istop)
   20 call rowio(ncol,c,-1,10,10,istop)
      if (istop.eq.1) go to 100
      do 30 i = 2, ncol-1
      if (c(i) .ge. dval .or. a(i) .ge. dval) go to 25
      if (b(i-1) .ge. dval .or. b(i+1) .ge. dval) go to 25
c
      x = -(c(i) - a(i)) / (2 * grd) !origin in lower left corner
c
c     x = -(a(i) - c(i)) / (2 * grd) !origin in upperr left corner
      y = (b(i-1) - b(i+1)) / (2 * grd)
      z(i) = 127 * ((1 + (cx * x) + (cy * y)))
      go to 30
   25 z(i) = ddval
   30 continue
      call rowio(ncol,z,0,10,11,istop)
      do 35 i = 1, ncol + 1
      a(i) = b(i)
      b(i) = c(i)
   35 continue
      go to 20
  100 continue
c
c     Close Output File
c
      call rowio(ncol,z,0,10,11,istop)
      close(10)
      close(11)
  999 stop 'Normal termination'
      end
          subroutine rowio(n,z,iop,idev,jdev,iend)
c  WHERE IOP<0 READ; IOP=0 WRITE; IOP>0 READ&WRITE
          dimension z(n)
          iend=0
          if(iop)1,2,1
1         read(idev,end=10) xo,z
          if(iop)9,9,2
2         write(jdev) xo,z
9         return
10        iend=1
          return
          end
