c
c @(#)agrad2.f	1.1
c
      subroutine agrad2(psi,dx2psi,dy2psi,m,n,hy)
c
c --- this routine calculates the fourth order laplacian of psi
c --- at the boundaries and in the interior field
c
c
       dimension psi(m*n),dx2psi(m*n),dy2psi(m*n)
c
       msq=m*n
       fac=1./12./hy**2
        ip=0
      do 100 j=1,n
      do 110 i=1,m
      ip=ip+1
      if(i.eq.1.or.i.eq.m) go to 50
      if(i.eq.2.or.i.eq.m-1) go to 60
      dx2psi(ip)=fac*(-psi(ip-2)+16.*psi(ip-1)-30.*psi(ip)
     +                 +16.*psi(ip+1)-psi(ip+2))
      go to 70
c
  50  continue
      nw=1
      if(i.eq.m)nw=-1
      dx2psi(ip)=fac*(45.*psi(ip)-10.*psi(ip+nw*5)+61.*psi(ip+nw*4)
     +                -156.*psi(ip+3*nw)+214.*psi(ip+2*nw)
     +                -154.*psi(ip+nw))
       go to 70
 60   continue
      nw=1
      if(i.eq.m-1)nw=-1
      dx2psi(ip)=fac*(-15.*psi(ip)+psi(ip+4*nw)-6.*psi(ip+nw*3)
     +                 +14.*psi(ip+2*nw)-4.*psi(ip+nw)
     +                 +10.*psi(ip-nw))
c
  70  continue
      if(j.eq.1.or.j.eq.n) go to 80
      if(j.eq.2.or.j.eq.n-1) go to 90
      dy2psi(ip)=fac*(-psi(ip-2*m)+16.*psi(ip-m)-30.*psi(ip)
     +                +16.*psi(ip+m)-psi(ip+2*m))
      go to 110
c
  80  continue
      nw=1
      if(j.eq.n)nw=-1
      dy2psi(ip)=fac*(45*psi(ip)-10.*psi(ip+5*nw*m)+61.*psi(ip+4*nw*m)
     +                -156.*psi(ip+3*nw*m)+214.*psi(ip+2*m*nw)
     +                -154.*psi(ip+nw*m))
      go to 110
 90   continue
      nw=1
      if(j.eq.n-1)nw=-1
      dy2psi(ip)=fac*(-15.*psi(ip)+psi(ip+4*nw*m)-6.*psi(ip+3*m*nw)
     +                +14.*psi(ip+2*nw*m)-4.*psi(ip+nw*m)
     +                +10.*psi(ip-nw*m))
 110   continue
 100  continue
      return
      end
