anchor.diff:       &      icorner,iblock)
         implicit none
   
  -      integer mx, my, mbc, meqn, icorner, iblock, refratio,maux
  -      double precision areacoarse(-mbc:mx+mbc+1,1-mbc:my+mbc+1)
  -      double precision areafine(-mbc:mx+mbc+1,1-mbc:my+mbc+1)
  +      integer mx, my, mbc, meqn, icorner, iblock, refratio
  +      double precision areacoarse(-mbc:mx+mbc+1,-mbc:my+mbc+1)
  +      double precision areafine(-mbc:mx+mbc+1,-mbc:my+mbc+1)
  +
         double precision qcoarse(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)
         double precision qfine(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)
   
  -      integer mq, ibc, jbc, ii, jj, r2, ifine, jfine
  -      double precision sumfine,qf,areaf,areac
  -
  -      r2 = refratio**2
  +      integer mq, ibc, jbc, ii, jj, ifine, jfine
  +      double precision sum,qf,kf,kc
  +
  +
         do mq = 1,meqn
            if (icorner .eq. 0) then
               do ibc = 1,mbc
                  do jbc = 1,mbc
  -                  sumfine = 0
  +                  sum = 0
                     do ii = 1,refratio
                        do jj = 1,refratio
                           ifine = (ibc-1)*refratio + ii
                           jfine = (jbc-1)*refratio + jj
                           qf = qfine(ifine,jfine,mq)
  -                        areaf = areafine(ifine,jfine)
  -                        sumfine = sumfine + areaf*qf
  -                     enddo
  -                  enddo
  -                  areac = areacoarse(1-ibc,1-ibc)
  -                  qcoarse(1-ibc,1-jbc,mq) = sumfine/areac
  +                        kf = areafine(ifine,jfine)
  +                        sum = sum + kf*qf
  +                     enddo
  +                  enddo
  +                  kc = areacoarse(1-ibc,1-ibc)
  +                  qcoarse(1-ibc,1-jbc,mq) = sum/kc
                  enddo
               enddo
            elseif (icorner .eq. 1) then
               do ibc = 1,mbc
                  do jbc = 1,mbc
  -                  sumfine = 0
  +                  sum = 0
                     do ii = 1,refratio
                        do jj = 1,refratio
                           ifine = (ibc-1)*refratio + ii
                           jfine = (jbc-1)*refratio + jj
                           qf = qfine(mx+1-ifine,jfine,mq)
  -                        areaf = areafine(mx+1-ifine,jfine)
  -                        sumfine = sumfine + areaf*qf
  -                     enddo
  -                  enddo
  -                  areac = areacoarse(mx+ibc,1-jbc)
  -                  qcoarse(mx+ibc,1-jbc,mq) = sumfine/areac
  +                        kf = areafine(mx+1-ifine,jfine)
  +                        sum = sum + kf*qf
  +                     enddo
  +                  enddo
  +                  kc = areacoarse(mx+ibc,1-jbc)
  +                  qcoarse(mx+ibc,1-jbc,mq) = sum/kc
                  enddo
               enddo
            elseif (icorner .eq. 2) then
               do ibc = 1,mbc
                  do jbc = 1,mbc
  -                  sumfine = 0
  +                  sum = 0
                     do ii = 1,refratio
                        do jj = 1,refratio
                           ifine = (ibc-1)*refratio + ii
                           jfine = (jbc-1)*refratio + jj
                           qf = qfine(ifine,my+1-jfine,mq)
  -                        areaf = areafine(ifine,my+1-jfine)
  -                        sumfine = sumfine + areaf*qf
  -                     enddo
  -                  enddo
  -                  areac = areacoarse(1-ibc,my+jbc)
  -                  qcoarse(1-ibc,my+jbc,mq) = sumfine/areac
  +                        kf = areafine(ifine,my+1-jfine)
  +                        sum = sum + kf*qf
  +                     enddo
  +                  enddo
  +                  kc = areacoarse(1-ibc,my+jbc)
  +                  qcoarse(1-ibc,my+jbc,mq) = sum/kc
                  enddo
               enddo
            elseif (icorner .eq. 3) then
               do ibc = 1,mbc
                  do jbc = 1,mbc
  -                  sumfine = 0
  +                  sum = 0
                     do ii = 1,refratio
                        do jj = 1,refratio
                           ifine = (ibc-1)*refratio + ii
                           jfine = (jbc-1)*refratio + jj
                           qf = qfine(mx+1-ifine,my+1-jfine,mq)
  -                        areaf = areafine(mx+1-ifine,my+1-jfine)
  -                        sumfine = sumfine + areaf*qf
  -                     enddo
  -                  enddo
  -                  areac = areacoarse(mx+ibc,my+jbc)
  -                  qcoarse(mx+ibc,my+jbc,mq) = sumfine/areac
  +                        kf = areafine(mx+1-ifine,my+1-jfine)
  +                        sum = sum + kf*qf
  +                     enddo
  +                  enddo
  +                  kc = areacoarse(mx+ibc,my+jbc)
  +                  qcoarse(mx+ibc,my+jbc,mq) = sum/kc
                  enddo
               enddo
            endif
         enddo
  +
         end
   
   
  
