c    @(#)arelvor.f	1.1
c
      subroutine arelvor(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(1),dx2psi(1),dy2psi(1)
c
       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
      nn=1
      if(i.eq.m)nn=-1
      dx2psi(ip)=fac*(45.*psi(ip)-10.*psi(ip+nn*5)+61.*psi(ip+nn*4)
     +                -156.*psi(ip+3*nn)+214.*psi(ip+2*nn)
     +                -154.*psi(ip+nn))
       go to 70
 60   continue
      nn=1
      if(i.eq.m-1)nn=-1
      dx2psi(ip)=fac*(-15.*psi(ip)+psi(ip+4*nn)-6.*psi(ip+nn*3)
     +                 +14.*psi(ip+2*nn)-4.*psi(ip+nn)
     +                 +10.*psi(ip-nn))
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
      nn=1
      if(j.eq.n)nn=-1
      dy2psi(ip)=fac*(45*psi(ip)-10.*psi(ip+5*nn*m)+61.*psi(ip+4*nn*m)
     +                -156.*psi(ip+3*nn*m)+214.*psi(ip+2*m*nn)
     +                -154.*psi(ip+nn*m))
      go to 110
 90   continue
      nn=1
      if(j.eq.n-1)nn=-1
      dy2psi(ip)=fac*(-15.*psi(ip)+psi(ip+4*nn*m)-6.*psi(ip+3*m*nn)
     +                +14.*psi(ip+2*nn*m)-4.*psi(ip+nn*m)
     +                +10.*psi(ip-nn*m))
 110   continue
 100  continue
      return
      end
