c*******************************************************************************
c* Coarse to fine grid interpolation using 3X3 bi-quadratic 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*  cint(4): The "interior" of the coarse grid.  Used to determine how much
c*           of the fine grid to fill (see implimentation notes below).
c*  cext(4): The "extent" of the coarse grid,
c*           coarse(cext(1):cext(3),cext(2):cext(4))
c*  cfil(4): The "fill" of the coarse grid.  The range of indices for which
c*           the coarse grid containts valid points for use in interpolation.
c*******************************************************************************
c* Implimentation notes:
c*    cint must always be a superset of fint, ie.
c*          fint(1) >= 2*cint(1) - 1,
c*          fint(2) >= 2*cint(2) - 1,
c*          fint(3) <= 2*cint(3)    ,
c*          fint(4) <= 2*cint(4)    .
c*    Usually, cint "=" fint.  When adaptive gridding is used, fint may be
c*    defined one fine zone interior to cint.  This buffer zone is filled
c*    because cint is used to define the fine grid interpolation fill region
c*    instead of fint.  These filled buffer zones constitute the Dirichlet
c*    boundary points for the fine grid.
c*******************************************************************************
      subroutine ctof3x3(fine,coarse,fext,cint,cext,cfil)
      implicit none
      integer fext(4),cint(4),cext(4),cfil(4)
      real*8 fine(fext(1):fext(3),fext(2):fext(4))
      real*8 coarse(cext(1):cext(3),cext(2):cext(4))
      real*8 tmp0,tmp1,tmp2,tmp3,tmp4
      integer c1i,c1f,c2i,c2f
      integer f1i,f1f,f2i,f2f
      integer c1,c2,f10,f11,f12,f13,f20,f21,f22,f23
      real*8 F2025o1024, F675o512,  F945o1024, F225o256,  F405o512,
     &       F315o512,   F135o256,  F441o1024, F105o256,  F189o512,
     &       F81o256,    F147o512,  F63o256,   F225o1024, F49o256,
     &       F75o512,    F135o1024, F105o1024, F45o512,   F35o512,
     &       F63o1024,   F27o512,   F21o512,   F25o1024,  F15o1024,
     &       F9o1024
      parameter(F2025o1024= 2025.D0/1024.D0,  F675o512= 675.D0 /512.D0,
     &           F945o1024=  945.D0/1024.D0,  F225o256= 225.D0 /256.D0,
     &            F405o512=  405.D0 /512.D0,  F315o512= 315.D0 /512.D0,
     &            F135o256=  135.D0 /256.D0, F441o1024= 441.D0/1024.D0,
     &            F105o256=  105.D0 /256.D0,  F189o512= 189.D0 /512.D0,
     &             F81o256=   81.D0 /256.D0,  F147o512= 147.D0 /512.D0,
     &             F63o256=   63.D0 /256.D0, F225o1024= 225.D0/1024.D0,
     &             F49o256=   49.D0 /256.D0,   F75o512=  75.D0 /512.D0,
     &           F135o1024=  135.D0/1024.D0, F105o1024= 105.D0/1024.D0,
     &             F45o512=   45.D0 /512.D0,   F35o512=  35.D0 /512.D0,
     &            F63o1024=   63.D0/1024.D0,   F27o512=  27.D0 /512.D0,
     &             F21o512=   21.D0 /512.D0,  F25o1024=  25.D0/1024.D0,
     &            F15o1024=   15.D0/1024.D0,   F9o1024=   9.D0/1024.D0 )
c
c Make sure the "fill" region of the coarse grid is at least 3X3.
      if(((cfil(3) - cfil(1)).lt.2).or.((cfil(4) - cfil(2)).lt.2)) then
        write(*,*) "ctof3x3: fill region smaller than 3X3!"
        stop
      end if
c
c Set the coarse interior which will be used to fill the fine grid using
c "centered" 3X3 interpolation.
      c1i = MAX(cint(1),cfil(1) + 1)
      f1i = 2*c1i
      c2i = MAX(cint(2),cfil(2) + 1)
      f2i = 2*c2i
      c1f = MIN(cint(3),cfil(3) - 1)
      f1f = 2*c1f
      c2f = MIN(cint(4),cfil(4) - 1)
      f2f = 2*c2f
