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:            if (dafm(i,j) * dafp(i,j) <= 0.0 .or. &
                  (a(i,j) - a(i-2,j))*(a(i+2,j) - a(i,j)) <= 0.0) then
  
                s = sign(1.0, d2ac(i,j))
                if (s == sign(1.0, d2ac(i-1,j)) .and. s == sign(1.0, d2ac(i+1,j)) .and. &
                     s == sign(1.0, d2af(i,j))) then
  
negative.code:   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
  