anchor.diff:  
        do j = jlo, jhi
           do i = ilo-1, ihi+2
  +           ! these live on cell-centers
              dafm(i,j) = a(i,j) - a_int(i,j)
              dafp(i,j) = a_int(i+1,j) - a(i,j)
   
  +           ! these live on cell-centers
              d2af(i,j) = 6.0*(a_int(i,j) - 2.0*a(i,j) + a_int(i+1,j))
              d2ac(i,j) = a(i-1,j) - 2.0*a(i,j) + a(i+1,j)
   
  -           d3a(i,j) = d2ac(i,j) - d2af(i-1,j)
  +           ! this lives on the interface
  +           d3a(i,j) = d2ac(i,j) - d2ac(i-1,j)
   
           enddo
        enddo
  
positive.code:               enddo
  
                if (rho < 1.0d0 - 1.d-12) then
                   ! we may need to limit
                   d3a_min = min(d3a(i-1,j), d3a(i,j), d3a(i+1,j), d3a(i+1,j))
                   d3a_max = max(d3a(i-1,j), d3a(i,j), d3a(i+1,j), d3a(i+1,j))
  
                   if (C3*max(abs(d3a_min), abs(d3a_max)) <= (d3a_max - d3a_min)) then
                      ! limit
  
negative.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
                         ! Eq. 32
                         al(i+1,j) = a(i,j) + 2.0d0*(1.0d0 - rho)*dafm(i,j) + rho*dafp(i,j)
                      endif
  
  