c-----------------------------------------------------------------------
c F90 version
c     fine(f1i-1: f1f-1: 2, f2i-1: f2f-1: 2) = 
c    &         F225o256*( coarse(c1i  :c1f  ,c2i  :c2f  )  )
c    &       +  F75o512*( coarse(c1i-1:c1f-1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i-1:c2f-1)  )
c    &       -  F45o512*( coarse(c1i+1:c1f+1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i+1:c2f+1)  )
c    &       + F25o1024*( coarse(c1i-1:c1f-1,c2i-1:c2f-1)  )
c    &       - F15o1024*( coarse(c1i+1:c1f+1,c2i-1:c2f-1) +
c    &                    coarse(c1i-1:c1f-1,c2i+1:c2f+1)  )
c    &       +  F9o1024*( coarse(c1i+1:c1f+1,c2i+1:c2f+1)  )
c     fine(f1i  : f1f  : 2, f2i-1: f2f-1: 2) = 
c    &         F225o256*( coarse(c1i  :c1f  ,c2i  :c2f  )  )
c    &       +  F75o512*( coarse(c1i+1:c1f+1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i-1:c2f-1)  )
c    &       -  F45o512*( coarse(c1i-1:c1f-1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i+1:c2f+1)  )
c    &       + F25o1024*( coarse(c1i+1:c1f+1,c2i-1:c2f-1)  )
c    &       - F15o1024*( coarse(c1i-1:c1f-1,c2i-1:c2f-1) +
c    &                    coarse(c1i+1:c1f+1,c2i+1:c2f+1)  )
c    &       +  F9o1024*( coarse(c1i-1:c1f-1,c2i+1:c2f+1)  )
c     fine(f1i-1: f1f-1: 2, f2i  : f2f  : 2) = 
c    &         F225o256*( coarse(c1i  :c1f  ,c2i  :c2f  )  )
c    &       +  F75o512*( coarse(c1i-1:c1f-1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i+1:c2f+1)  )
c    &       -  F45o512*( coarse(c1i+1:c1f+1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i-1:c2f-1)  )
c    &       + F25o1024*( coarse(c1i-1:c1f-1,c2i+1:c2f+1)  )
c    &       - F15o1024*( coarse(c1i-1:c1f-1,c2i-1:c2f-1) +
c    &                    coarse(c1i+1:c1f+1,c2i+1:c2f+1)  )
c    &       +  F9o1024*( coarse(c1i+1:c1f+1,c2i-1:c2f-1)  )
c     fine(f1i  : f1f  : 2, f2i  : f2f  : 2) = 
c    &         F225o256*( coarse(c1i  :c1f  ,c2i  :c2f  )  )
c    &       +  F75o512*( coarse(c1i+1:c1f+1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i+1:c2f+1)  )
c    &       -  F45o512*( coarse(c1i-1:c1f-1,c2i  :c2f  ) +
c    &                    coarse(c1i  :c1f  ,c2i-1:c2f-1)  )
c    &       + F25o1024*( coarse(c1i+1:c1f+1,c2i+1:c2f+1)  )
c    &       - F15o1024*( coarse(c1i+1:c1f+1,c2i-1:c2f-1) +
c    &                    coarse(c1i-1:c1f-1,c2i+1:c2f+1)  )
c    &       +  F9o1024*( coarse(c1i-1:c1f-1,c2i-1:c2f-1)  )
c-----------------------------------------------------------------------
c F77 version
      f21 = f2i
      do c2 = c2i,c2f
        f11 = f1i
        do c1 = c1i,c1f
          tmp0 =
     &         F225o256*( coarse(c1  ,c2  )  )
          tmp1 = tmp0
     &       - F15o1024*( coarse(c1+1,c2-1) +
     &                    coarse(c1-1,c2+1)  )
          tmp2 = tmp0
     &       - F15o1024*( coarse(c1-1,c2-1) +
     &                    coarse(c1+1,c2+1)  )
          fine(f11-1, f21-1) = tmp1
     &       +  F75o512*( coarse(c1-1,c2  ) +
     &                    coarse(c1  ,c2-1)  )
     &       -  F45o512*( coarse(c1+1,c2  ) +
     &                    coarse(c1  ,c2+1)  )
     &       + F25o1024*( coarse(c1-1,c2-1)  )
     &       +  F9o1024*( coarse(c1+1,c2+1)  )
          fine(f11  , f21  ) = tmp1
     &       +  F75o512*( coarse(c1+1,c2  ) +
     &                    coarse(c1  ,c2+1)  )
     &       -  F45o512*( coarse(c1-1,c2  ) +
     &                    coarse(c1  ,c2-1)  )
     &       + F25o1024*( coarse(c1+1,c2+1)  )
     &       +  F9o1024*( coarse(c1-1,c2-1)  )
          fine(f11  , f21-1) = tmp2
     &       +  F75o512*( coarse(c1+1,c2  ) +
     &                    coarse(c1  ,c2-1)  )
     &       -  F45o512*( coarse(c1-1,c2  ) +
     &                    coarse(c1  ,c2+1)  )
     &       + F25o1024*( coarse(c1+1,c2-1)  )
     &       +  F9o1024*( coarse(c1-1,c2+1)  )
          fine(f11-1, f21  ) = tmp2
     &       +  F75o512*( coarse(c1-1,c2  ) +
     &                    coarse(c1  ,c2+1)  )
     &       -  F45o512*( coarse(c1+1,c2  ) +
     &                    coarse(c1  ,c2-1)  )
     &       + F25o1024*( coarse(c1-1,c2+1)  )
     &       +  F9o1024*( coarse(c1+1,c2-1)  )
c     write(*,1000) "C ",c1,c2,fine(f11-1,f21-1),fine(f11  ,f21  ),
c    .                         fine(f11  ,f21-1),fine(f11-1,f21  )
          f11 = f11 + 2
        end do
        f21 = f21 + 2
      end do
