c  normal.for
c
c  To calculate and output in three grids the three elements of
c  the unit normal vector for a grid, to be used for shaded relief.
c  Reference is to the x=east, y=north, z=up axes.  In general, n1
c  and n2 in [-1,1], n3 in [0,1).  Often n3<<1, and needs to be
c  augmented to give a pleasing shadow effect.
c
c  Lin Cordell 22 Aug 90.
      common x0,dx,y0,dy
      character*50 ifile,ofile
      dimension f(2000),fx(2000),w1(2000),w2(2000),w3(2000)
      dimension id(14),pgm(2)
      write(6,1)
1     format(1x,'Enter input grid name:'/' *'$)
      read(5,100)ifile
100   format(a)
      open(10,file=ifile,form='unformatted',status='old')
      read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      if(ncol.gt.2000.or.nz.ne.1.or.dx.eq.0.e0)then
        print*,'Cant handle it.'
        stop
      endif
      print*,'Enter n1 file name:'
      print 101
101   format(' *'$)
      read(5,100)ofile
      open(12,file=ofile,form='unformatted',status='unknown')
      write(12)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      print*,'Enter n2 file name:'
      print 101
      read(5,100)ofile
      open(13,file=ofile,form='unformatted',status='unknown')
      write(13)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      print*,'Enter n3 file name:'
      print 101
      read(5,100)ofile
      open(14,file=ofile,form='unformatted',status='unknown')
      write(14)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      wmin=1.e38
      do 10 j=1,nrow
      read(10) y1,(w1(i),i=1,ncol)
      do 10 i=1,ncol
   10 if(w1(i).lt.wmin) wmin=w1(i)
      rewind(10)
      read(10) id
      c=1.0
      call work(0,wmin,c,avslp,f,fx,w1,w2,w3,ncol,nrow)
      print*,'Average slope = ',avslp
      print*,'Enter desired average slope:'
      read(5,*) c
      c=c/avslp
      rewind(10)
      read(10) id
      call work(1,wmin,c,avslp,f,fx,w1,w2,w3,ncol,nrow)
      print*,'Average slope = ',avslp
      stop
      end
      subroutine work(iflag,wmin,c,avslp,gy,gx,w1,w2,w3,ncol,nrow)
      dimension gy(ncol),gx(ncol),w1(ncol),w2(ncol),w3(ncol)
      common x0,dx,y0,dy
      ddval=1.e37
      dval=1.e38
      bx=2.*dx/c
      by=2.*dy/c
      avslp=0.
      numslp=0
c First two rows.
      j=2
      read(10)y1,w1
      read(10)y2,w2
      read(10)y3,w3
c      do 2 i=1,ncol
c      if(w1(i).ge.ddval) w1(i)=wmin
c      if(w2(i).ge.ddval) w2(i)=wmin
c    2 if(w3(i).ge.ddval) w3(i)=wmin
      do 3 i=2,ncol-1
      if(w3(i).ge.ddval.or.w1(i).ge.ddval) go to 2
      if(w2(i+1).ge.ddval.or.w2(i-1).ge.ddval) go to 2
      gy(i)=(w3(i)-w1(i))/by
      gx(i)=(w2(i+1)-w2(i-1))/bx
      denom=gx(i)*gx(i)+gy(i)*gy(i)+1.
      avslp=avslp+sqrt(denom-1.0)
      numslp=numslp+1
      if(iflag.eq.1) then
        w1(i)=1./sqrt(denom)
        gx(i)=-gx(i)*w1(i)
        gy(i)=-gy(i)*w1(i)
      endif
      go to 3
2     w1(i)=dval
      gx(i)=dval
      gy(i)=dval
3     continue
      if(iflag.eq.1) then
        gx(1)=gx(2)
        gx(ncol)=gx(ncol-1)
        gy(1)=gy(2)
        gy(ncol)=gy(ncol-1)
        w1(1)=w1(2)
        w1(ncol)=w1(ncol-1)
        write(12)y1,gx
        write(13)y1,gy
        write(14)y1,w1
        write(12)y2,gx
        write(13)y2,gy
        write(14)y2,w1
      endif
c Rest of the rows.
4     j=j+1
      do 7 i=1,ncol
      w1(i)=w2(i)
7     w2(i)=w3(i)
      read(10)y3,w3
c      do 6 i=1,ncol
c    6 if(w3(i).ge.ddval) w3(i)=wmin
      do 5 i=2,ncol-1
      if(w3(i).ge.ddval.or.w1(i).ge.ddval) go to 6
      if(w2(i+1).ge.ddval.or.w2(i-1).ge.ddval) go to 6
      gy(i)=(w3(i)-w1(i))/by
      gx(i)=(w2(i+1)-w2(i-1))/bx
      denom=gx(i)*gx(i)+gy(i)*gy(i)+1.
      avslp=avslp+sqrt(denom-1.0)
      numslp=numslp+1
      if(iflag.eq.1) then
        w1(i)=1./sqrt(denom)
        gx(i)=-gx(i)*w1(i)
        gy(i)=-gy(i)*w1(i)
      endif
      go to 5
6     w1(i)=dval
      gx(i)=dval
      gy(i)=dval
5     continue
      if(iflag.eq.1) then
        gx(1)=gx(2)
        gx(ncol)=gx(ncol-1)
        gy(1)=gy(2)
        gy(ncol)=gy(ncol-1)
        w1(1)=w1(2)
        w1(ncol)=w1(ncol-1)
        write(12)y2,gx
        write(13)y2,gy
        write(14)y2,w1
      endif
      if(j.lt.(nrow-1))go to 4
c  Last row.
      if(iflag.eq.1) then
        write(12)y3,gx
        write(13)y3,gy
        write(14)y3,w1
      endif
      avslp=avslp/float(numslp)
      return
      end
