c  gradient.for;no subroutines
c slight modification of  R. Simpson's program
c  calculates the magnitude of the horizontal gradient by fitting a parabola
c  to three adjacent points (which turns out to be the same as taking
c differences between the z and x values of the adjacent points).
c
      character*50 infil,outfil
      character id*56,pgm1*8,pgm2*8
      dimension a(10000),b(10000),c(10000),grad(10000)
      data dcval/1.7e38/,pgm2/'gradient'/
	maxcol=10000
c
1     print*,'*Give input grid:'
      read 801,infil
801   format(a50)
      open(10,file=infil,status='old',form='unformatted',
     &err=5)
      go to 7
5     print*,'Input file not found or wrong type.  Try again'
      go to 1
c
7     if(ncol.gt.maxcol) then
	print 807,maxcol
807	format(' No. cols exceeds max. of ',i5)
	go to 300
      endif
      print*,'*Give output grid:'
      read 801,outfil
      read(10) id,pgm1,ncol,nrow,nz,xo,dx,yo,dy
      print*,id
      print*,pgm1
      print*,'ncol=',ncol,'     nrow=',nrow
      print*,'xo=  ',xo,'     yo=  ',yo
      print*,'dx=  ',dx,'     dy=  ',dy
      irecl=ncol+1
      irecl=max(irecl,30)
	print*
	print*,'Enter output title'
	read(5,802) id
802   format(a56)
      open(11,file=outfil,status='unknown',form='unformatted')
      write(11) id,pgm2,ncol,nrow,nz,xo,dx,yo,dy
c
      read(10) dum,(b(i),i=1,ncol)
      read(10) dum,(c(i),i=1,ncol)
c
      do 10 i=1,ncol
 10   grad(i)=dcval
      write(11) dum,(grad(i),i=1,ncol)
c
      do 200 j=2,nrow-1
      do 120 i=1,ncol
      a(i)=b(i)
 120  b(i)=c(i)
      read(10) dum,(c(i),i=1,ncol)
      do 40 i=1,ncol
 40   grad(i)=dcval
      do 100 i=2,ncol-1
      if(b(i-1).ge.1.0e+38.or.b(i+1).ge.1.0e+38) go to 100
      if(a(i).ge.1.0e+38.or.c(i).ge.1.0e+38) go to 100
      dzdx=(b(i+1)-b(i-1))/(2.0*dx)
      dzdy=(c(i)-a(i))/(2.0*dy)
      grad(i)=sqrt(dzdx*dzdx+dzdy*dzdy)
 100  continue
c
      write(11) dum,(grad(i),i=1,ncol)
 200  continue
c
      do 510 i=1,ncol
 510  grad(i)=dcval
      write(11) dum,(grad(i),i=1,ncol)
c
      close(11)
300   close(10)
      stop
      end