1000  format(a2,1p,2i4,4e20.13)
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     fine(f10, f2i-1: f2f-1: 2) = 
c    &         F675o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +   F75o512*( coarse(c1+2, c2i  :c2f  )  )
c    &      -  F135o256*( coarse(c1+1, c2i  :c2f  )  )
c    &      + F225o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      - F135o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      -   F45o512*( coarse(c1+1, c2i-1:c2f-1)  )
c    &      +   F27o512*( coarse(c1+1, c2i+1:c2f+1)  )
c    &      +  F25o1024*( coarse(c1+2, c2i-1:c2f-1)  )
c    &      -  F15o1024*( coarse(c1+2, c2i+1:c2f+1)  )
c     fine(f11, f2i-1: f2f-1: 2) = 
c    &         F315o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +  F105o256*( coarse(c1+1, c2i  :c2f  )  )
c    &      -   F45o512*( coarse(c1+2, c2i  :c2f  )  )
c    &      + F105o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      -   F21o512*( coarse(c1+1, c2i+1:c2f+1)  )
c    &      -  F63o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      +   F35o512*( coarse(c1+1, c2i-1:c2f-1)  )
c    &      +   F9o1024*( coarse(c1+2, c2i+1:c2f+1)  )
c    &      -  F15o1024*( coarse(c1+2, c2i-1:c2f-1)  )
c     fine(f10, f2i  : f2f  : 2) = 
c    &         F675o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +   F75o512*( coarse(c1+2, c2i  :c2f  )  )
c    &      -  F135o256*( coarse(c1+1, c2i  :c2f  )  )
c    &      + F225o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      - F135o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      +   F27o512*( coarse(c1+1, c2i-1:c2f-1)  )
c    &      +  F25o1024*( coarse(c1+2, c2i+1:c2f+1)  )
c    &      -   F45o512*( coarse(c1+1, c2i+1:c2f+1)  )
c    &      -  F15o1024*( coarse(c1+2, c2i-1:c2f-1)  )
c     fine(f11, f2i  : f2f  : 2) = 
c    &         F315o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +  F105o256*( coarse(c1+1, c2i  :c2f  )  )
c    &      -   F45o512*( coarse(c1+2, c2i  :c2f  )  )
c    &      + F105o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      +   F35o512*( coarse(c1+1, c2i+1:c2f+1)  )
c    &      -  F63o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      -   F21o512*( coarse(c1+1, c2i-1:c2f-1)  )
c    &      +   F9o1024*( coarse(c1+2, c2i-1:c2f-1)  )
c    &      -  F15o1024*( coarse(c1+2, c2i+1:c2f+1)  )
c-----------------------------------------------------------------------
c F77 version
        f21 = f2i
        do c2 = c2i,c2f
          tmp1 = 
     &         F675o512*( coarse(c1  , c2  )  )
     &      +   F75o512*( coarse(c1+2, c2  )  )
     &      -  F135o256*( coarse(c1+1, c2  )  )
          tmp2 = 
     &         F315o512*( coarse(c1  , c2  )  )
     &      +  F105o256*( coarse(c1+1, c2  )  )
     &      -   F45o512*( coarse(c1+2, c2  )  )
          tmp3 = 
     &      -  F15o1024*( coarse(c1+2, c2+1)  )
          tmp4 = 
     &      -  F15o1024*( coarse(c1+2, c2-1)  )
          fine(f10, f21-1) = tmp1 + tmp3
     &      + F225o1024*( coarse(c1  , c2-1)  )
     &      - F135o1024*( coarse(c1  , c2+1)  )
     &      -   F45o512*( coarse(c1+1, c2-1)  )
     &      +   F27o512*( coarse(c1+1, c2+1)  )
     &      +  F25o1024*( coarse(c1+2, c2-1)  )
          fine(f11, f21-1) = tmp2 + tmp4
     &      + F105o1024*( coarse(c1  , c2-1)  )
     &      -   F21o512*( coarse(c1+1, c2+1)  )
     &      -  F63o1024*( coarse(c1  , c2+1)  )
     &      +   F35o512*( coarse(c1+1, c2-1)  )
     &      +   F9o1024*( coarse(c1+2, c2+1)  )
          fine(f10, f21  ) = tmp1 + tmp4
     &      + F225o1024*( coarse(c1  , c2+1)  )
     &      - F135o1024*( coarse(c1  , c2-1)  )
     &      +   F27o512*( coarse(c1+1, c2-1)  )
     &      +  F25o1024*( coarse(c1+2, c2+1)  )
     &      -   F45o512*( coarse(c1+1, c2+1)  )
          fine(f11, f21  ) = tmp2 + tmp3
     &      + F105o1024*( coarse(c1  , c2+1)  )
     &      +   F35o512*( coarse(c1+1, c2+1)  )
     &      -  F63o1024*( coarse(c1  , c2-1)  )
     &      -   F21o512*( coarse(c1+1, c2-1)  )
     &      +   F9o1024*( coarse(c1+2, c2-1)  )
