c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx c c Flux smoothness based detection c c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine cles_dcflag(ux,vx,dcflag,ncomps,nvars, $ ixlo,ixhi,iylo,iyhi, nx, ny, 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 dx,dy integer i, j, m, slow, shigh, loop, size, enoOrder double precision smoothness, dmdx, maxv(4), minv(4), variance(4) integer span double precision TriggerValue, smooth_eps TriggerValue = 0.0d0 smooth_eps = 0.001 span = 3 enoOrder = extra(1) size = enoOrder-1 ! ---- test criteria at each point if ( direction .eq. 1 ) then do j = 1, ny, 1 do i = 1, nx+1, 1 DO loop =1,4 do m=i-size, i+size-1 dmdx = abs(fx(loop,m+1,j) $ -fx(loop,m-1,j)) if ( m .eq. i-size ) then maxv(loop) = dmdx minv(loop) = dmdx else maxv(loop) = MAX(maxv(loop), dmdx) minv(loop) = MIN(minv(loop), dmdx) endif enddo variance(loop) = maxv(loop)/(minv(loop)+smooth_eps) if ( loop .eq. 1 ) then smoothness = variance(loop) else smoothness = MAX(smoothness, variance(loop)) endif END DO IF ( smoothness .gt. TriggerValue ) 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 do j = 1, ny+1, 1 DO i = 1, nx, 1 DO loop =1,4 do m=j-size, j+size-1 dmdx = abs(fx(loop,i,m+1) $ -fx(loop,i,m-1)) if ( m .eq. j-size ) then maxv(loop) = dmdx minv(loop) = dmdx else maxv(loop) = MAX(maxv(loop), dmdx) minv(loop) = MIN(minv(loop), dmdx) endif enddo variance(loop) = maxv(loop)/(minv(loop)+smooth_eps) if ( loop .eq. 1 ) then smoothness = variance(loop) else smoothness = MAX(smoothness, variance(loop)) endif END DO IF ( smoothness .gt. TriggerValue ) 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