c*******************************************************************************
c* Fine to coarse grid interpolation using 4X4 bi-cubic interpolation.
c*     fine: The fine grid.
c*   coarse: The coarse grid.
c*  fext(4): The "extent" of the fine grid,
c*           fine(fext(1):fext(3),fext(2):fext(4))
c*  ffil(4): The "fill" of the fine grid.  The range of indices for which
c*           the fine grid containts valid points for use in interpolation.
c*  cint(4): The "interior" of the coarse grid.  The range of indices which
c*           are to be filled by the interpolation.
c*  cext(4): The "extent" of the coarse grid,
c*           coarse(cext(1):cext(3),cext(2):cext(4))
c*******************************************************************************
c* Implimentation notes:
c*    ffil must always be a superset of cint, ie.
c*          ffil(1) <= 2*cint(1) - 1,
c*          ffil(2) <= 2*cint(2) - 1,
c*          ffil(3) >= 2*cint(3)    ,
c*          ffil(4) >= 2*cint(4)    .
c*******************************************************************************
      subroutine ftoc4x4(fine,coarse,fext,ffil,cint,cext)
      implicit none
      integer fext(4),ffil(4),cint(4),cext(4)
      real*8 fine(fext(1):fext(3),fext(2):fext(4))
      real*8 coarse(cext(1):cext(3),cext(2):cext(4))
      integer c1i,c1f,c2i,c2f
      integer f1i,f1f,f2i,f2f
      integer c1,c2,f10,f11,f12,f13,f20,f21,f22,f23
      real*8 F225o256, F135o256, F81o256, F75o256, F45o256,
     &       F25o256,  F15o256,  F9o256,  F5o256,  F1o256
      parameter( F225o256 = 225.D0/256.D0, F135o256 = 135.D0/256.D0,
     &            F81o256 =  81.D0/256.D0,  F75o256 =  75.D0/256.D0,
     &            F45o256 =  45.D0/256.D0,  F25o256 =  25.D0/256.D0,
     &            F15o256 =  15.D0/256.D0,   F9o256 =   9.D0/256.D0,
     &             F5o256 =   5.D0/256.D0,   F1o256 =   1.D0/256.D0 )
c
c Make sure the "fill" region of the fine grid is at least 4X4.
      if(((ffil(3) - ffil(1)).lt.3).or.((ffil(4) - ffil(2)).lt.3)) then
        write(*,*) "ftoc4x4: fill region smaller than 4X4!"
        stop
      end if
c
c Set the coarse interior which can be filled by "centered" 4X4 interpolation.
		c1i = cint(1)
		if(ffil(1).eq.(2*c1i-1)) c1i = c1i + 1
c     c1i = MAX(cint(1),(ffil(1) + 1)/2 + 1)
      f1i = 2*c1i
		c2i = cint(2)
		if(ffil(2).eq.(2*c2i-1)) c2i = c2i + 1
c     c2i = MAX(cint(2),(ffil(2) + 1)/2 + 1)
      f2i = 2*c2i
		c1f = cint(3)
		if(ffil(3).eq.(2*c1f)) c1f = c1f - 1
c     c1f = MIN(cint(3),(ffil(3) - 1)/2)
      f1f = 2*c1f
		c2f = cint(4)
		if(ffil(4).eq.(2*c2f)) c2f = c2f - 1
c     c2f = MIN(cint(4),(ffil(4) - 1)/2)
      f2f = 2*c2f