positive.code:      &      icorner, is_block_bdry)
        implicit none
  
        integer mx,my,mbc,meqn,refratio,icorner, maux
        integer is_block_bdry(0:3)
        double precision areacoarse(-mbc:mx+mbc+1,-mbc:my+mbc+1)
        double precision   areafine(-mbc:mx+mbc+1,-mbc:my+mbc+1)
        double precision qcoarse(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)
        double precision qfine(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)
        double precision sumfine
  
        integer i,j,ibc,jbc,i1,j1,ii,jj,mq,r2
        integer ifine, jfine
        double precision areaf, areac, qf
  
        r2 = refratio*refratio
        do mq = 1,meqn
           if (icorner == 0) then
              do ibc = 1,mbc
                 do jbc = 1,mbc
  c                 # Average fine grid corners onto coarse grid ghost corners
                    sumfine = 0
                    do ii = 1,refratio
                       do jj = 1,refratio
                          ifine = (ibc-1)*refratio + ii
                          jfine = (jbc-1)*refratio + jj
                          if (is_block_bdry(0) .eq. 1) then
                             areaf = areafine(ifine,my+1-jfine)
                             qf =   qfine(ifine,my+1-jfine,mq)
                          elseif (is_block_bdry(2) .eq. 1) then
                             areaf = areafine(mx+1-ifine,jfine)
                             qf =   qfine(mx+1-ifine,jfine,mq)
                          endif
                          sumfine = sumfine + areaf*qf
                       enddo
                    enddo
                    areac = areacoarse(1-ibc,1-jbc)
                    qcoarse(1-ibc,1-jbc,mq) = sumfine/areac
                 enddo
              enddo
           elseif (icorner .eq. 1) then
              do ibc = 1,mbc
                 do jbc = 1,mbc
  c                 # Average fine grid corners onto coarse grid ghost corners
                    sumfine = 0
                    do ii = 1,refratio
                       do jj = 1,refratio
                          ifine = (ibc-1)*refratio + ii
                          jfine = (jbc-1)*refratio + jj
                          if (is_block_bdry(1) .eq. 1) then
                             areaf = areafine(mx+1-ifine,my+1-jfine)
                             qf =   qfine(mx+1-ifine,my+1-jfine,mq)
                          elseif (is_block_bdry(2) .eq. 1) then
                             areaf = areafine(ifine,jfine)
                             qf =   qfine(ifine,jfine,mq)
                          endif
                          sumfine = sumfine + areaf*qf
                       enddo
                    enddo
                    areac = areacoarse(mx+ibc,1-jbc)
                    qcoarse(mx+ibc,1-jbc,mq) = sumfine/areac
                 enddo
              enddo
           elseif (icorner .eq. 2) then
              do ibc = 1,mbc
                 do jbc = 1,mbc
  c                 # Average fine grid corners onto coarse grid ghost corners
                    sumfine = 0
                    do ii = 1,refratio
                       do jj = 1,refratio
                          ifine = (ibc-1)*refratio + ii
                          jfine = (jbc-1)*refratio + jj
                          if (is_block_bdry(0) .eq. 1) then
                             areaf = areafine(ifine,jfine)
                             qf =   qfine(ifine,jfine,mq)
                          elseif (is_block_bdry(3) .eq. 1) then
                             areaf = areafine(mx+1-ifine,my+1-jfine)
                             qf =   qfine(mx+1-ifine,my+1-jfine,mq)
                          endif
                          sumfine = sumfine + areaf*qf
                       enddo
                    enddo
                    areac = areacoarse(1-ibc,my+jbc)
                    qcoarse(1-ibc,my+jbc,mq) = sumfine/areac
                 enddo
              enddo
           elseif (icorner .eq. 3) then
              do ibc = 1,mbc
                 do jbc = 1,mbc
  c                 # Average fine grid corners onto coarse grid ghost corners
                    sumfine = 0
                    do ii = 1,refratio
                       do jj = 1,refratio
                          ifine = (ibc-1)*refratio + ii
                          jfine = (jbc-1)*refratio + jj
                          if (is_block_bdry(1) .eq. 1) then
                             areaf = areafine(mx+1-ifine,jfine)
                             qf =   qfine(mx+1-ifine,jfine,mq)
                          elseif (is_block_bdry(3) .eq. 1) then
                             areaf = areafine(ifine,my+1-jfine)
                             qf =   qfine(ifine,my+1-jfine,mq)
                          endif
                          sumfine = sumfine + areaf*qf
                       enddo
                    enddo
                    areac = areacoarse(mx+ibc,my+jbc)
                    qcoarse(mx+ibc,my+jbc,mq) = sumfine/areac
                 enddo
              enddo
           endif
        enddo
  
        end
  
  
negative.code:             do ibc = 1,mbc
                 do jbc = 1,mbc
                    if (is_block_bdry(0) .eq. 1) then
                       qthis(1-ibc,my+jbc,mq) =
       &                     qneighbor(ibc,jbc,mq)
                       qneighbor(1-ibc,1-jbc,mq) =
       &                     qthis(ibc,my+1-jbc,mq)
  