anchor.diff:          do i = ilo-1, ihi+1
   
              ! interpolate to the edges
  -           a_int(i,j) = (7.0/12.0)*(a(i,j-1) + a(i,j)) - (1.0/12.0)*(a(i,j-2) + a(i,j+1))
  +           a_int(i,j) = (7.0d0/12.0d0)*(a(i,j-1) + a(i,j)) - (1.0d0/12.0d0)*(a(i,j-2) + a(i,j+1))
   
              al(i,j) = a_int(i,j)
              ar(i,j) = a_int(i,j)
  
positive.code: 
                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
                   rho = d2a_lim/d2af(i,j)
  
negative.code:                        ! Eqs. 29, 30
                         ar(i,j) = a(i,j) - rho*dafm(i,j)  ! note: typo in Eq 29
                         al(i,j+1) = a(i,j) + rho*dafp(i,j)
                      else if (abs(dafm(i,j)) >= 2.0d0*abs(dafp(i,j))) then
                         ! 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
  