c     write(*,1000) "B1",c1,c2,fine(f10  ,f21-1),fine(f11  ,f21-1),
c    .                         fine(f10  ,f21  ),fine(f11  ,f21  )
          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
          tmp1 =
     &        F945o1024*( coarse(c1  , c2  )  )
     &     -    F63o256*( coarse(c1+1, c2-1)  )
     &     -   F15o1024*( coarse(c1+2, c2-2)  )
          fine(f10, f21) = tmp1
     &     -   F189o512*( coarse(c1+1, c2  )  )
     &     +   F315o512*( coarse(c1  , c2-1)  )
     &     +  F105o1024*( coarse(c1+2, c2  )  )
     &     -  F135o1024*( coarse(c1  , c2-2)  )
     &     +    F35o512*( coarse(c1+2, c2-1)  )
     &     +    F27o512*( coarse(c1+1, c2-2)  )
          fine(f11, f21) =
     &        F441o1024*( coarse(c1  , c2  )  )
     &     +    F49o256*( coarse(c1+1, c2-1)  )
     &     +    F9o1024*( coarse(c1+2, c2-2)  )
     &     +   F147o512*( coarse(c1+1, c2  ) +
     &                    coarse(c1  , c2-1)  )
     &     -   F63o1024*( coarse(c1  , c2-2) +
     &                    coarse(c1+2, c2  )  )
     &     -    F21o512*( coarse(c1+1, c2-2) +
     &                    coarse(c1+2, c2-1)  )
          fine(f10, f20) =
     &       F2025o1024*( coarse(c1  , c2  )  )
     &     +    F81o256*( coarse(c1+1, c2-1)  )
     &     +   F25o1024*( coarse(c1+2, c2-2)  )
     &     -   F405o512*( coarse(c1+1, c2  ) +
     &                    coarse(c1  , c2-1)  )
     &     +  F225o1024*( coarse(c1  , c2-2) +
     &                    coarse(c1+2, c2  )  )
     &     -    F45o512*( coarse(c1+1, c2-2) +
     &                    coarse(c1+2, c2-1)  )
          fine(f11, f20) = tmp1
     &     -   F189o512*( coarse(c1  , c2-1)  )
     &     +   F315o512*( coarse(c1+1, c2  )  )
     &     +  F105o1024*( coarse(c1  , c2-2)  )
     &     -  F135o1024*( coarse(c1+2, c2  )  )
     &     +    F35o512*( coarse(c1+1, c2-2)  )
     &     +    F27o512*( coarse(c1+2, c2-1)  )
c     write(*,1000) "14",c1,c2,fine(f10  ,f21  ),fine(f11  ,f21  ),
c    .                         fine(f10  ,f20  ),fine(f11  ,f20  )
        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
          tmp1 =
     &        F945o1024*( coarse(c1  , c2  )  )
     &     -    F63o256*( coarse(c1+1, c2+1)  )
     &     -   F15o1024*( coarse(c1+2, c2+2)  )
          fine(f10, f20) =
     &       F2025o1024*( coarse(c1  , c2  )  )
     &     +    F81o256*( coarse(c1+1, c2+1)  )
     &     +   F25o1024*( coarse(c1+2, c2+2)  )
     &     -   F405o512*( coarse(c1+1, c2  ) +
     &                    coarse(c1  , c2+1)  )
     &     +  F225o1024*( coarse(c1  , c2+2) +
     &                    coarse(c1+2, c2  )  )
     &     -    F45o512*( coarse(c1+1, c2+2) +
     &                    coarse(c1+2, c2+1)  )
          fine(f11, f20) = tmp1
     &     -   F189o512*( coarse(c1  , c2+1)  )
     &     +   F315o512*( coarse(c1+1, c2  )  )
     &     +  F105o1024*( coarse(c1  , c2+2)  )
     &     -  F135o1024*( coarse(c1+2, c2  )  )
     &     +    F35o512*( coarse(c1+1, c2+2)  )
     &     +    F27o512*( coarse(c1+2, c2+1)  )
          fine(f10, f21) = tmp1
     &     -   F189o512*( coarse(c1+1, c2  )  )
     &     +   F315o512*( coarse(c1  , c2+1)  )
     &     +  F105o1024*( coarse(c1+2, c2  )  )
     &     -  F135o1024*( coarse(c1  , c2+2)  )
     &     +    F35o512*( coarse(c1+2, c2+1)  )
     &     +    F27o512*( coarse(c1+1, c2+2)  )
          fine(f11, f21) =
     &        F441o1024*( coarse(c1  , c2  )  )
     &     +    F49o256*( coarse(c1+1, c2+1)  )
     &     +    F9o1024*( coarse(c1+2, c2+2)  )
     &     +   F147o512*( coarse(c1+1, c2  ) +
     &                    coarse(c1  , c2+1)  )
     &     -   F63o1024*( coarse(c1  , c2+2) +
     &                    coarse(c1+2, c2  )  )
     &     -    F21o512*( coarse(c1+1, c2+2) +
     &                    coarse(c1+2, c2+1)  )
c     write(*,1000) "21",c1,c2,fine(f10  ,f20  ),fine(f11  ,f20  ),
c    .                         fine(f10  ,f21  ),fine(f11  ,f21  )
        end if
