c
c @(#)agrad.f	1.1
c
      subroutine agrad(r,dxr,dyr,m,n,hy)
c
c --- this routine calculates fourth order derivatives
c --- at the boundaries and at the interior
c
c
      dimension r(m*n),dxr(m*n),dyr(m*n)
c
       fac=-1./12./hy
      ip=0
      do 25 j=1,n
      do 20 i=1,m
      ip=ip+1
      if(i.eq.1.or.i.eq.m) go to 5
      if(i.eq.2.or.i.eq.m-1) go to 6
      dxr(ip)=-fac*(-r(ip+2)+8*r(ip+1)-8*r(ip-1)+ r(ip-2))
      go to 10
c
   5  continue
      nw=1
      if(i.eq.m) nw=-1
      dxr(ip)=fac*nw*(25.*r(ip)+3.*r(ip+nw*4)-16.*r(ip+3*nw)+
     +               36.*r(ip+2*nw)-48.*r(ip+nw))
c
      go to 10
 6    continue
      nw=1
      if(i.eq.m-1)nw=-1
      dxr(ip)=fac*nw*(10.*r(ip)-r(ip+3*nw)+6.*r(ip+2*nw)-18.*r(ip+nw)
     +               +3.*r(ip-nw))
c
   10 continue
      if(j.eq.1.or.j.eq.n) go to 15
      if(j.eq.2.or.j.eq.n-1) go to 16
      dyr(ip)=-fac*(-r(ip+2*m)+8.*r(ip+m)-8.*r(ip-m)+r(ip-2*m))
      go to 20
c
   15 continue
      nw=1
      if(j.eq.n) nw=-1
      dyr(ip)=fac*nw*(25.*r(ip)+3.*r(ip+4*m*nw)-16.*r(ip+3*m*nw)
     +               +36.*r(ip+2*m*nw)-48.*r(ip+m*nw))
      go to 20
c
 16   continue
       nw=1
       if(j.eq.n-1) nw=-1
       dyr(ip)=fac*nw*(10.*r(ip)-r(ip+3*m*nw)+6.*r(ip+2*m*nw)-
     +                     18.*r(ip+nw*m)+3.*r(ip-nw*m))
   20 continue
   25 continue
      return
      end
