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