!rho:  0.10000000E-02  0.10000000E-02  0.14988374E-02  0.14988374E-02  0.14988374E-02  0.14988374E-02  0.10000000E-02
!u:    0.00000000E+00  0.00000000E+00  0.00000000E+00 -0.32206959E+04  0.32206959E+04  0.00000000E+00  0.00000000E+00
!p:    0.84087411E+06  0.84087411E+06  0.26033973E+07  0.26033973E+07  0.26033973E+07  0.26033973E+07  0.84087411E+06
                                                               


!     flaten zone structure in regions where shocks are too thin

    implicit none
    integer, parameter :: q=16

    integer:: nzn, ishkbn      
    real, DIMENSION(q) :: u, p,flatn,flatn1,dp,du,shockd,shockd1
    real :: epsiln, omg1, omg2,smallu,smallp,igodu

    integer :: i, nzn5, nzn6, nzn7, nzn8
    real :: shkbrn, utest, dutest, dp2, dptest, ptest
    real, DIMENSION(q) :: scrch1(q), scrch2(q), scrch3(q)
    real :: sig1,sig2,ak1,ak2,wig1,wig2,wig3
    logical :: repeat = .FALSE.

u=0.d0
u(1) = 0.00000000E+00
u(2) =  0.00000000E+00
u(3) =  0.00000000E+00
u(4) = -0.32206959E+04
u(5) =  0.32206959E+04
u(6) =  0.00000000E+00
u(7) =  0.00000000E+00

p =  0.84087411E+06
p(1) = 0.84087411E+06
p(2) =  0.84087411E+06
p(3) =  0.26033973E+07
p(4) =  0.26033973E+07
p(5) =  0.26033973E+07
p(6) =  0.26033973E+07 
p(7) = 0.84087411E+06
     
10 continue

           epsiln = 0.33                                             
                                                                
ishkbn =0.

          igodu=0.d0
          smallu = 1.d-10
          smallp =1.d-20
                                                
          omg1   = 0.75                                             
          omg2   = 10.0                                             
          sig1   = 0.50                                             
          sig2   = 1.00                                             
          ak1    = 2.00                                             
          ak2    = 0.01                                             
                                                                     
          wig1   = 2.00                                             
          wig2   = 0.00                             
          wig3   = 0.3333 - wig2                                    

    nzn = 8

    nzn5 = nzn + 5
    nzn6 = nzn + 6
    nzn7 = nzn + 7
    nzn8 = nzn + 8

    do i = 1, nzn8
       flatn (i) = 0.0e00
       shockd(i) = 0.0e00
    end do

    do i = 2, nzn7
       shkbrn     = ishkbn

! compute the w_j parameter in Eq. A.1 in Colella & Woodward.  w_j is equal to
! 1 if the jth zone is inside a pressure and velocity jump in the sweep direction,
! in a manner consistent with a shock -- store it in shockd1

       dp(i)      = p(i+1) - p(i-1)
       du(i)      = u(i+1) - u(i-1)
       scrch1(i)  = epsiln * min (p(i+1), p(i-1)) - abs( dp(i) )
       utest      = smallu - abs (du(i))

       if (utest .LT. 0.e0) then
          dutest = du(i)
       else
          dutest = 0.e0
       endif

       if (scrch1(i) .LT. 0.e0) then 
          scrch1(i) = 1.e0
       else
          scrch1(i) = 0.e0
       endif

       if (du(i) .GE. 0.e0) scrch1(i) = 0.e0

       if (dutest .EQ. 0.e0) scrch1(i) = 0.e0

       if (shkbrn .NE. 0.e0) then 
          shockd1(i) = 0.e0
       else
          shockd1(i) = scrch1(i)
       endif

    end do

    do i = 3, nzn6
       shockd(i) = max (shockd1(i-1), shockd1(i), shockd1(i+1))
    end do


! compute the dissipative flux, using Eq. A.2 in Colella & Woodward

    do i = 3, nzn6
       dp2 = p(i+2) - p(i-2)

       if (dp2 .EQ. 0.e0) dp2 = smallp

       dptest = dp(i)
       ptest = dp2 - smallp

       if (ptest .EQ. 0.e0) dptest = 0.e0

       scrch2(i) = dptest / dp2 - omg1
       scrch3(i) = scrch1(i) * max (0.e00, scrch2(i) * omg2)
    end do

    do i = 2, nzn7
!
!
! this is the at.hoc modification
       if (dp(i) .LT. 0.e0) scrch2(i) = scrch3(i+1)
       if (dp(i) .GT. 0.e0) scrch2(i) = scrch3(i-1)
       if (dp(i) .eq. 0.e0) scrch2(i) = scrch3(i  )
!
!this is the original
!!$       if (dp(i) .LT. 0.e0) then
!!$          scrch2(i) = scrch3(i+1)
!!$       else
!!$          scrch2(i) = scrch3(i-1)
!!$       endif
!
    end do

    do i = 1, nzn6
       flatn(i) = max (scrch3(i), scrch2(i))
       flatn(i) = max (0.0e00, min (1.0e00, flatn(i)))
    end do

    do i = 1, nzn8
       flatn (i) = flatn(i) * (1.0e00 - igodu) + igodu
       flatn1(i) = 1.0e00 - flatn(i)
    end do

do i=1,16
write(6,fmt='(I,7D10.3)') i,flatn(i),scrch1(i),scrch2(i),scrch3(i),dp(i),du(i),p(i)
end do

p(:)=p(16:1:-1)
u(:)=-u(16:1:-1)

if(.not.repeat) THEN
 repeat = .TRUE.
 write(6,*)
 goto 10
end if


  end