c-----------------------------------------------------------------------
c F90 version
c     coarse(c1i:c1f,c2i:c2f) =
c    &          F81o256*( fine( f1i-1: f1f-1: 2, f2i-1: f2f-1: 2) +
c    &                    fine( f1i  : f1f  : 2, f2i-1: f2f-1: 2) +
c    &                    fine( f1i-1: f1f-1: 2, f2i  : f2f  : 2) +
c    &                    fine( f1i  : f1f  : 2, f2i  : f2f  : 2)  )
c    &         - F9o256*( fine( f1i-1: f1f-1: 2, f2i-2: f2f-2: 2) +
c    &                    fine( f1i  : f1f  : 2, f2i-2: f2f-2: 2) +
c    &                    fine( f1i-2: f1f-2: 2, f2i-1: f2f-1: 2) +
c    &                    fine( f1i+1: f1f+1: 2, f2i-1: f2f-1: 2) +
c    &                    fine( f1i-2: f1f-2: 2, f2i  : f2f  : 2) +
c    &                    fine( f1i+1: f1f+1: 2, f2i  : f2f  : 2) +
c    &                    fine( f1i-1: f1f-1: 2, f2i+1: f2f+1: 2) +
c    &                    fine( f1i  : f1f  : 2, f2i+1: f2f+1: 2)  )
c    &         + F1o256*( fine( f1i-2: f1f-2: 2, f2i-2: f2f-2: 2) +
c    &                    fine( f1i+1: f1f+1: 2, f2i-2: f2f-2: 2) +
c    &                    fine( f1i-2: f1f-2: 2, f2i+1: f2f+1: 2) +
c    &                    fine( f1i+1: f1f+1: 2, f2i+1: f2f+1: 2)  )
c-----------------------------------------------------------------------
c F77 version
      f21 = f2i
      do c2 = c2i,c2f
        f11 = f1i
        do c1 = c1i,c1f
          coarse(c1,c2) =
     &          F81o256*( fine( f11-1, f21-1) +
     &                    fine( f11  , f21-1) +
     &                    fine( f11-1, f21  ) +
     &                    fine( f11  , f21  )  )
     &         - F9o256*( fine( f11-1, f21-2) +
     &                    fine( f11  , f21-2) +
     &                    fine( f11-2, f21-1) +
     &                    fine( f11+1, f21-1) +
     &                    fine( f11-2, f21  ) +
     &                    fine( f11+1, f21  ) +
     &                    fine( f11-1, f21+1) +
     &                    fine( f11  , f21+1)  )
     &         + F1o256*( fine( f11-2, f21-2) +
     &                    fine( f11+1, f21-2) +
     &                    fine( f11-2, f21+1) +
     &                    fine( f11+1, f21+1)  )
          f11 = f11 + 2
        end do
        f21 = f21 + 2
      end do
c-----------------------------------------------------------------------
c
c Set cint(1) boundary
      if(c1i.gt.cint(1)) then
        c1 = cint(1)
        f10 = 2*c1 -1
        f11 = f10 + 1
        f12 = f11 + 1
        f13 = f12 + 1
c-----------------------------------------------------------------------
c F90 version
c       coarse(c1,c2i:c2f) =
c    &         F135o256*( fine( f11, f2i-1: f2f-1: 2) +
c    &                    fine( f11, f2i  : f2f  : 2)  )
c    &        + F45o256*( fine( f10, f2i-1: f2f-1: 2) +
c    &                    fine( f10, f2i  : f2f  : 2) -
c    &                    fine( f12, f2i-1: f2f-1: 2) -
c    &                    fine( f12, f2i  : f2f  : 2)  )
c    &        - F15o256*( fine( f11, f2i-2: f2f-2: 2) +
c    &                    fine( f11, f2i+1: f2f+1: 2)  )
c    &         + F9o256*( fine( f13, f2i-1: f2f-1: 2) +
c    &                    fine( f13, f2i  : f2f  : 2)  )
c    &         - F5o256*( fine( f10, f2i-2: f2f-2: 2) +
c    &                    fine( f10, f2i+1: f2f+1: 2) -
c    &                    fine( f12, f2i-2: f2f-2: 2) -
c    &                    fine( f12, f2i+1: f2f+1: 2)  )
c    &         - F1o256*( fine( f13, f2i-2: f2f-2: 2) +
c    &                    fine( f13, f2i+1: f2f+1: 2)  )
c-----------------------------------------------------------------------
c F77 version
        f21 = f2i
        do c2 = c2i,c2f
          coarse(c1,c2) =
     &         F135o256*( fine( f11, f21-1) +
     &                    fine( f11, f21  )  )
     &        + F45o256*( fine( f10, f21-1) +
     &                    fine( f10, f21  ) -
     &                    fine( f12, f21-1) -
     &                    fine( f12, f21  )  )
     &        - F15o256*( fine( f11, f21-2) +
     &                    fine( f11, f21+1)  )
     &         + F9o256*( fine( f13, f21-1) +
     &                    fine( f13, f21  )  )
     &         - F5o256*( fine( f10, f21-2) +
     &                    fine( f10, f21+1) -
     &                    fine( f12, f21-2) -
     &                    fine( f12, f21+1)  )
     &         - F1o256*( fine( f13, f21-2) +
     &                    fine( f13, f21+1)  )
          f21 = f21 + 2
        end do
