subroutine hgrad(grid,nx,ny,dx,dy,store)
c
c  Subroutine HGRAD calculates the maximum horizontal 
c  gradient of a two-dimensional function using simple finite 
c  difference relations.  Function is specified on a 
c  rectangular grid with x and y axes directed north and east, 
c  respectively.  North is arbitrary.
c
c  Input parameters:
c    nx - number of elements in the west-to-east direction.  
c    ny - number of elements in south-to-north direction.  
c    grid - a singly dimensioned real array containing the 
c         two-dimensional function.  Elements should 
c         be in order of west to east, then south to north (i.e., 
c         element 1 is southwest corner, element ny is 
c         southeast corner, element (nx-1)*ny+1 is northwest 
c         corner, and element ny*nx is northeast corner. 
c    dx - sample interval in the x direction, units irrelevant.
c    dy - sample interval in the y direction, units irrelevant.
c    store - singly-dimensioned real array used internally.
c         Should be dimensioned at least nx*ny in the calling 
c         program.
c
c  Output parameters:
c    grid - upon output, grid contains the maximum horizontal 
c           gradient with same orientation as above.
c
      dimension grid(nx*ny),store(nx*ny)
      index(i,j,ncol)=(j-1)*ncol+i
      dx2=2.*dx
      dy2=2.*dy
      do 10 j=1,nx
         jm1=j-1
         if(jm1.lt.1)jm1=1
         jp1=j+1
         if(jp1.gt.nx)jp1=nx
         do 10 i=1,ny
            im1=i-1
            if(im1.lt.1)im1=1
            ip1=i+1
            if(ip1.gt.ny)ip1=ny
            dfdx=(grid(index(ip1,j,ny))-grid(index(im1,j,ny)))/dx2
            dfdy=(grid(index(i,jp1,ny))-grid(index(i,jm1,ny)))/dy2
            store(index(i,j,ny))=sqrt(dfdx**2+dfdy**2)
   10       continue
      do 20 i=1,nx*ny
   20 grid(i)=store(i)
      return
      end