c-----------------------------------------------------------------------
c F90 version
c     fine(f1i-1: f1f-1: 2, f20) = 
c    &         F675o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +   F75o512*( coarse(c1i  :c1f  , c2+2)  )
c    &      -  F135o256*( coarse(c1i  :c1f  , c2+1)  )
c    &      + F225o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      - F135o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      +   F27o512*( coarse(c1i+1:c1f+1, c2+1)  )
c    &      +  F25o1024*( coarse(c1i-1:c1f-1, c2+2)  )
c    &      -   F45o512*( coarse(c1i-1:c1f-1, c2+1)  )
c    &      -  F15o1024*( coarse(c1i+1:c1f+1, c2+2)  )
c     fine(f1i  : f1f  : 2, f20) = 
c    &         F675o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +   F75o512*( coarse(c1i  :c1f  , c2+2)  )
c    &      -  F135o256*( coarse(c1i  :c1f  , c2+1)  )
c    &      + F225o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      - F135o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      -   F45o512*( coarse(c1i+1:c1f+1, c2+1)  )
c    &      +   F27o512*( coarse(c1i-1:c1f-1, c2+1)  )
c    &      +  F25o1024*( coarse(c1i+1:c1f+1, c2+2)  )
c    &      -  F15o1024*( coarse(c1i-1:c1f-1, c2+2)  )
c     fine(f1i-1: f1f-1: 2, f21) = 
c    &         F315o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +  F105o256*( coarse(c1i  :c1f  , c2+1)  )
c    &      -   F45o512*( coarse(c1i  :c1f  , c2+2)  )
c    &      + F105o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      +   F35o512*( coarse(c1i-1:c1f-1, c2+1)  )
c    &      -  F63o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      -   F21o512*( coarse(c1i+1:c1f+1, c2+1)  )
c    &      +   F9o1024*( coarse(c1i+1:c1f+1, c2+2)  )
c    &      -  F15o1024*( coarse(c1i-1:c1f-1, c2+2)  )
c     fine(f1i  : f1f  : 2, f21) = 
c    &         F315o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +  F105o256*( coarse(c1i  :c1f  , c2+1)  )
c    &      -   F45o512*( coarse(c1i  :c1f  , c2+2)  )
c    &      + F105o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      -   F21o512*( coarse(c1i-1:c1f-1, c2+1)  )
c    &      -  F63o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      +   F35o512*( coarse(c1i+1:c1f+1, c2+1)  )
c    &      +   F9o1024*( coarse(c1i-1:c1f-1, c2+2)  )
c    &      -  F15o1024*( coarse(c1i+1:c1f+1, c2+2)  )
c-----------------------------------------------------------------------
c F77 version
        f11 = f1i
        do c1 = c1i,c1f
          tmp1 = 
     &         F675o512*( coarse(c1  , c2  )  )
     &      +   F75o512*( coarse(c1  , c2+2)  )
     &      -  F135o256*( coarse(c1  , c2+1)  )
          tmp2 = 
     &         F315o512*( coarse(c1  , c2  )  )
     &      +  F105o256*( coarse(c1  , c2+1)  )
     &      -   F45o512*( coarse(c1  , c2+2)  )
          tmp3 = 
     &      -  F15o1024*( coarse(c1-1, c2+2)  )
          tmp4 = 
     &      -  F15o1024*( coarse(c1+1, c2+2)  )
          fine(f11-1, f20) = tmp1 + tmp4
     &      + F225o1024*( coarse(c1-1, c2  )  )
     &      - F135o1024*( coarse(c1+1, c2  )  )
     &      +   F27o512*( coarse(c1+1, c2+1)  )
     &      +  F25o1024*( coarse(c1-1, c2+2)  )
     &      -   F45o512*( coarse(c1-1, c2+1)  )
          fine(f11  , f20) = tmp1 + tmp3
     &      + F225o1024*( coarse(c1+1, c2  )  )
     &      - F135o1024*( coarse(c1-1, c2  )  )
     &      -   F45o512*( coarse(c1+1, c2+1)  )
     &      +   F27o512*( coarse(c1-1, c2+1)  )
     &      +  F25o1024*( coarse(c1+1, c2+2)  )
          fine(f11-1, f21) = tmp2 + tmp3
     &      + F105o1024*( coarse(c1-1, c2  )  )
     &      +   F35o512*( coarse(c1-1, c2+1)  )
     &      -  F63o1024*( coarse(c1+1, c2  )  )
     &      -   F21o512*( coarse(c1+1, c2+1)  )
     &      +   F9o1024*( coarse(c1+1, c2+2)  )
          fine(f11  , f21) = tmp2 + tmp4
     &      + F105o1024*( coarse(c1+1, c2  )  )
     &      -   F21o512*( coarse(c1-1, c2+1)  )
     &      -  F63o1024*( coarse(c1-1, c2  )  )
     &      +   F35o512*( coarse(c1+1, c2+1)  )
     &      +   F9o1024*( coarse(c1-1, c2+2)  )