c-----------------------------------------------------------------------
c
c Set cint(1)-cint(4) corner
        if(c2f.lt.cint(4)) then
          c2 = cint(4)
          f20 = 2*c2
          f21 = f20 - 1
          f22 = f21 - 1
          f23 = f22 - 1
          coarse(c1,c2) =
     &         F225o256*( fine( f11, f21)  )
     &       +  F75o256*( fine( f11, f20) +
     &                    fine( f10, f21) -
     &                    fine( f12, f21) -
     &                    fine( f11, f22)  )
     &       +  F25o256*( fine( f10, f20) -
     &                    fine( f12, f20) -
     &                    fine( f10, f22) +
     &                    fine( f12, f22)  )
     &       +  F15o256*( fine( f13, f21) +
     &                    fine( f11, f23)  )
     &       +   F5o256*( fine( f13, f20) -
     &                    fine( f13, f22) +
     &                    fine( f10, f23) -
     &                    fine( f12, f23)  )
     &       +   F1o256*( fine( f13, f23)  )
        end if
      end if
c
c Set cint(2) boundary
      if(c2i.gt.cint(2)) then
        c2 = cint(2)
        f20 = 2*c2 -1
        f21 = f20 + 1
        f22 = f21 + 1
        f23 = f22 + 1
c
c Set cint(2)-cint(1) corner
        if(c1i.gt.cint(1)) then
          c1 = cint(1)
          f10 = 2*c1 -1
          f11 = f10 + 1
          f12 = f11 + 1
          f13 = f12 + 1
          coarse(c1,c2) =
     &         F225o256*( fine( f11, f21)  )
     &       +  F75o256*( fine( f11, f20) +
     &                    fine( f10, f21) -
     &                    fine( f12, f21) -
     &                    fine( f11, f22)  )
     &       +  F25o256*( fine( f10, f20) -
     &                    fine( f12, f20) -
     &                    fine( f10, f22) +
     &                    fine( f12, f22)  )
     &       +  F15o256*( fine( f13, f21) +
     &                    fine( f11, f23)  )
     &       +   F5o256*( fine( f13, f20) -
     &                    fine( f13, f22) +
     &                    fine( f10, f23) -
     &                    fine( f12, f23)  )
     &       +   F1o256*( fine( f13, f23)  )
        end if
c-----------------------------------------------------------------------
c F90 version
c       coarse(c1i:c1f,c2) =
c    &         F135o256*( fine( f1i-1: f1f-1: 2, f21) +
c    &                    fine( f1i  : f1f  : 2, f21)  )
c    &        + F45o256*( fine( f1i-1: f1f-1: 2, f20) +
c    &                    fine( f1i  : f1f  : 2, f20) -
c    &                    fine( f1i-1: f1f-1: 2, f22) -
c    &                    fine( f1i  : f1f  : 2, f22)  )
c    &        - F15o256*( fine( f1i-2: f1f-2: 2, f21) +
c    &                    fine( f1i+1: f1f+1: 2, f21)  )
c    &         + F9o256*( fine( f1i-1: f1f-1: 2, f23) +
c    &                    fine( f1i  : f1f  : 2, f23)  )
c    &         - F5o256*( fine( f1i-2: f1f-2: 2, f20) +
c    &                    fine( f1i+1: f1f+1: 2, f20) -
c    &                    fine( f1i-2: f1f-2: 2, f22) -
c    &                    fine( f1i+1: f1f+1: 2, f22)  )
c    &         - F1o256*( fine( f1i-2: f1f-2: 2, f23) +
c    &                    fine( f1i+1: f1f+1: 2, f23)  )
c-----------------------------------------------------------------------
c F77 version
        f11 = f1i
        do c1 = c1i,c1f
          coarse(c1,c2) =
     &         F135o256*( fine( f11-1, f21) +
     &                    fine( f11  , f21)  )
     &        + F45o256*( fine( f11-1, f20) +
     &                    fine( f11  , f20) -
     &                    fine( f11-1, f22) -
     &                    fine( f11  , f22)  )
     &        - F15o256*( fine( f11-2, f21) +
     &                    fine( f11+1, f21)  )
     &         + F9o256*( fine( f11-1, f23) +
     &                    fine( f11  , f23)  )
     &         - F5o256*( fine( f11-2, f20) +
     &                    fine( f11+1, f20) -
     &                    fine( f11-2, f22) -
     &                    fine( f11+1, f22)  )
     &         - F1o256*( fine( f11-2, f23) +
     &                    fine( f11+1, f23)  )
          f11 = f11 + 2
        end do
c-----------------------------------------------------------------------
      end if
c
c Set cint(3) boundary
      if(c1f.lt.cint(3)) then
        c1 = cint(3)
        f10 = 2*c1
        f11 = f10 - 1
        f12 = f11 - 1
        f13 = f12 - 1
