anchor.diff:             dafp(i,j) = a_int(i,j+1) - 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,j+1))
  +           d2af(i,j) = 6.0d0*(a_int(i,j) - 2.0d0*a(i,j) + a_int(i,j+1))
           enddo
        enddo
   
        do j = jlo-3, jhi+3
           do i = ilo-1, ihi+1
  -           d2ac(i,j) = a(i,j-1) - 2.0*a(i,j) + a(i,j+1)
  +           d2ac(i,j) = a(i,j-1) - 2.0d0*a(i,j) + a(i,j+1)
           enddo
        enddo
   
  
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
  
             ! interpolate to the edges
             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)
  