c     write(*,1000) "B2",c1,c2,fine(f11-1,f20  ),fine(f11  ,f20  ),
c    .                         fine(f11-1,f21  ),fine(f11  ,f21  )
          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
          tmp1 =
     &        F945o1024*( coarse(c1  , c2  )  )
     &     -    F63o256*( coarse(c1-1, c2+1)  )
     &     -   F15o1024*( coarse(c1-2, c2+2)  )
          fine(f11, f20) = tmp1
     &     -   F189o512*( coarse(c1  , c2+1)  )
     &     +   F315o512*( coarse(c1-1, c2  )  )
     &     +  F105o1024*( coarse(c1  , c2+2)  )
     &     -  F135o1024*( coarse(c1-2, c2  )  )
     &     +    F35o512*( coarse(c1-1, c2+2)  )
     &     +    F27o512*( coarse(c1-2, c2+1)  )
          fine(f10, f20) =
     &       F2025o1024*( coarse(c1  , c2  )  )
     &     +    F81o256*( coarse(c1-1, c2+1)  )
     &     +   F25o1024*( coarse(c1-2, c2+2)  )
     &     -   F405o512*( coarse(c1-1, c2  ) +
     &                    coarse(c1  , c2+1)  )
     &     +  F225o1024*( coarse(c1  , c2+2) +
     &                    coarse(c1-2, c2  )  )
     &     -    F45o512*( coarse(c1-1, c2+2) +
     &                    coarse(c1-2, c2+1)  )
          fine(f11, f21) =
     &        F441o1024*( coarse(c1  , c2  )  )
     &     +    F49o256*( coarse(c1-1, c2+1)  )
     &     +    F9o1024*( coarse(c1-2, c2+2)  )
     &     +   F147o512*( coarse(c1-1, c2  ) +
     &                    coarse(c1  , c2+1)  )
     &     -   F63o1024*( coarse(c1  , c2+2) +
     &                    coarse(c1-2, c2  )  )
     &     -    F21o512*( coarse(c1-1, c2+2) +
     &                    coarse(c1-2, c2+1)  )
          fine(f10, f21) = tmp1
     &     -   F189o512*( coarse(c1-1, c2  )  )
     &     +   F315o512*( coarse(c1  , c2+1)  )
     &     +  F105o1024*( coarse(c1-2, c2  )  )
     &     -  F135o1024*( coarse(c1  , c2+2)  )
     &     +    F35o512*( coarse(c1-2, c2+1)  )
     &     +    F27o512*( coarse(c1-1, c2+2)  )
c     write(*,1000) "32",c1,c2,fine(f10  ,f20  ),fine(f11  ,f21  ),
c    .                         fine(f11  ,f21  ),fine(f10  ,f21  )
        end if
c-----------------------------------------------------------------------
c F90 version
c     fine(f11, f2i-1: f2f-1: 2) = 
c    &         F315o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +  F105o256*( coarse(c1-1, c2i  :c2f  )  )
c    &      -   F45o512*( coarse(c1-2, c2i  :c2f  )  )
c    &      + F105o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      -   F21o512*( coarse(c1-1, c2i+1:c2f+1)  )
c    &      -  F63o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      +   F35o512*( coarse(c1-1, c2i-1:c2f-1)  )
c    &      +   F9o1024*( coarse(c1-2, c2i+1:c2f+1)  )
c    &      -  F15o1024*( coarse(c1-2, c2i-1:c2f-1)  )
c     fine(f10, f2i-1: f2f-1: 2) = 
c    &         F675o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +   F75o512*( coarse(c1-2, c2i  :c2f  )  )
c    &      -  F135o256*( coarse(c1-1, c2i  :c2f  )  )
c    &      + F225o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      - F135o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      -   F45o512*( coarse(c1-1, c2i-1:c2f-1)  )
c    &      +   F27o512*( coarse(c1-1, c2i+1:c2f+1)  )
c    &      +  F25o1024*( coarse(c1-2, c2i-1:c2f-1)  )
c    &      -  F15o1024*( coarse(c1-2, c2i+1:c2f+1)  )
c     fine(f11, f2i  : f2f  : 2) = 
c    &         F315o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +  F105o256*( coarse(c1-1, c2i  :c2f  )  )
c    &      -   F45o512*( coarse(c1-2, c2i  :c2f  )  )
c    &      + F105o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      +   F35o512*( coarse(c1-1, c2i+1:c2f+1)  )
c    &      -  F63o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      -   F21o512*( coarse(c1-1, c2i-1:c2f-1)  )
c    &      +   F9o1024*( coarse(c1-2, c2i-1:c2f-1)  )
c    &      -  F15o1024*( coarse(c1-2, c2i+1:c2f+1)  )
c     fine(f10, f2i  : f2f  : 2) = 
c    &         F675o512*( coarse(c1  , c2i  :c2f  )  )
c    &      +   F75o512*( coarse(c1-2, c2i  :c2f  )  )
c    &      -  F135o256*( coarse(c1-1, c2i  :c2f  )  )
c    &      + F225o1024*( coarse(c1  , c2i+1:c2f+1)  )
c    &      - F135o1024*( coarse(c1  , c2i-1:c2f-1)  )
c    &      +   F27o512*( coarse(c1-1, c2i-1:c2f-1)  )
c    &      +  F25o1024*( coarse(c1-2, c2i+1:c2f+1)  )
c    &      -   F45o512*( coarse(c1-1, c2i+1:c2f+1)  )
c    &      -  F15o1024*( coarse(c1-2, c2i-1:c2f-1)  )
c-----------------------------------------------------------------------
c F77 version
        f21 = f2i
        do c2 = c2i,c2f
          tmp1 = 
     &         F675o512*( coarse(c1  , c2  )  )
     &      +   F75o512*( coarse(c1-2, c2  )  )
     &      -  F135o256*( coarse(c1-1, c2  )  )
          tmp2 = 
     &         F315o512*( coarse(c1  , c2  )  )
     &      +  F105o256*( coarse(c1-1, c2  )  )
     &      -   F45o512*( coarse(c1-2, c2  )  )
          tmp3 = 
     &      -  F15o1024*( coarse(c1-2, c2+1)  )
          tmp4 = 
     &      -  F15o1024*( coarse(c1-2, c2-1)  )
          fine(f11, f21-1) = tmp2 + tmp4
     &      + F105o1024*( coarse(c1  , c2-1)  )
     &      -   F21o512*( coarse(c1-1, c2+1)  )
     &      -  F63o1024*( coarse(c1  , c2+1)  )
     &      +   F35o512*( coarse(c1-1, c2-1)  )
     &      +   F9o1024*( coarse(c1-2, c2+1)  )
          fine(f10, f21-1) = tmp1 + tmp3
     &      + F225o1024*( coarse(c1  , c2-1)  )
     &      - F135o1024*( coarse(c1  , c2+1)  )
     &      -   F45o512*( coarse(c1-1, c2-1)  )
     &      +   F27o512*( coarse(c1-1, c2+1)  )
     &      +  F25o1024*( coarse(c1-2, c2-1)  )
          fine(f11, f21  ) = tmp2 + tmp3
     &      + F105o1024*( coarse(c1  , c2+1)  )
     &      +   F35o512*( coarse(c1-1, c2+1)  )
     &      -  F63o1024*( coarse(c1  , c2-1)  )
     &      -   F21o512*( coarse(c1-1, c2-1)  )
     &      +   F9o1024*( coarse(c1-2, c2-1)  )
          fine(f10, f21  ) = tmp1 + tmp4
     &      + F225o1024*( coarse(c1  , c2+1)  )
     &      - F135o1024*( coarse(c1  , c2-1)  )
     &      +   F27o512*( coarse(c1-1, c2-1)  )
     &      +  F25o1024*( coarse(c1-2, c2+1)  )
     &      -   F45o512*( coarse(c1-1, c2+1)  )
