c    @(#)darajac.f	1.1
      subroutine darajac(psi,q,r,m,n,sf)
c
c --- subroutine to compute the arakawa jacobian j(psi,q)=r
c
****** at input
* psi stream function
* q   2nd field
* m,n field size
*
***** at output
* r   jacobian of psi and q
*     
*****
      dimension psi(m,n),q(m,n),r(m,n)
      sfi=sf/12.
      sfs=sf/6.
      sfc=sf/3.
c
c --- interior
c
      m1=m-1
      n1=n-1
      do 15 i=2,m1
        do 10 j=2,n1
          r(i,j)=sfi*((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))
   10   continue
   15 continue
c
c --- left edge
c
      do 20 i=2,m1
        r(i,1)=sfs*((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))
   20 continue
c
c --- bottom
c
      do 30 i=2,n1
        r(m,i)=sfs*((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))
   30 continue
c
c --- right edge
c
      do 40 i=2,m1
        r(i,n)=sfs*((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))
   40 continue
c
c --- top
c
      do 50 i=2,n1
        r(1,i)=sfs*((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))
   50 continue
c
c --- corners
c
      r(1,1)=sfc*((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))
c
      r(m,1)=sfc*((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))
c
      r(m,n)=sfc*((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))
c
      r(1,n)=sfc*((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))
c

      end
