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