anchor.diff:  end subroutine limit2
   
   
  +!-----------------------------------------------------------------------------
  +! limit4
  +!-----------------------------------------------------------------------------
   subroutine limit4(idir, a, lda, qx, qy, ng)
   
     ! 4th order limited centered difference
  
positive.code: 
  end subroutine nolimit
  
  subroutine limit2(idir, a, lda, qx, qy, ng)
  
    ! 2nd order limited centered difference
  
negative.code: end subroutine limit4
  
  
  !-----------------------------------------------------------------------------
  ! flatten
  !-----------------------------------------------------------------------------
  subroutine flatten(idir, p, u, xi, qx, qy, ng, smallp, delta, z0, z1)
  
    ! 1-d flattening near shocks
    !
    ! flattening kicks in behind strong shocks and reduces the
    ! reconstruction to using piecewise constant slopes, making things
    ! first-order.  See Saltzman (1994) page 159 for this
    ! implementation.
    implicit none
  
    integer, intent(in) :: idir
    integer, intent(in) :: qx, qy, ng
  
    ! 0-based indexing to match python
    double precision, intent(inout) :: p(0:qx-1, 0:qy-1), u(0:qx-1, 0:qy-1)
    double precision, intent(  out) :: xi(0:qx-1, 0:qy-1)
    double precision, intent(in   ) :: smallp, delta, z0, z1
  
  !f2py intent(in) :: idir
  !f2py depend(qx, qy) :: p, u, xi
  !f2py intent(in) :: p, u
  !f2py intent(out) :: xi
  
  ! call this as: xi = reconstruction_f.flatten(1,p,u,qx,qy,ng,smallp,delta,z0,z1)
  
    integer :: ilo, ihi, jlo, jhi
    integer :: nx, ny
    integer :: i, j
  
    double precision :: test1, test2
    double precision :: dp, dp2, z
  
    nx = qx - 2*ng; ny = qy - 2*ng
    ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1
  
    xi(:,:) = 1.0d0
  
    select case (idir)
  
    case (1)
  
       do j = jlo-2, jhi+2
          do i = ilo-2, ihi+2
  
             dp  = abs(p(i+1,j) - p(i-1,j))
             dp2 = abs(p(i+2,j) - p(i-2,j))
             z = dp/max(dp2,smallp)
  
             test1 = u(i-1,j) - u(i+1,j)
             test2 = dp/min(p(i+1,j),p(i-1,j))
  
             if (test1 > 0.0d0 .and. test2 > delta) then
                xi(i,j) = min(1.0d0, max(0.0d0, 1.0d0 - (z - z0)/(z1 - z0)))
             else
                xi(i,j) = 1.0
             endif
  
          enddo
       enddo
  
    case (2)
  
       do j = jlo-2, jhi+2
          do i = ilo-2, ihi+2
  
             dp  = abs(p(i,j+1) - p(i,j-1))
             dp2 = abs(p(i,j+2) - p(i,j-2))
             z = dp/max(dp2,smallp)
  
             test1 = u(i,j-1) - u(i,j+1)
             test2 = dp/min(p(i,j+1),p(i,j-1))
  
             if (test1 > 0.0d0 .and. test2 > delta) then
                xi(i,j) = min(1.0d0, max(0.0d0, 1.0d0 - (z - z0)/(z1 - z0)))
             else
                xi(i,j) = 1.0
             endif
  
          enddo
       enddo
  
    end select
  
  end subroutine flatten
  
      
  !-----------------------------------------------------------------------------
  ! flatten_multid
  !-----------------------------------------------------------------------------
  subroutine flatten_multid(xi_x, xi_y, p, xi, qx, qy, ng)
  
    ! multi-dimensional flattening
    implicit none
  
    integer, intent(in) :: qx, qy, ng
  
    ! 0-based indexing to match python
    double precision, intent(inout) :: xi_x(0:qx-1, 0:qy-1)
    double precision, intent(inout) :: xi_y(0:qx-1, 0:qy-1)
    double precision, intent(inout) :: p(0:qx-1, 0:qy-1)
    double precision, intent(  out) :: xi(0:qx-1, 0:qy-1)
  
  !f2py depend(qx, qy) :: xi_x, xi_y, p, xi
  !f2py intent(in) :: xi_x, xi_y, p
  !f2py intent(out) :: xi
  
  ! call this as: xi = reconstruction_f.flatten_multid(xi_x,xi_y,,p,qx,qy,ng)
  
    integer :: ilo, ihi, jlo, jhi
    integer :: nx, ny
    integer :: i, j
  
    integer :: sx, sy
  
    nx = qx - 2*ng; ny = qy - 2*ng
    ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1
  
    xi(:,:) = 1.0d0
  
    do j = jlo-2, jhi+2
       do i = ilo-2, ihi+2
  
          sx = int(sign(1.0d0, p(i+1,j) - p(i-1,j)))
          sy = int(sign(1.0d0, p(i,j+1) - p(i,j-1)))
  
          xi(i,j) = min(min(xi_x(i,j), xi_x(i-sx,j)), &
                        min(xi_y(i,j), xi_y(i,j-sy)))
  
       enddo
    enddo
  
  end subroutine flatten_multid
  