c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
c Curvature based detection
c
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine cles_dcflag(ux,vx,dcflag,ncomps,nvars,
$ ixlo,ixhi,iylo,iyhi, nx, ny, xl, xr, yl, yr, dx, dy,
$ direction,extra)
implicit none
integer ncomps, nvars
integer ixlo, ixhi, iylo, iyhi
integer direction, extra(*)
INTEGER nx,ny
integer dcflag(1:nx+1,1:ny+1,2)
double precision ux(ncomps,ixlo:ixhi,iylo:iyhi)
double precision vx(nvars,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION xl,xr,yl,yr
DOUBLE PRECISION dx,dy
integer i, j, m, slow, shigh
double precision curvP, curvRho, cP, cRho
integer span
double precision CurvThressholdP, CurvThressholdRho
! ---- initialize dcflag to zero ( eg, smooth everywhere)
do j=1,ny+1
do i=1,nx+1
dcflag(i,j,direction) = 0
enddo
enddo
CurvThressholdP = 0.015d0
CurvThressholdRho = 0.01d0
span = 4
! ---- test criteria at each point
if ( direction .eq. 1 ) then
cP = CurvThressholdP * dx * dx * 0.25d0
cRho = CurvThressholdRho * dx * dx * 0.25d0
do j = 1, ny, 1
do i = 1, nx+1, 1
curvP = (vx(5,i-1,j)-2.*vx(5,i,j)+vx(5,i+1,j))
curvP = curvP/(vx(5,i-1,j)+2.*vx(5,i,j)+vx(5,i+1,j))
curvrho = (vx(1,i-1,j)-2.*vx(1,i,j)+vx(1,i+1,j))
curvrho = curvrho/(vx(1,i-1,j)+2.*vx(1,i,j)+vx(1,i+1,j))
IF ((abs(curvP).gt. cP).and.(abs(curvrho).gt.cRho)
$ .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,j,direction) = 1
enddo
END IF
enddo
enddo
else
cP = CurvThressholdP * dy * dy * 0.25d0
cRho = CurvThressholdRho * dy * dy * 0.25d0
do j = 1, ny+1, 1
DO i = 1, nx, 1
curvP = (vx(5,i,j-1)-2.*vx(5,i,j)+vx(5,i,j+1))
curvP = curvP/(vx(5,i,j-1)+2.*vx(5,i,j)+vx(5,i,j+1))
curvrho = (vx(1,i,j-1)-2.*vx(1,i,j)+vx(1,i,j+1))
curvrho = curvrho/(vx(1,i,j-1)+2.*vx(1,i,j)+vx(1,i,j+1))
IF ((curvP.gt. cP).and.(curvrho.gt.cRho)
$ .and.(curvP*curvrho.gt.0.0)) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(j-span,1)
shigh = min(j+span,ny+1)
do m=slow,shigh
dcflag(i,m,direction) = 1
enddo
END IF
enddo
enddo
endif
return
end