c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
c Curvature based detection
c
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine cles_dcflag(ux,vx,dcflag,ncomps,nvars,
$ ixlo,ixhi, nx, dx, direction,extra)
implicit none
integer ncomps, nvars
integer ixlo, ixhi
integer direction, extra(*)
INTEGER nx
integer dcflag(1:nx+1,1)
double precision ux(ncomps,ixlo:ixhi)
double precision vx(nvars,ixlo:ixhi)
DOUBLE PRECISION dx
integer i, m, slow, shigh
double precision curvP, curvRho, cP, cRho
integer span
double precision CurvThressholdP, CurvThressholdRho
CurvThressholdP = 0.0d0
CurvThressholdRho = 0.0d0
span = 3
cP = CurvThressholdP *dx*dx * 0.25d0
cRho = CurvThressholdRho *dx*dx * 0.25d0
! ---- test criteria at each point
DO i = 1, nx+1, 1
curvP = (vx(5,i-1)-2.*vx(5,i)+vx(5,i+1))
curvP = curvP/(vx(5,i-1)+2.*vx(5,i)+vx(5,i+1))
curvrho = (vx(1,i-1)-2.*vx(1,i)+vx(1,i+1))
curvrho = curvrho/(vx(1,i-1)+2.*vx(1,i)+vx(1,i+1))
IF ((abs(curvP).gt. CurvThressholdP).and.
$ (abs(curvrho).gt.CurvThressholdRho)
$ .and.(curvP*curvrho.gt.0.0)) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(i-span,1)
shigh = min(i+span,nx+1)
do m=slow,shigh
dcflag(m,direction) = 1
enddo
END IF
END DO
return
end