anchor.diff:                          do jj = 1,refratio
                              ifine = (i-1)*refratio + ii
                              jfine = (j-1)*refratio + jj
  -                           kf = areafine(ifine,jfine,k)
  +                           kf = areafine(ifine,jfine)
  +                           volf = kf*dz
                              qf = qfine(ifine,jfine,k,m)
  -                           sum = sum + kf*qf
  +                           sum = sum + volf*qf
                           enddo
                       enddo
   
  -                    kc = areacoarse(i+ic_add,j+jc_add,k)
  +                    kc = areacoarse(i+ic_add,j+jc_add)
  +                    volc = kc*dz
                       qc = qcoarse(i+ic_add, j+jc_add,k,m)
  -                    cons_diff = (qc*kc - sum)/r2
  +                    cons_diff = (qc*volc - sum)/r2
   
                       do ii = 1,refratio
                           do jj = 1,refratio
                              ifine  = (i-1)*refratio + ii
                              jfine  = (j-1)*refratio + jj
  -                           kf = areafine(ifine,jfine,k)
  -                           qfine(ifine,jfine,k,m) = qfine(ifine,jfine,k,m) + cons_diff/kf
  +                           kf = areafine(ifine,jfine)
  +                           volf = kf*dz
  +                           qfine(ifine,jfine,k,m) = qfine(ifine,jfine,k,m) + cons_diff/volf
                          end do
                       end do
                   end do !! j loop
  
positive.code:     double precision :: qcoarse(1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc,meqn)
      double precision ::   qfine(1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc,meqn)
  
      double precision :: areacoarse(-mbc:mx+mbc+1,-mbc:my+mbc+1,-mbc:mz+mbc+1)
      double precision ::   areafine(-mbc:mx+mbc+1,-mbc:my+mbc+1,-mbc:mz+mbc+1)
  
      integer :: ii, jj, i,j, i1, i2, j1, j2, ig, jg, mq, mth, k
      integer :: ic,jc,ic_add, jc_add, ifine, jfine
  
negative.code: 
      double precision ::  qcoarse(1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc,meqn)
      double precision ::    qfine(1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc,meqn)
      double precision ::  areacoarse(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision ::    areafine(-mbc:mx+mbc+1,-mbc:my+mbc+1)
  
      integer :: i,j,k,ii, jj, ifine, jfine, m, ig, jg, ic_add, jc_add
      double precision :: kf, kc, r2, sum, cons_diff, qf, qc, volf, dz, volc
  
      p4est_refineFactor = 2
      refratio = 2
  
      !! This is a bogus value, since we shouldn't end up here (yet). 
      dz = 1
  
      !! # Get (ig,jg) for grid from linear (igrid) coordinates
      ig = mod(igrid,refratio)
  