anchor.diff:                         ! 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
   
  
positive.code:                  d2a_lim = 0.0d0
                endif
  
                if (abs(d2af(i,j)) < 1.e-12*max(abs(a(i-2,j)), abs(a(i-1,j)), &
                                               abs(a(i,j)), abs(a(i+1,j)), abs(a(i+2,j)))) then
                   rho = 0.0
                else
                   ! MC Eq. 27
  
negative.code: 
       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)
  
             ! this lives on the interface
             d3a(i,j) = d2ac(i,j) - d2ac(i-1,j)
  
          enddo
       enddo
  