c
c Set cint(3)-cint(2) corner
        if(c2i.gt.cint(2)) then
          c2 = cint(2)
          f20 = 2*c2 -1
          f21 = f20 + 1
          f22 = f21 + 1
          f23 = f22 + 1
          coarse(c1,c2) =
     &         F225o256*( fine( f11, f21)  )
     &       +  F75o256*( fine( f11, f20) +
     &                    fine( f10, f21) -
     &                    fine( f12, f21) -
     &                    fine( f11, f22)  )
     &       +  F25o256*( fine( f10, f20) -
     &                    fine( f12, f20) -
     &                    fine( f10, f22) +
     &                    fine( f12, f22)  )
     &       +  F15o256*( fine( f13, f21) +
     &                    fine( f11, f23)  )
     &       +   F5o256*( fine( f13, f20) -
     &                    fine( f13, f22) +
     &                    fine( f10, f23) -
     &                    fine( f12, f23)  )
     &       +   F1o256*( fine( f13, f23)  )
        end if
c-----------------------------------------------------------------------
c F90 version
c       coarse(c1,c2i:c2f) =
c    &         F135o256*( fine( f11, f2i-1: f2f-1: 2) +
c    &                    fine( f11, f2i  : f2f  : 2)  )
c    &        + F45o256*( fine( f10, f2i-1: f2f-1: 2) +
c    &                    fine( f10, f2i  : f2f  : 2) -
c    &                    fine( f12, f2i-1: f2f-1: 2) -
c    &                    fine( f12, f2i  : f2f  : 2)  )
c    &        - F15o256*( fine( f11, f2i-2: f2f-2: 2) +
c    &                    fine( f11, f2i+1: f2f+1: 2)  )
c    &         + F9o256*( fine( f13, f2i-1: f2f-1: 2) +
c    &                    fine( f13, f2i  : f2f  : 2)  )
c    &         - F5o256*( fine( f10, f2i-2: f2f-2: 2) +
c    &                    fine( f10, f2i+1: f2f+1: 2) -
c    &                    fine( f12, f2i-2: f2f-2: 2) -
c    &                    fine( f12, f2i+1: f2f+1: 2)  )
c    &         - F1o256*( fine( f13, f2i-2: f2f-2: 2) +
c    &                    fine( f13, f2i+1: f2f+1: 2)  )
c-----------------------------------------------------------------------
c F77 version
        f21 = f2i
        do c2 = c2i,c2f
          coarse(c1,c2) =
     &         F135o256*( fine( f11, f21-1) +
     &                    fine( f11, f21  )  )
     &        + F45o256*( fine( f10, f21-1) +
     &                    fine( f10, f21  ) -
     &                    fine( f12, f21-1) -
     &                    fine( f12, f21  )  )
     &        - F15o256*( fine( f11, f21-2) +
     &                    fine( f11, f21+1)  )
     &         + F9o256*( fine( f13, f21-1) +
     &                    fine( f13, f21  )  )
     &         - F5o256*( fine( f10, f21-2) +
     &                    fine( f10, f21+1) -
     &                    fine( f12, f21-2) -
     &                    fine( f12, f21+1)  )
     &         - F1o256*( fine( f13, f21-2) +
     &                    fine( f13, f21+1)  )
          f21 = f21 + 2
        end do
c-----------------------------------------------------------------------
      end if
c
c Set cint(4) boundary
      if(c2f.lt.cint(4)) then
        c2 = cint(4)
        f20 = 2*c2
        f21 = f20 - 1
        f22 = f21 - 1
        f23 = f22 - 1