c     write(*,1000) "B3",c1,c2,fine(f11  ,f21-1),fine(f10  ,f21-1),
c    .                         fine(f11  ,f21  ),fine(f10  ,f21  )
          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     fine(f1i-1: f1f-1: 2, f21) = 
c    &         F315o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +  F105o256*( coarse(c1i  :c1f  , c2-1)  )
c    &      -   F45o512*( coarse(c1i  :c1f  , c2-2)  )
c    &      + F105o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      +   F35o512*( coarse(c1i-1:c1f-1, c2-1)  )
c    &      -  F63o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      -   F21o512*( coarse(c1i+1:c1f+1, c2-1)  )
c    &      +   F9o1024*( coarse(c1i+1:c1f+1, c2-2)  )
c    &      -  F15o1024*( coarse(c1i-1:c1f-1, c2-2)  )
c     fine(f1i  : f1f  : 2, f21) = 
c    &         F315o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +  F105o256*( coarse(c1i  :c1f  , c2-1)  )
c    &      -   F45o512*( coarse(c1i  :c1f  , c2-2)  )
c    &      + F105o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      -   F21o512*( coarse(c1i-1:c1f-1, c2-1)  )
c    &      -  F63o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      +   F35o512*( coarse(c1i+1:c1f+1, c2-1)  )
c    &      +   F9o1024*( coarse(c1i-1:c1f-1, c2-2)  )
c    &      -  F15o1024*( coarse(c1i+1:c1f+1, c2-2)  )
c     fine(f1i-1: f1f-1: 2, f20) = 
c    &         F675o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +   F75o512*( coarse(c1i  :c1f  , c2-2)  )
c    &      -  F135o256*( coarse(c1i  :c1f  , c2-1)  )
c    &      + F225o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      - F135o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      +   F27o512*( coarse(c1i+1:c1f+1, c2-1)  )
c    &      +  F25o1024*( coarse(c1i-1:c1f-1, c2-2)  )
c    &      -   F45o512*( coarse(c1i-1:c1f-1, c2-1)  )
c    &      -  F15o1024*( coarse(c1i+1:c1f+1, c2-2)  )
c     fine(f1i  : f1f  : 2, f20) = 
c    &         F675o512*( coarse(c1i  :c1f  , c2  )  )
c    &      +   F75o512*( coarse(c1i  :c1f  , c2-2)  )
c    &      -  F135o256*( coarse(c1i  :c1f  , c2-1)  )
c    &      + F225o1024*( coarse(c1i+1:c1f+1, c2  )  )
c    &      - F135o1024*( coarse(c1i-1:c1f-1, c2  )  )
c    &      -   F45o512*( coarse(c1i+1:c1f+1, c2-1)  )
c    &      +   F27o512*( coarse(c1i-1:c1f-1, c2-1)  )
c    &      +  F25o1024*( coarse(c1i+1:c1f+1, c2-2)  )
c    &      -  F15o1024*( coarse(c1i-1:c1f-1, c2-2)  )
c-----------------------------------------------------------------------
c F77 version
        f11 = f1i
        do c1 = c1i,c1f
          tmp1 = 
     &         F675o512*( coarse(c1  , c2  )  )
     &      +   F75o512*( coarse(c1  , c2-2)  )
     &      -  F135o256*( coarse(c1  , c2-1)  )
          tmp2 = 
     &         F315o512*( coarse(c1  , c2  )  )
     &      +  F105o256*( coarse(c1  , c2-1)  )
     &      -   F45o512*( coarse(c1  , c2-2)  )
          tmp3 = 
     &      -  F15o1024*( coarse(c1-1, c2-2)  )
          tmp4 = 
     &      -  F15o1024*( coarse(c1+1, c2-2)  )
          fine(f11-1, f21) = tmp2 + tmp3
     &      + F105o1024*( coarse(c1-1, c2  )  )
     &      +   F35o512*( coarse(c1-1, c2-1)  )
     &      -  F63o1024*( coarse(c1+1, c2  )  )
     &      -   F21o512*( coarse(c1+1, c2-1)  )
     &      +   F9o1024*( coarse(c1+1, c2-2)  )
          fine(f11  , f21) = tmp2 + tmp4
     &      + F105o1024*( coarse(c1+1, c2  )  )
     &      -   F21o512*( coarse(c1-1, c2-1)  )
     &      -  F63o1024*( coarse(c1-1, c2  )  )
     &      +   F35o512*( coarse(c1+1, c2-1)  )
     &      +   F9o1024*( coarse(c1-1, c2-2)  )
          fine(f11-1, f20) = tmp1 + tmp4
     &      + F225o1024*( coarse(c1-1, c2  )  )
     &      - F135o1024*( coarse(c1+1, c2  )  )
     &      +   F27o512*( coarse(c1+1, c2-1)  )
     &      +  F25o1024*( coarse(c1-1, c2-2)  )
     &      -   F45o512*( coarse(c1-1, c2-1)  )
          fine(f11  , f20) = tmp1 + tmp3
     &      + F225o1024*( coarse(c1+1, c2  )  )
     &      - F135o1024*( coarse(c1-1, c2  )  )
     &      -   F45o512*( coarse(c1+1, c2-1)  )
     &      +   F27o512*( coarse(c1-1, c2-1)  )
     &      +  F25o1024*( coarse(c1+1, c2-2)  )
