CCC# c/*@@
CCC# c  @file prolong2.F
CCC# c  @author Paul Walker
CCC# c  @date Halloween 1996!
CCC# c  @desc
CCC# c  The new 2D prolongation operator based on the idea of calling
CCC# c  sweeps of the 1D interpolator with temporaries. 
CCC# c  @enddesc
CCC# c  @rcsid $Id$
CCC# c@@*/
CCC# 
CCC# c/*@@
CCC# c  @routine prolong2
CCC# c  @author Paul Walker
CCC# c  @date Halloween 1996!
CCC# c  @desc The new 2D prolongation operator
CCC# c  This operator is based on the idea that 2D interpolation
CCC# c  can be done as sweeps of a 1D interpolator in one direction
CCC# c  followed by sweeps again in the second direction. This
CCC# c  philosophy allows us very easy access to high order
CCC# c  interpolation, since we only need to write it in 1-D once
CCC# c  and then it is black-boxed external in higher dimensions.
CCC# c  <p>
CCC# c  This requires fortran 90. No getting out of it. Sorry!
CCC# c  @enddesc
CCC# c  @history
CCC# c  
CCC# c  @endhistory
CCC# c @@*/

      subroutine prolong2(datac, lbc, ubc, shc, 
     $     dataf, lbf, ubf, shf,
     $     lbr, ubr, shr,
     $     argv, argc)

      implicit none
      integer lbc(2),ubc(2),shc(2),lbf(2),ubf(2),shf(2)
      integer lbr(2),ubr(2),shr(2), argc
      real*8 argv(argc)
      real*8 datac(shc(1),shc(2)), dataf(shf(1),shf(2))
      
c     OK, so the plan of attack is to make temporaries in
c     each direction and use them with a 1D interpolator.
      real*8 tmpline(shf(1),shf(2))
      real*8 tmpc1(shc(1)), tmpc2(shc(2))
      real*8 tmpf1(shf(1)), tmpf2(shf(2))
      integer lbct(1), ubct(1), shct(1)
      integer lbft(1), ubft(1), shft(1)
      integer lbrt(1), ubrt(1), shrt(1)
      
c     deltas so we can locate ourselves
      integer dif(2), dic(2)
      
c     Starting points
      integer startf(2), endf(2), starttmp, endtmp
      integer nlines, cbelow, cabove
        
c     Loopers and coordinates
      integer i, j, ic, jc, lfi, lfj
        
c     Get the nearest index BELOW or ON our position in the coarse
c     grid (but still within the bbox)
      integer getcoarseidx
        
c     Set up the deltas
      do i=1,2
         dif(i) = (ubf(i) - lbf(i)) / (shf(i) - 1)
         dic(i) = (ubc(i) - lbc(i)) / (shc(i) - 1)
      enddo

c Set up the starting and ending points
      do i=1,2
         startf(i) = (lbr(i) - lbf(i))/dif(i)+1
         endf(i) = (ubr(i) - lbf(i))/dif(i)+1
      enddo

c     Handle the case of one of the dimensions is 1 as a special
c     case for efficiency sake. But do this later...
      if (shr(1) .eq. 1) then
c     We can only be more efficient if it is coincident with a coarse
c     grid line. Leave this...
      endif

      if (shr(2) .eq. 1) then
c     We can only be more efficient if it is coincident with a coarse
c     grid line. Leave this...
      endif
      
c     Now locate the coarse grid line directly below and above our 
c     array. These become our starting points.
      cbelow = lbf(1)-mod(lbf(1),dic(1))
      if (cbelow .lt. lbf(1)) then
         cbelow = lbf(1)
      endif
      cabove = ubf(1)+dic(1)-mod(ubf(1),dic(1))
      if (cabove .gt.ubf(1)) then
         cabove = ubf(1)
      endif
      
      nlines = (cabove-cbelow)/(dic(1)) + 1
      if (nlines .gt. shc(1)) then
         write (*,*) "ERROR: Too many lines!"
      endif
      
c     OK, loop over the i index interpolating in the j dir
c     Set up the temporary bounding boxes outside the loop.
      lbct(1) = lbc(2)
      ubct(1) = ubc(2)
      shct(1) = shc(2)
      lbft(1) = lbf(2)
      ubft(1) = ubf(2)
      shft(1) = shf(2)
      lbrt(1) = lbf(2)
      ubrt(1) = ubf(2)
      shrt(1) = shf(2)
      
c     Loop over the WHOLE array on the coarse-fine intersection
      do i=1,nlines
         lfi = (i-1)*dic(1)+cbelow
         ic = getcoarseidx(lfi,lbc(1),ubc(1),dic(1))
c     OK, so use f90 notation to grab the temporaries
         if (ic .lt. 1 .or. ic .gt. shc(1)) then 
            write (*,*) "FAILED BOUND IN PROLONG2 1",ic,nlines
         endif
         tmpc2 = datac(ic,:)
c     Call the line interpolater
         call lineInterp(tmpc2,lbct,ubct,shct,
     $        tmpf2, lbft, ubft, shft,
     $        lbrt, ubrt, shrt)
c Copy the data back into the new array. Be careful with 
c bounds here (later)
         do j=1,shf(2)
            if (ic .lt. 1 .or. ic .gt. shf(1)) then 
               write (*,*) "FAILED BOUND IN PROLONG2 2",i,nlines
            endif
            tmpline(i,j) = tmpf2(j)
         enddo
      enddo

c OK, now loop over the whole y direction  filling in with
c the pre-interpolated temporaries. One thing here, though,
c is that there are some points which could use the coarse
c grid, and therefore get more boundary points. We\'ll put
c that in later if needs be.

c Note this CHANGES since we pass a temporary which is coincident
c with the fine array with lower resolution.
      lbct(1) = cbelow
      ubct(1) = cabove
      shct(1) = nlines
c     These are the same, of course
      lbft(1) = lbf(1)
      ubft(1) = ubf(1)
      shft(1) = shf(1)
      lbrt(1) = lbr(1)
      ubrt(1) = ubr(1)
      shrt(1) = shr(1)

      do j=startf(2),endf(2)
c OK, so use f90 notation to grab the temporaries. This will
c have the wrong size at the target, note, but is is
c guaranteed to be too BIG since the coarse is at least
c bigger than the fine in integer space.
         tmpc1(1:nlines) = tmpline(1:nlines,j)
c Call the line interpolater
         call lineInterp(tmpc1,lbct,ubct,shct,
     $        tmpf1, lbft, ubft, shft,
     $        lbrt, ubrt, shrt)
c Copy the data back into the new array. Also note bounds!
         do i=startf(1),endf(1)
              dataf(i,j) = tmpf1(i)
           enddo
        enddo
        
c     And that should be the complete interpolation
        return
        end
      
c/*@@
c  @routine getcoarseidx
c  @author Paul Walker
c  @date Haloween 1996
c  @desc
c  Returns the coarse grid index closest to (but lower than)
c  the fine grid coordinate we pass in
c  @enddesc
c@@*/
      integer function getcoarseidx(lfi,lbc,ubc,dic)
      implicit none
      integer lfi, lbc, ubc, dic
      integer check		
      getcoarseidx = (lfi - lbc)/dic + 1			
     
      return
      end