function yderiv_shift, infield, dx, y1d

; This subroutine uses IDL's DERIV.PRO to compute the 
; the 2D y-derivative of the 2D input array.
; Optionally returns 1D array of y-coordinate.

if not(keyword_set(dx)) then dx = 1.

dims = size(infield, /dim)
nx = dims(0)	; define dimensionality
ny = dims(1)

y1d = findgen(ny)*dx

; define arrays
outfield = fltarr(nx,ny)  	; this will be output array

outfield(*,     0) = (infield(*,     1) - infield(*,     0))/dx
outfield(*,  ny-1) = (infield(*,  ny-1) - infield(*,  ny-2))/dx
outfield(*,1:ny-2) = (infield(*,2:ny-1) - infield(*,0:ny-3))/(2*dx)

return, outfield

end