anchor.diff:    nx = qx - 2*ng; ny = qy - 2*ng
     ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1
   
  +  ! our convention here is that:
  +  !     al(i,j)   will be al_{i-1/2,j),
  +  !     al(i+1,j) will be al_{i+1/2,j)
   
     ! we need interface values on all faces of the domain
     if (idir == 1) then
  
positive.code:                        ! Eq. 31
                         ar(i,j) = a(i,j) - 2.0d0*(1.0d0 - rho)*dafp(i,j) - rho*dafm(i,j)
                      else if (abs(dafp(i,j)) >= 2.0*abs(dafm(i,j))) then
                         al(i+1,j) = a(i,j) + 2.0d0*(1.0d0 - rho)*dafm(i,j) + rho*dafp(i,j)
                      endif
  
  
negative.code:               enddo
  
                if (rho < 1.0d0 - 1.d-12) then
                   ! we may need to limit -- these quantities are at cell-centers
                   d3a_min = min(d3a(i-1,j), d3a(i,j), d3a(i+1,j), d3a(i+2,j))
                   d3a_max = max(d3a(i-1,j), d3a(i,j), d3a(i+1,j), d3a(i+2,j))
  
                   if (C3*max(abs(d3a_min), abs(d3a_max)) <= (d3a_max - d3a_min)) then
                      ! limit
  