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 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,dy integer i, m, slow, shigh double precision crvP, crvRho, crvS, x, y double precision ThressCrvP,ThressCrvRho,ThressCrvS integer dspan dspan = 4 ! ---- test criteria at each point if ( direction .eq. 1 ) then ThressCrvP = CurvP *dx*dx * 0.25d0 ThressCrvRho = CurvRho *dx*dx * 0.25d0 ThressCrvS = CurvS *dx*dx * 0.25d0 do i = 1, nx+1, 1 crvP = (vx(5,i-1)-2.*vx(5,i)+vx(5,i+1)) crvP = crvP/(vx(5,i-1)+2.*vx(5,i)+vx(5,i+1)) crvRho = (vx(1,i-1)-2.*vx(1,i)+vx(1,i+1)) crvRho = crvRho/(vx(1,i-1)+2.*vx(1,i)+vx(1,i+1)) crvS = (vx(6,i+1)-vx(6,i-1))/2.0d0 IF (((abs(crvP).gt. ThressCrvP).and. $ (abs(crvrho).gt.ThressCrvRho) $ .and.(crvP*crvRho.gt.0.0))) THEN ! we are at a shock ! ---- we will use WENO here slow = max(i-dspan,1) shigh = min(i+dspan,nx+1) do m=slow,shigh dcflag(m,direction) = CLES_SWITCH_WENO enddo END IF IF (abs(crvS).gt.CurvS) THEN ! we are at a shock ! ---- we will use WENO here slow = max(i-dspan,1) shigh = min(i+dspan,nx+1) do m=slow,shigh dcflag(m,direction) = CLES_SWITCH_WENO enddo END IF ! Flag the boundaries call cles_xLocation(i,x) if ( x.lt.(GeomMin(1)+2.0d0*dx) ) then if ( trunc .ne. 1 ) then do m=1, i+1 dcflag(m,direction) = CLES_SWITCH_WENO enddo endif else if (x.gt.(GeomMax(1)-2.0d0*dx)) then do m=i,nx+1 dcflag(m,direction) = CLES_SWITCH_WENO enddo endif enddo end if return end