anchor.diff:        double precision favg(-mbc:mx+mbc+1,-mbc:my+mbc+1)
         integer compute_avg, ghost_only
   
  -      deg = 4
  +      deg = 2
         call integrate_exact(mx,my,mbc,xlower,ylower,dx,dy,
        &      blockno,deg,area,f,favg,compute_avg,ghost_only)
   
  
positive.code:       double precision favg(-mbc:mx+mbc+1,-mbc:my+mbc+1)
  
        integer i,j,ii,jj, k
        double precision xlow, ylow, a,b,c,d,dx,dy, area1
        double precision quad_square, quad_square_area
  
        do i = -mbc,mx+mbc+1
           do j = -mbc,my+mbc+1
  
negative.code: c     # ------------------------------------------------------
  c     # Misc. functions needed for all quad calculations
  c     # ------------------------------------------------------
        double precision function quad_square_favg(f,a,b,c,d,deg,
       &      blockno,area_phys)
        implicit none
  
        double precision f
        integer deg
        integer blockno
        double precision a,b, c,d, area_phys
        double precision x(8),wx(8), xc,yc
        double precision y(8),wy(8)
        double precision area_element_exact, ae, h
        double precision area_phys1, area_phys2
        double precision area_comp1, area_comp2
        double precision favg1, favg2
        integer i,j
  
        call gltable2(x,wx,deg,a,b)
        call gltable2(y,wy,deg,c,d)
  
  c     # area_comp = computational area (b-a)*(d-c) of mesh cell
  c     # area_phys = physical area of mesh cell
  
        favg2 = 0
        area_comp2 = 0
        area_phys2 = 0
        do i = 1,deg
           favg1 = 0
           area_comp1 = 0
           area_phys1 = 0
           xc = x(i)
           do j = 1,deg
              yc = y(j)
              h = wy(j)
              area_comp1 = area_comp1 + h
              ae = area_element_exact(xc,yc,blockno)
              area_phys1 = area_phys1 + ae*h
              favg1 = favg1 + f(xc,yc)*ae*h
           enddo
           area_comp2 = area_comp2 + wx(i)*area_comp1
           area_phys2 = area_phys2 + wx(i)*area_phys1
           favg2 = favg2 + wx(i)*favg1
        enddo
  
  c     # Note that this is the average in computational space.
        quad_square_favg = favg2/area_comp2
        area_phys = area_phys2
  
        end
  
  