c    @(#)arajac.f	1.1
c
c
c
      subroutine arajac(psi,q,r,m,n)
      dimension psi(m,n),q(m,n),r(m,n)
c
c --- subroutine to compute the arakawa jacobian j(psi,q)=r
c
c --- interior
c
      m1=m-1
      n1=n-1
      do 15 i=2,m1
        do 10 j=2,n1
          r(i,j)=((psi(i-1,j)-psi(i,j-1))*q(i-1,j-1)+
     /      (psi(i-1,j-1)+psi(i-1,j)-psi(i+1,j-1)-psi(i+1,j))*q(i,j-1)+
     /      (psi(i,j-1)-psi(i+1,j))*q(i+1,j-1)+
     /      (psi(i-1,j+1)+psi(i,j+1)-psi(i-1,j-1)-psi(i,j-1))*q(i-1,j)+
     /      (psi(i,j-1)+psi(i+1,j-1)-psi(i,j+1)-psi(i+1,j+1))*q(i+1,j)+
     /      (psi(i,j+1)-psi(i-1,j))*q(i-1,j+1)+
     /      (psi(i+1,j)+psi(i+1,j+1)-psi(i-1,j)-psi(i-1,j+1))*q(i,j+1)+
     /      (psi(i+1,j)-psi(i,j+1))*q(i+1,j+1))/12.
   10   continue
   15 continue
c
c --- left edge
c
      do 20 i=2,m1
        r(i,1)=((psi(i-1,2)+psi(i,2)-2.*psi(i,1))*q(i-1,1)+
     /     (2.*(psi(i-1,1)-psi(i+1,1)))*q(i,1)+
     /     (2.*psi(i,1)-psi(i,2)-psi(i+1,2))*q(i+1,1)+
     /     (psi(i,2)-psi(i-1,1))*q(i-1,2)+
     /     (psi(i+1,1)+psi(i+1,2)-psi(i-1,1)-psi(i-1,2))*q(i,2)+
     /     (psi(i+1,1)-psi(i,2))*q(i+1,2))/12.
   20 continue
c
c --- bottom
c
      do 30 i=2,n1
        r(m,i)=((psi(m-1,i-1)+psi(m-1,i)-2.*psi(m,i))*q(m,i-1)+
     /     (2.*(psi(m,i-1)-psi(m,i+1)))*q(m,i)+
     /     (2.*psi(m,i)-psi(m-1,i)-psi(m-1,i+1))*q(m,i+1)+
     /     (psi(m-1,i)-psi(m,i-1))*q(m-1,i-1)+
     /     (psi(m,i+1)+psi(m-1,i+1)-psi(m,i-1)-psi(m-1,i-1))*q(m-1,i)+
     /     (psi(m,i+1)-psi(m-1,i))*q(m-1,i+1))/12.
   30 continue
c
c --- right edge
c
      do 40 i=2,m1
        r(i,n)=((psi(i+1,n1)+psi(i,n1)-2.*psi(i,n))*q(i+1,n)+
     /     (2.*(psi(i+1,n)-psi(i-1,n)))*q(i,n)+
     /     (2.*psi(i,n)-psi(i,n1)-psi(i-1,n1))*q(i-1,n)+
     /     (psi(i,n1)-psi(i+1,n))*q(i+1,n1)+
     /     (psi(i-1,n)+psi(i-1,n1)-psi(i+1,n)-psi(i+1,n1))*q(i,n1)+
     /     (psi(i-1,n)-psi(i,n1))*q(i-1,n1))/12.
   40 continue
c
c --- top
c
      do 50 i=2,n1
        r(1,i)=((psi(2,i+1)+psi(2,i)-2.*psi(1,i))*q(1,i+1)+
     /     (2.*(psi(1,i+1)-psi(1,i-1)))*q(1,i)+
     /     (2.*psi(1,i)-psi(2,i)-psi(2,i-1))*q(1,i-1)+
     /     (psi(2,i)-psi(1,i+1))*q(2,i+1)+
     /     (psi(1,i-1)+psi(2,i-1)-psi(1,i+1)-psi(2,i+1))*q(2,i)+
     /     (psi(1,i-1)-psi(2,i))*q(2,i-1))/12.
   50 continue
c
c --- corners
c
      r(1,1)=((2.*(psi(1,2)-psi(2,1)))*q(1,1)+
     /   (2.*psi(1,1)-psi(2,2)-psi(1,2))*q(2,1)+
     /   (psi(2,1)-psi(1,2))*q(2,2)+
     /   (psi(2,1)+psi(2,2)-2.*psi(1,1))*q(1,2))/12.
c
      r(m,1)=((2.*(psi(m-1,1)-psi(m,2)))*q(m,1)+
     /   (2.*psi(m,1)-psi(m-1,2)-psi(m-1,1))*q(m,2)+
     /   (psi(m,2)-psi(m-1,1))*q(m-1,2)+
     /   (psi(m,2)+psi(m-1,2)-2.*psi(m,1))*q(m-1,1))/12.
c
      r(m,n)=((2.*(psi(m,n1)-psi(m-1,n)))*q(m,n)+
     /   (2.*psi(m,n)-psi(m-1,n1)-psi(m,n1))*q(m-1,n)+
     /   (psi(m-1,n)-psi(m,n1))*q(m-1,n1)+
     /   (psi(m-1,n)+psi(m-1,n1)-2.*psi(m,n))*q(m,n1))/12.
c
      r(1,n)=((2.*(psi(2,n)-psi(1,n1)))*q(1,n)+
     /   (2.*psi(1,n)-psi(2,n1)-psi(2,n))*q(1,n1)+
     /   (psi(1,n1)-psi(2,n))*q(2,n1)+
     /   (psi(1,n1)+psi(2,n1)-2.*psi(1,n))*q(2,n))/12.
c
      return
      end
