anchor.diff:        integer deg
         integer blockno
         double precision a,b, c,d
  -      double precision x(8),wx(8), xc,yc,qn
  -      double precision y(8),wy(8), qn2,fc
  -      double precision fint2,fint1,area2,h
  -      double precision area_element_exact, ae
  +      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
         integer i,j
   
         call gltable2(x,wx,deg,a,b)
         call gltable2(y,wy,deg,c,d)
   
  -      fint2 = 0
  +      area_phys2 = 0
         do i = 1,deg
            xc = x(i)
  -         fint1 = 0
  +         area_phys1 = 0
            do j = 1,deg
               yc = y(j)
               h = wy(j)
               ae = area_element_exact(xc,yc,blockno)
  -            fint1 = fint1 + ae*h
  +            area_phys1 = area_phys1 + ae*h
            enddo
  -         fint2 = fint2 + wx(i)*fint1
  -      enddo
  -
  -      quad_square_area = fint2
  +         area_phys2 = area_phys2 + wx(i)*area_phys1
  +      enddo
  +
  +      quad_square_area = area_phys2
   
         end
   
  
positive.code:             c = ylow
              d = ylow + dy
              if (compute_avg .eq. 1) then
                 favg(i,j) = quad_square(f,a,b,c,d,deg,
       &               blockno,area1,compute_avg)
                 area(i,j) = area1
              else
                 area(i,j) = quad_square_area(a,b,c,d,deg,
       &               blockno)
  
negative.code:       double precision favg(-mbc:mx+mbc+1,-mbc:my+mbc+1)
        integer compute_avg, ghost_only
  
        deg = 2
        call integrate_exact(mx,my,mbc,xlower,ylower,dx,dy,
       &      blockno,deg,area,f,favg,compute_avg,ghost_only)
  
  