c     write(*,1000) "B4",c1,c2,fine(f11-1,f21  ),fine(f11  ,f21  ),
c    .                         fine(f11-1,f20  ),fine(f11  ,f20  )
          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
          tmp1 =
     &        F945o1024*( coarse(c1  , c2  )  )
     &     -    F63o256*( coarse(c1-1, c2-1)  )
     &     -   F15o1024*( coarse(c1-2, c2-2)  )
          fine(f11, f21) =
     &        F441o1024*( coarse(c1  , c2  )  )
     &     +    F49o256*( coarse(c1-1, c2-1)  )
     &     +    F9o1024*( coarse(c1-2, c2-2)  )
     &     +   F147o512*( coarse(c1-1, c2  ) +
     &                    coarse(c1  , c2-1)  )
     &     -   F63o1024*( coarse(c1  , c2-2) +
     &                    coarse(c1-2, c2  )  )
     &     -    F21o512*( coarse(c1-1, c2-2) +
     &                    coarse(c1-2, c2-1)  )
          fine(f10, f21) = tmp1
     &     -   F189o512*( coarse(c1-1, c2  )  )
     &     +   F315o512*( coarse(c1  , c2-1)  )
     &     +  F105o1024*( coarse(c1-2, c2  )  )
     &     -  F135o1024*( coarse(c1  , c2-2)  )
     &     +    F35o512*( coarse(c1-2, c2-1)  )
     &     +    F27o512*( coarse(c1-1, c2-2)  )
          fine(f11, f20) = tmp1
     &     -   F189o512*( coarse(c1  , c2-1)  )
     &     +   F315o512*( coarse(c1-1, c2  )  )
     &     +  F105o1024*( coarse(c1  , c2-2)  )
     &     -  F135o1024*( coarse(c1-2, c2  )  )
     &     +    F35o512*( coarse(c1-1, c2-2)  )
     &     +    F27o512*( coarse(c1-2, c2-1)  )
          fine(f10, f20) =
     &       F2025o1024*( coarse(c1  , c2  )  )
     &     +    F81o256*( coarse(c1-1, c2-1)  )
     &     +   F25o1024*( coarse(c1-2, c2-2)  )
     &     -   F405o512*( coarse(c1-1, c2  ) +
     &                    coarse(c1  , c2-1)  )
     &     +  F225o1024*( coarse(c1  , c2-2) +
     &                    coarse(c1-2, c2  )  )
     &     -    F45o512*( coarse(c1-1, c2-2) +
     &                    coarse(c1-2, c2-1)  )
c     write(*,1000) "43",c1,c2,fine(f11  ,f21  ),fine(f10  ,f21  ),
c    .                         fine(f11  ,f20  ),fine(f10  ,f20  )
        end if
      end if
c
      return
      end
