1      REAL*8 AA(513,513), DD(513,513)
2      REAL*8 X (513,513), Y (513,513)
3      REAL*8 RX(513,513), RY(513,513)
4      DO J = 2,N-1
5         DO I = 2,N-1
6            XX = X(I+1,J)-X(I-1,J)
7            YX = Y(I+1,J)-Y(I-1,J)
8            XY = X(I,J+1)-X(I,J-1)
9            YY = Y(I,J+1)-Y(I,J-1)
10            A = 0.25 * (XY*XY+YY*YY)
11            B = 0.25* (XX*XX+YX*YX)
12            C = 0.125 * (XX*XY+YX*YY)
13            AA(I,J) = -B
14            DD(I,J) = B+B+A*REL
15            PXX = X(I+1,J)-2.*X(I,J)+X(I-1,J)
16            QXX = Y(I+1,J)-2.*Y(I,J)+Y(I-1,J)
17            PYY = X(I,J+1)-2.*X(I,J)+X(I,J-1)
18            QYY = Y(I,J+1)-2.*Y(I,J)+Y(I,J-1)
19            PXY = X(I+1,J+1)-X(I+1,J-1)-X(I-1,J+1)+X(I-1,J-1)
20            QXY = Y(I+1,J+1)-Y(I+1,J-1)-Y(I-1,J+1)+Y(I-1,J-1)
21            RX(I,J) = A*PXX+B*PYY-C*PXY
22            RY(I,J) = A*QXX+B*QYY-C*QXY
23         END DO
24      END DO
25