anchor.diff:        enddo 
        
   
  -
         do i=2-mbc, mx+mbc-1
            do j=2-mbc, my+mbc-1
               do ii = -1,1
  
positive.code: c     # Use cell-averaged values instead of pointwise values at 
  c     # last layers of ghost cells (where we cannot compute the
  c     # Laplacian.)
        do j = 1-mbc,my+mbc
            aux(1-mbc,j,1) = q(1-mbc,j,1)
            aux(mx+mbc,j,1) = q(mx+mbc,j,1)
  
negative.code:                   fp(ii,jj) = 0.5*qp**2
                  end do
              end do
  c           # Compute the average flux from pointwise values            
              flap = (fp(-1,0) + fp(1,0) + fp(0,-1) + 
       &                  fp(0,1) - 4*fp(0,0))
              aux(i,j,2) = fp(0,0) + flap/24.
  