anchor.diff:          do i = ilo-2, ihi+3
   
              ! interpolate to the edges
  -           a_int(i,j) = (7.0/12.0)*(a(i-1,j) + a(i,j)) - (1.0/12.0)*(a(i-2,j) + a(i+1,j))
  +           a_int(i,j) = (7.0d0/12.0d0)*(a(i-1,j) + a(i,j)) - (1.0d0/12.0d0)*(a(i-2,j) + a(i+1,j))
   
              al(i,j) = a_int(i,j)
              ar(i,j) = a_int(i,j)
  
positive.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.0*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
  
negative.code:         do i = ilo-1, ihi+1
  
             ! limit? MC Eq. 24 and 25
             if (dafm(i,j) * dafp(i,j) <= 0.0d0 .or. &
                  (a(i,j) - a(i,j-2))*(a(i,j+2) - a(i,j)) <= 0.0d0) then
  
                ! we are at an extrema
  
  