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