c-----------------------------------------------------------------------
c F90 version
c       coarse(c1i:c1f,c2) =
c    &         F135o256*( fine( f1i-1: f1f-1: 2, f21) +
c    &                    fine( f1i  : f1f  : 2, f21)  )
c    &        + F45o256*( fine( f1i-1: f1f-1: 2, f20) +
c    &                    fine( f1i  : f1f  : 2, f20) -
c    &                    fine( f1i-1: f1f-1: 2, f22) -
c    &                    fine( f1i  : f1f  : 2, f22)  )
c    &        - F15o256*( fine( f1i-2: f1f-2: 2, f21) +
c    &                    fine( f1i+1: f1f+1: 2, f21)  )
c    &         + F9o256*( fine( f1i-1: f1f-1: 2, f23) +
c    &                    fine( f1i  : f1f  : 2, f23)  )
c    &         - F5o256*( fine( f1i-2: f1f-2: 2, f20) +
c    &                    fine( f1i+1: f1f+1: 2, f20) -
c    &                    fine( f1i-2: f1f-2: 2, f22) -
c    &                    fine( f1i+1: f1f+1: 2, f22)  )
c    &         - F1o256*( fine( f1i-2: f1f-2: 2, f23) +
c    &                    fine( f1i+1: f1f+1: 2, f23)  )
c-----------------------------------------------------------------------
c F77 version
        f11 = f1i
        do c1 = c1i,c1f
          coarse(c1,c2) =
     &         F135o256*( fine( f11-1, f21) +
     &                    fine( f11  , f21)  )
     &        + F45o256*( fine( f11-1, f20) +
     &                    fine( f11  , f20) -
     &                    fine( f11-1, f22) -
     &                    fine( f11  , f22)  )
     &        - F15o256*( fine( f11-2, f21) +
     &                    fine( f11+1, f21)  )
     &         + F9o256*( fine( f11-1, f23) +
     &                    fine( f11  , f23)  )
     &         - F5o256*( fine( f11-2, f20) +
     &                    fine( f11+1, f20) -
     &                    fine( f11-2, f22) -
     &                    fine( f11+1, f22)  )
     &         - F1o256*( fine( f11-2, f23) +
     &                    fine( f11+1, f23)  )
          f11 = f11 + 2
        end do
c-----------------------------------------------------------------------
c
c Set cint(4)-cint(3) corner
        if(c1f.lt.cint(3)) then
          c1 = cint(3)
          f10 = 2*c1
          f11 = f10 - 1
          f12 = f11 - 1
          f13 = f12 - 1
          coarse(c1,c2) =
     &         F225o256*( fine( f11, f21)  )
     &       +  F75o256*( fine( f11, f20) +
     &                    fine( f10, f21) -
     &                    fine( f12, f21) -
     &                    fine( f11, f22)  )
     &       +  F25o256*( fine( f10, f20) -
     &                    fine( f12, f20) -
     &                    fine( f10, f22) +
     &                    fine( f12, f22)  )
     &       +  F15o256*( fine( f13, f21) +
     &                    fine( f11, f23)  )
     &       +   F5o256*( fine( f13, f20) -
     &                    fine( f13, f22) +
     &                    fine( f10, f23) -
     &                    fine( f12, f23)  )
     &       +   F1o256*( fine( f13, f23)  )
        end if
      end if
c
      return
      end
c*******************************************************************************
c* Fine to coarse grid residual transfer operation
c*     fine: The fine grid.
c*   coarse: The coarse grid.
c*  fext(4): The "extent" of the fine grid,
c*           fine(fext(1):fext(3),fext(2):fext(4))
c*  cint(4): The "interior" of the coarse grid.  The range of indices which
c*           are to be filled by the interpolation.
c*  cext(4): The "extent" of the coarse grid,
c*           coarse(cext(1):cext(3),cext(2):cext(4))
c*******************************************************************************
      subroutine resid_trans(fine,coarse,fext,cint,cext)
      implicit none
      integer fext(4),cint(4),cext(4)
      real*8 fine(fext(1):fext(3),fext(2):fext(4))
      real*8 coarse(cext(1):cext(3),cext(2):cext(4))
      integer c1i,c1f,c2i,c2f
      integer f1i,f1f,f2i,f2f
      integer c1,c2,f10,f11,f12,f13,f20,f21,f22,f23
c
c Set the coarse interior to be filled
      c1i = cint(1)
      f1i = 2*c1i
      c2i = cint(2)
      f2i = 2*c2i
      c1f = cint(3)
      f1f = 2*c1f
      c2f = cint(4)
      f2f = 2*c2f
c
c-----------------------------------------------------------------------
c F90 version
c     coarse(c1i:c1f,c2i:c2f) =
c    &                 fine( f1i-1: f1f-1: 2, f2i-1: f2f-1: 2) +
c    &                 fine( f1i  : f1f  : 2, f2i-1: f2f-1: 2) +
c    &                 fine( f1i-1: f1f-1: 2, f2i  : f2f  : 2) +
c    &                 fine( f1i  : f1f  : 2, f2i  : f2f  : 2)
c-----------------------------------------------------------------------
c F77 version
      f21 = f2i
      do c2 = c2i,c2f
        f11 = f1i
        do c1 = c1i,c1f
          coarse(c1,c2) =
     &                 fine( f11-1, f21-1) +
     &                 fine( f11  , f21-1) +
     &                 fine( f11-1, f21  ) +
     &                 fine( f11  , f21  )
          f11 = f11 + 2
        end do
        f21 = f21 + 2
      end do
c
      return
      end
