#define geti(a) 1
CCC#  /*@@
CCC#    @file      prolong3.F
CCC#    @date      Tue Jan 30 11:47:12 1996
CCC#    @author    Paul Walker
CCC#    @desc 
CCC#       DAGH wrapper for Karen's data parallel interpolator
CCC#       which will make these useable as dagh prolongation functions.
CCC#       <p>
CCC#    @enddesc 
CCC#  @@*/
CCC# 
CCC#  /*@@
CCC#    @routine    prolong3
CCC#    @date       Tue Jan 30 11:48:31 1996
CCC#    @author     Paul Walker
CCC#    @desc 
CCC#       This is the fortran-interface to the dagh prolong routine as 
CCC#       described on p. 16 of the DAGH Users Ref (draft).
CCC#    @enddesc 
CCC# 
CCC#    @calledby   RefineSystem, TruncByThresh
CCC#    @calls      mymin
CCC#    @history 
CCC#  
CCC#    @endhistory 
CCC#    @comment 
CCC#       Note that even though this is "calledby" @seeroutien RefineSystem
CCC#       that is not really true.  @seeroutine main sets this up as the
CCC#       prolong function with a call to the AMR library, and then it is 
CCC#       called elsewhere by DAGH internals, but as Prolong in C++ not as 
CCC#       f_prolong.
CCC#    @endcomment
CCC# 
CCC#    @var     DAGH_ARGS
CCC#    @vdesc   DAGH defined args
CCC#    @vtype   misc
CCC#    @vio     inout
CCC#    @vcomment 
CCC#       Basically x_from is the coarse grid, x_to is the fine grid, and x_r is the 
CCC#       region on which we should interpolate
CCC#    @endvar 
CCC# 
CCC#    @par     prolong_order
CCC#    @pdesc   interpolation order for prolongation
CCC#    @ptype   int
CCC#    @pvalues 1,2
CCC#    @pcomment 
CCC#       Linear or quadratic.  For more details, ask karen!
CCC#    @endpar 
CCC# 
CCC# 
CCC# @@*/

      subroutine prolong3(data_from, lb_from, ub_from, 
     $  shape_from,
     $  data_to, lb_to, ub_to, shape_to,
     $  lb_r, ub_r, shape_r)
      implicit none

c DAGH Variables
      integer lb_from(3), ub_from(3), shape_from(3)
      real*8 data_from(shape_from(1), shape_from(2), shape_from(3))
      integer lb_to(3), ub_to(3), shape_to(3)
      real*8 data_to(shape_to(1), shape_to(2), shape_to(3))
      integer lb_r(3), ub_r(3), shape_r(3)

c Local vars
      real*8 tmp(2,2), tmp2(2), frac(3)
      integer dit(3), dif(3), ts(3), te(3)
      integer ct(3), cf(3), idxf(3)
      integer i,j,k

c Orient ourselves
      do i=1,3
         dit(i) = (ub_to(i)  -lb_to(i))  / (shape_to(i)-1)
         dif(i) = (ub_from(i)-lb_from(i))/ (shape_from(i)-1)
         
         ts(i) = (lb_r(i)-lb_to(i))/dit(i) + 1
         te(i) = (ub_r(i)-lb_to(i))/dit(i) + 1
      enddo

      if (geti("prolong_order") .eq. 2) then
c         write (*,*) "WARNING: Write this interpolator!"
      endif

c This is based heavily on the 1d interpolator
      if (geti("prolong_order") .eq. 1 .or. .true.) then
         do i=ts(1),te(1)
            ct(1) = (i-1)*dit(1) + lb_to(1)
            idxf(1) = (ct(1)-lb_from(1))/dif(1) + 1
            if (idxf(1)+1 .gt. shape_from(1)) then 
               idxf(1) = idxf(1)-1
            endif
            cf(1) = (idxf(1)-1)*dif(1) + lb_from(1)
            frac(1) = 1.0D0*(ct(1)-cf(1))/(1.0D0*dif(1))
            do j=ts(2),te(2)
               ct(2) = (j-1)*dit(2) + lb_to(2)
               idxf(2) = (ct(2)-lb_from(2))/dif(2) + 1
               if (idxf(2)+1 .gt. shape_from(2)) then 
                  idxf(2) = idxf(2)-1
               endif
               cf(2) = (idxf(2)-1)*dif(2) + lb_from(2)
               frac(2) = 1.0D0*(ct(2)-cf(2))/(1.0D0*dif(2))
               do k=ts(3),te(3)
                  ct(3) = (k-1)*dit(3) + lb_to(3)
                  idxf(3) = (ct(3)-lb_from(3))/dif(3) + 1
                  if (idxf(3)+1 .gt. shape_from(3)) then 
                     idxf(3) = idxf(3)-1
                  endif
                  cf(3) = (idxf(3)-1)*dif(3) + lb_from(3)
                  frac(3) = 1.0D0*(ct(3)-cf(3))/(1.0D0*dif(3))

c     Actually do the interpolation
                  tmp(1,1) = (1.0-frac(2))*data_from(idxf(1),idxf(2),
     $                  idxf(3)) +
     $                  frac(2)*data_from(idxf(1),idxf(2)+1,idxf(3))
                  tmp(2,1) = (1.0-frac(2))*data_from(idxf(1)+1,idxf(2),
     $                  idxf(3)) +
     $                  frac(2)*data_from(idxf(1)+1,idxf(2)+1,idxf(3))
                  tmp(1,2) = (1.0-frac(2))*data_from(idxf(1),idxf(2),
     $                  idxf(3)) + 
     $                  frac(2)*data_from(idxf(1),idxf(2)+1,idxf(3))
                  tmp(2,2) = (1.0-frac(2))*data_from(idxf(1)+1,idxf(2),
     $                  idxf(3)) + 
     $                  frac(2)*data_from(idxf(1)+1,idxf(2)+1,idxf(3))

                  tmp2(1) = (1.0-frac(1))*tmp(1,1) +
     $                 frac(1)*tmp(2,1)
                  tmp2(2) = (1.0-frac(1))*tmp(1,2) +
     $                 frac(1)*tmp(2,2)
                  
                  data_to(i,j,k) = (1.0-frac(3))*tmp2(1) +
     $                 frac(3)*tmp2(2)
               enddo
            enddo
         enddo
      endif

      return
      end

     
CCC#  /*@@
CCC#    @routine    mymin
CCC#    @date       Fri May 24 13:25:28 1996
CCC#    @author     Paul Walker
CCC#    @desc 
CCC#       A replacement for min, which appears to be broken on
CCC#       AIX.
CCC#    @enddesc 
CCC#  @@*/

      function mymin(a,b)
      real*8 a
      integer b
      integer mymin
      if (a .lt. b*1.0) then
         mymin = int(a)
      else
         mymin = b
      endif
      return
      end