anchor.diff:  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
  
positive.code: 
  end subroutine nolimit
  
  subroutine limit2(idir, a, lda, qx, qy, ng)
  
    ! 2nd order limited centered difference
  
negative.code: end subroutine limit2
  
  
  !-----------------------------------------------------------------------------
  ! limit4
  !-----------------------------------------------------------------------------
  subroutine limit4(idir, a, lda, qx, qy, ng)
  
    ! 4th order limited centered difference
  