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
include 'cles.i'
include 'cuser.i'
integer ncomps, nvars
integer ixlo, ixhi
integer nx, direction, extra(*)
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 sum, sumth, sumth2, sumy, sumyt, th
double precision x, d, det, a, b, c, dp
integer span, n, imin, imax
parameter (span=3)
double precision thvec(-span:span)
c double precision minc, minjump, dwidth
c thickness of the profile in terms of grid sizes
c a value of 1 means a discontinuity of 6 grids points
c because of the thickness of the tanh()
c dwidth = 0.75d0
c minimum value of a correlation coefficient considered
c to be a good match
c minc = 0.9d0
c to discriminate small oscillations in the density that
c have no meaning, we stipulate that any jump in density
c that is below 2% of the mean value is not a jump.
c minjump = 0.2d-1
c ---- test criteria at each point
sumth = 0.0d0
sumth2 = 0.0d0
do m=-span,span, 1
x = m/dwidth
th = tanh(x)
thvec(m) = th
sumth = sumth + th
sumth2 = sumth2 + th*th
enddo
n = span*2+1
imin = max(ixlo+span,1)
imax = min(ixhi-span,nx)
DO i = imin, imax, 1
sumy = 0.0d0
sumyt = 0.0d0
do m=-span,span, 1
sumy = sumy + vx(1,i+m)
sumyt = sumyt + vx(1,i+m)*thvec(m)
enddo
c Optimal values
det = n*sumth2-sumth*sumth
a = (sumy*sumth2-sumth*sumyt)/det
b = (n*sumyt-sumth*sumy)/det
c some tests conditions
if ( a .lt. 0.0d0 ) cycle
if ( a + b .lt. 0.0d0 ) cycle
if ( a - b .lt. 0.0d0 ) cycle
if ( abs(b) .lt. minjump*a ) cycle
c Compute correlation
sum = 0.0d0
sumy = 0.0d0
do m=-span,span, 1
dp = vx(1,i+m)-a
sum = sum + dp*thvec(m)
sumy = sumy + dp*dp
enddo
if ( abs(sumy) .lt. 1.0d-12 ) cycle
c = sum/sqrt(sumy*sumth2)
if ( abs(c) .gt. minc ) then
slow = max(i-span,1)
shigh = min(i+1+span,nx+1)
do m=slow,shigh
dcflag(m,direction) = CLES_SWITCH_WENO
enddo
endif
enddo
return
end