INTEGER bi, bj, i, j, xnode, ynode, xq, yq, m, n, p, kx, ky
DO xq = 1,2
gradx(xq) = Xquad(3-xq) * recip_dxG (i,j,bi,bj) +
& Xquad(xq) * recip_dxG (i+1,j,bi,bj)
grady(xq) = Xquad(3-xq) * recip_dyG (i,j,bi,bj) +
& Xquad(xq) * recip_dyG (i,j+1,bi,bj)
xq = 2 - mod(n,2)
if (xq.eq.xnode) kx = 2
& (2*ynode-3) * Xquad(kx) * grady(xq)
& (Xquad(3-xq)*dyG(i,j,bi,bj) + Xquad(xq)*dyG(i+1,j,bi,bj)) *