c-----------------------------------------------------------------------
c grid.f  
c 
c Grid-to-Grid operations for AMR
c  - Prolongation operators
c  - Restriction operators
c-----------------------------------------------------------------------
c  Routine to convert from global array coordinates to local array coordinates
      integer function getindx(loc, lb, stride )
         implicit none
   
         integer   loc, lb, stride

         getindx = (loc - lb)/stride + 1
         return
      end
c-----------------------------------------------------------------------
c Prolongation operator: prolong2d2 (2D 2nd order)
c  USES INTEGER Lower and Upper bounds for use with DAGH
c
c  - Uses linear interpolation where necessary
c Interface:
c   shapef(2) := shape of fine grid function
c   shapec(2) := shape of coarse grid function
c
c   uf(,) := fine grid function
c   uc(,) := coarse grid function
c
c   lbc(2) := lower bound for coarse grid
c   ubc(2) := upper bound for coarse grid
c   lbf(2) := lower bound for fine grid
c   ubf(2) := upper bound for fine grid
c   lbr(2) := lower bound for region prolongation desired
c   ufr(2) := upper bound for region prolongation desired

      subroutine prolong2d2(
     .                                uc, lbc, ubc, shapec,
     .                                uf, lbf, ubf, shapef, 
     .                                lbr, ubr, args, argc)
         implicit none

         integer   shapef(2), shapec(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),
     .            getindx


                        integer argc
                        real*8 args(argc) 

c      Local variables

         integer   i, j, ic,jc,
     .            ilbf(2), iubf(2),
     .            ilbc(2), iubc(2),
     .            refine,
     .            ifine, icoarse,
     .            jfine, jcoarse, stridec, stridef

         real*8   hf, hc, rem(2), x, y, 
     .            a11, a12, a21, a22



c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c         Get strides coarse grid
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

         stridec = (ubc(1) - lbc(1))/(shapec(1)-1)
         stridef = (ubf(1) - lbf(1))/(shapef(1)-1)
         refine = stridec/stridef

c	write(45,*)'In prolong: stridec=',stridec,' stridef=',stridef
c	write(46,*)'In prolong: stridec=',
c     .               (ubc(1) - lbc(1))/(shapec(1)-1),
c     .                         ' stridef=',
c     .               (ubf(1) - lbf(1))/(shapef(1)-1)

c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c      Prolongation region is defined on the domain of the fine grid. 
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c      Loop through points in prolongation region

c	 write(55,*)shapec(1), shapec(2)
c	 write(55,*)shapef(1), shapef(2)
c        write(55,*)'stride coarse ',stridec
c        write(55,*)'stride fine ',stridef
c        write(55,*)'Region l/u x l/u y-----------------------'
c        write(55,*)lbr(1), ubr(1),lbr(2), ubr(2)
c        write(55,*)'Coarse l/u x l/u y-----------------------'
c        write(55,*)lbc(1), ubc(1),lbc(2), ubc(2)
c        write(55,*)'Fine l/u x l/u y-----------------------'
c        write(55,*)lbf(1), ubf(1),lbf(2), ubf(2)
c        write(55,*)'********************'

         do i = lbr(1), ubr(1), stridef
           do j = lbr(2), ubr(2), stridef

c             calculate coarse grid integer coordinates difference
            ic = i - lbc(1)
            jc = j - lbc(2)

            ifine = getindx(i, lbf(1), stridef)
            jfine = getindx(j, lbf(2), stridef)

            icoarse = getindx(i, lbc(1), stridec)
            jcoarse = getindx(j, lbc(2), stridec)

c            if(trace)then
c            write(6,*)icoarse, jcoarse, ifine, jfine
c            write(6,*)ic, jc, mod(ic,refine), mod(jc,refine)
c            write(6,*)'**********'
c            end if

c            if(trace)then
c            write(6,*)i,j,ic,jc
c            write(6,*)'m ',mod(i-1,refine),mod(j-1,refine)
c            write(6,*)'********************************** '
c            end if

            if(mod(ic,stridec) .eq. 0 .and. mod(jc,stridec) .eq.0)then
               uf(ifine,jfine) = uc(icoarse,jcoarse)
            end if

            if(mod(ic,stridec) .eq. 0 .and. mod(jc,stridec) .ne.0)then
               uf(ifine,jfine) = 0.5 * uc(icoarse,jcoarse) + 
     .                           0.5 * uc(icoarse,jcoarse+1)
            end if

            if(mod(ic,stridec) .ne. 0 .and. mod(jc,stridec) .eq.0)then
               uf(ifine,jfine) = 0.5 * uc(icoarse,jcoarse) + 
     .                           0.5 * uc(icoarse+1,jcoarse)
            end if

            if(mod(ic,stridec) .ne. 0 .and. mod(jc,stridec).ne. 0)then
               uf(ifine,jfine) = 0.25 * uc(icoarse,jcoarse) + 
     .                           0.25 * uc(icoarse+1,jcoarse) +
     .                           0.25 * uc(icoarse,jcoarse+1) +
     .                           0.25 * uc(icoarse+1,jcoarse+1)
            end if

           end do
         end do

         return
      end
c-----------------------------------------------------------------------
c Restriction operator: restrict2d (pure injections)
c Interface:
c   shapef(2) := shape of fine grid function
c   shapec(2) := shape of coarse grid function
c
c   uf(,) := fine grid function
c   uc(,) := coarse grid function
c
c   lbc(2) := lower bound for coarse grid
c   ubc(2) := upper bound for coarse grid
c   lbf(2) := lower bound for fine grid
c   ubf(2) := upper bound for fine grid
c   lbr(2) := lower bound for region restriction desired
c   ufr(2) := upper bound for region restriction desired

      subroutine restrict2d2(uf, lbf, ubf, shapef, 
     .                                 uc, lbc, ubc, shapec,
     .                                 lbr, ubr, args, argc)
         implicit none

         integer   shapef(2), shapec(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
         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...
         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. jcoarse .gt. shapec(2))then
		write(0,*)'ERROR in restriction: ',icoarse,jcoarse
           end if
               uc(icoarse,jcoarse) = uf(ifine,jfine)

            end do
         end do


         return
      end
c-----------------------------------------------------------------------