#define geti(a) 1
CCC#  /*@@
CCC#    @file      restrict2.F
CCC#    @date      Wed Oct 30 11:02:13 1996
CCC#    @author    Paul Walker
CCC#    @desc 
CCC#       A 2-dimensional restrictor, based on the one which
CCC#       was written by mijan, but with some cleanups and improvements
CCC#    @enddesc 
CCC#  @@*/
CCC# 
CCC#  /*@@
CCC#    @routine    restrict2
CCC#    @date       Wed Oct 30 11:02:32 1996
CCC#    @author     Paul Walker
CCC#    @desc 
CCC#       A 2-dimensional restrictor which does direct
CCC#       injection or full-weighted averaging.
CCC#       <p>
CCC#       The variable list is the standard DAGH interface,
CCC#       so it is not documented here.
CCC#    @enddesc 
CCC#    @calledby   RecursiveAMRStep
CCC#    @history 
CCC#    @hdate Wed Oct 30 11:03:02 1996 @hauthor Paul Walker
CCC#    @hdesc Cleaned up, put in grdoc and restrict order. 
CCC#    @endhistory 
CCC# 
CCC#    @par     restrict_order
CCC#    @pdesc   What sort of restriction to do
CCC#    @ptype   int
CCC#    @pvalues 0,1
CCC#    @pcomment 
CCC#       0 = direct, 1 = full_weighting.
CCC#    @endpar 
CCC#  @@*/

cc#define DEBUG

      subroutine restrict2(uf, lbf, ubf, shapef, 
     $     uc, lbc, ubc, shapec,
     $     lbr, ubr, shaper, args, argc)
      implicit none
      
      integer   shapef(2), shapec(2), shaper(2)
      
      real*8   uf(shapef(1), shapef(2)), 
     $     uc(shapec(1), shapec(2))
      
      integer  lbf(2), ubf(2), 
     $     lbc(2), ubc(2),
     $     lbr(2), ubr(2)
      
      integer argc
      real*8 args(argc) 
         
c      Local variables

      integer   i, j, imin, imax, jmin, jmax,
     .     ifine, icoarse,
     .     jfine, jcoarse, refine, stridec, stridef,
     .     ilbc(2), iubc(2), getindx 
      
      real*8   hf, hc, x, y

      
      stridec = (ubc(1) - lbc(1))/(shapec(1) - 1)
      stridef = (ubf(1) - lbf(1))/(shapef(1) - 1)
CC#ifdef DEBUG
CC      write (*,*) "STRIDES ",stridec, stridef
CC      write (*,*) "BOUND C", ubc(1),lbc(1),shapec(1) 
CC#endif
      refine = stridec/stridef
      

c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c     Find coarse domain over which to refine
c     Take three regions and select out intersection

      imin = max(lbf(1), lbc(1), lbr(1))
      imax = min(ubf(1), ubc(1), ubr(1))

      jmin = max(lbf(2), lbc(2), lbr(2))
      jmax = min(ubf(2), ubc(2), ubr(2))
      
      if (mod(imin-lbc(1),stridec) .ne. 0) then
         imin = imin + mod(imin-lbc(1),stridec) 
      endif
      if (mod(jmin-lbc(2),stridec) .ne. 0) then
         jmin = jmin + mod(jmin-lbc(2),stridec) 
      endif
      
c      Inject points to coarse grid from fine grid
c     Loop from lower bound to upper bound with stride of refine.
c     Convert the integer coordinates to fine and coarse grid absolute
c     coordinates...

C I know this replicates code, but it is faster since it takes the iff out of the
c inner loop.
      if (geti("restrict_order") .eq. 0) then 
CC#ifdef DEBUG
CC         write (*,*) "Is : ",imin,imax,jmin,jmax,stridec
CC         write (*,*) "Sz : ",shapec,shapef
CC#endif
         do i = imin, imax,stridec
            do j = jmin, jmax,stridec
               
               ifine = getindx(i, lbf(1), stridef)
               jfine = getindx(j, lbf(2), stridef)
               
               icoarse = getindx(i, lbc(1), stridec)
               jcoarse = getindx(j, lbc(2), stridec)
CC#ifdef DEBUG
CC               if (icoarse .gt. shapec(1) .or.
CC     $              icoarse .lt. 1 .or.
CC     $              jcoarse .gt. shapec(2) .or.
CC     $              icoarse .lt. 1) then
CC                  write (*,*) "AWFUL BOUNDS ERROR 1 ",i,j,icoarse,jcoarse 
CC               endif 
CC#endif
               uc(icoarse,jcoarse) = uf(ifine,jfine)
            enddo
         enddo
      else
c This could be made more efifcient by doing the boundaries first,
c then the inner regions. If we use this in production it should be
c changed.
         do i = imin, imax,stridec
            do j = jmin, jmax,stridec
               
               ifine = getindx(i, lbf(1), stridef)
               jfine = getindx(j, lbf(2), stridef)
               
               icoarse = getindx(i, lbc(1), stridec)
               jcoarse = getindx(j, lbc(2), stridec) 
               if (icoarse .gt. shapec(1) .or.
     $              icoarse .lt. 1 .or.
     $              jcoarse .gt. shapec(2) .or.
     $              icoarse .lt. 1) then
                  write (*,*) "AWFUL BOUNDS ERROR 2 ",i,j 
               endif
               if ((ifine .eq. 1) .or.
     $              (ifine .eq. shapef(1)) .or.
     $              (jfine .eq. 1) .or.
     $              (jfine .eq. shapef(2))) then 
                  uc(icoarse,jcoarse) = uf(ifine,jfine)
               else
                  uc(icoarse,jcoarse) = 0.25 * uf(ifine,jfine) +
     $                 0.125 * (uf(ifine+1,jfine) +
     $                 uf(ifine-1,jfine) +
     $                 uf(ifine,jfine+1) +
     $                 uf(ifine,jfine-1)) +
     $                 0.0625 * (uf(ifine+1,jfine+1) +
     $                 uf(ifine+1,jfine-1) +
     $                 uf(ifine-1,jfine+1) +
     $                 uf(ifine-1,jfine-1))
               endif
            end do
         end do
      endif
      
      return
      end
c-----------------------------------------------------------------------


CCC# /*@@
CCC#    @routine    getindx
CCC#    @date       Mon Jun 17 12:27:29 1996
CCC#    @author     Mijan Huq
CCC#    @desc 
CCC#       Figures out indices.
CCC#    @enddesc 
CCC#    @calls     
CCC#    @calledby  restrict2
CCC# @@*/

      integer function getindx(loc, lb, stride )
      implicit none
      
      integer   loc, lb, stride
      
      getindx = (loc - lb)/stride + 1
      return
      end