!================================================================
!                      v0(j),var(j)
! ------------#----------   ----- j ----------------#------------
!      fm(j-1)|fp(j-1)               fm(j)|fp(j)
!           flux(j-1)                 flux(j)
!
! fp(j): ---(-3)---(-2)---(-1)---(0)---(1)---(2)---(3)---(4)---
!         1: |--------------------|
!                 2: |------------------|
!                       3: |------------------|
!                              4: |------------------|
!
! fm(j): ---(-3)---(-2)---(-1)---(0)---(1)---(2)---(3)---(4)---
!                                    1: |-----------------|
!                              2: |-----------------|
!                       3: |------------------|
!                 4: |------------------|
!================================================================


subroutine wenozn7_p( fps, fpx )
    use commondata
    implicit none
    
    real*8 :: fps(-3:3), fpx
    real*8 :: omega(4), fu(4), beta(4)
    real*8 :: tao, tao7, tao12, ac, aa
    real*8 :: totalpha, alpha(4), fomega(4)
    
    omega(1) = 1.d0/35.d0;   omega(2) = 12.d0/35.d0
    omega(3) = 18.d0/35.d0;  omega(4) = 4.d0/35.d0
    
    fu(1) = -1.d0/4.d0*fps(-3) + 13.d0/12.d0*fps(-2) - 23.d0/12.d0*fps(-1) + 25.d0/12.d0*fps(0)
    fu(2) = 1.d0/12.d0*fps(-2) - 5.d0/12.d0*fps(-1) + 13.d0/12.d0*fps(0) + 1.d0/4.d0*fps(1)
    fu(3) = -1.d0/12.d0*fps(-1) + 7.d0/12.d0*fps(0) + 7.d0/12.d0*fps(1) - 1.d0/12.d0*fps(2)
    fu(4) = 1.d0/4.d0*fps(0) + 13.d0/12.d0*fps(1) - 5.d0/12.d0*fps(2) + 1.d0/12.d0*fps(3)
    
    beta(1) = fps(-3)*( 547.d0*fps(-3) - 3882.d0*fps(-2) + 4642.d0*fps(-1) - 1854.d0*fps(0) ) &
            + fps(-2)*( 7043.d0*fps(-2) - 17246.d0*fps(-1) + 7042.d0*fps(0) ) &
            + fps(-1)*( 11003.d0*fps(-1) - 9402.d0*fps(0) ) + 2107.d0*fps(0)**2
    beta(2) = fps(-2)*( 267.d0*fps(-2) - 1642.d0*fps(-1) + 1602.d0*fps(0) - 494.d0*fps(1) ) &
            + fps(-1)*( 2843.d0*fps(-1) - 5966.d0*fps(0) + 1922.d0*fps(1) ) &
            + fps(0)*( 3443.d0*fps(0) - 2522.d0*fps(1) ) + 547.d0*fps(1)**2
    beta(3) = fps(-1)*( 547.d0*fps(-1) - 2522.d0*fps(0) + 1922.d0*fps(1) - 494.d0*fps(2) ) &
            + fps(0)*( 3443.d0*fps(0) - 5966.d0*fps(1) + 1602.d0*fps(2) ) &
            + fps(1)*( 2843.d0*fps(1) - 1642.d0*fps(2) ) + 267.d0*fps(2)**2
    beta(4) = fps(0)*( 2107.d0*fps(0) - 9402.d0*fps(1) + 7042.d0*fps(2) - 1854.d0*fps(3) ) &
            + fps(1)*( 11003.d0*fps(1) - 17246.d0*fps(2) + 4642.d0*fps(3) ) &
            + fps(2)*( 7043.d0*fps(2) - 3882.d0*fps(3) ) + 547.d0*fps(3)**2
    beta(:) = beta(:)/1000.d0
    
    tao7 = abs(beta(1) - beta(4))
    tao = tao7
    aa = 100.d0
    ac = aa*( ( beta(1) + beta(4) - tao + eps )/( tao + eps ) )**2
    tao12 = ( fps(-3) - 6.d0*fps(-2) + 15.d0*fps(-1) - 20.d0*fps(0) &
            + 15.d0*fps(1) - 6.d0*fps(2) + fps(3) )**2
    
    alpha(:) = omega(:)*( ac + (tao12/(beta(:)+eps))**2 )
    !alpha(:) = omega(:)*( 1.d0 + (tao12/(beta(:)+eps))**2 )
    totalpha = alpha(1) + alpha(2) + alpha(3)+ alpha(4)
    fomega(:) = alpha(:)/totalpha
    fpx = fu(1)*fomega(1) + fu(2)*fomega(2) + fu(3)*fomega(3)+ fu(4)*fomega(4)
    
    return
    end subroutine wenozn7_p


subroutine wenozn7_m( fms, fmx )
    use commondata
    implicit none
    
    real*8 :: fms(-2:4), fmx
    real*8 :: omega(4), fu(4), beta(4)
    real*8 :: tao, tao7, tao12, ac, aa
    real*8 :: totalpha, alpha(4), fomega(4)
    
    omega(1) = 1.d0/35.d0;   omega(2) = 12.d0/35.d0
    omega(3) = 18.d0/35.d0;  omega(4) = 4.d0/35.d0
    
    fu(1) = -1.d0/4.d0*fms(4) + 13.d0/12.d0*fms(3) - 23.d0/12.d0*fms(2) + 25.d0/12.d0*fms(1)
    fu(2) = 1.d0/12.d0*fms(3) - 5.d0/12.d0*fms(2) + 13.d0/12.d0*fms(1) + 1.d0/4.d0*fms(0)
    fu(3) = -1.d0/12.d0*fms(2) + 7.d0/12.d0*fms(1) + 7.d0/12.d0*fms(0) - 1.d0/12.d0*fms(-1)
    fu(4) = 1.d0/4.d0*fms(1) + 13.d0/12.d0*fms(0) - 5.d0/12.d0*fms(-1) + 1.d0/12.d0*fms(-2)
    
    beta(1) = fms(4)*( 547.d0*fms(4) - 3882.d0*fms(3) + 4642.d0*fms(2) - 1854.d0*fms(1) ) &
            + fms(3)*( 7043.d0*fms(3) - 17246.d0*fms(2) + 7042.d0*fms(1) ) &
            + fms(2)*( 11003.d0*fms(2) - 9402.d0*fms(1) ) + 2107.d0*fms(1)**2
    beta(2) = fms(3)*( 267.d0*fms(3) - 1642.d0*fms(2) + 1602.d0*fms(1) - 494.d0*fms(0) ) &
            + fms(2)*( 2843.d0*fms(2) - 5966.d0*fms(1) + 1922.d0*fms(0) ) &
            + fms(1)*( 3443.d0*fms(1) - 2522.d0*fms(0) ) + 547.d0*fms(0)**2
    beta(3) = fms(2)*( 547.d0*fms(2) - 2522.d0*fms(1) + 1922.d0*fms(0) - 494.d0*fms(-1) ) &
            + fms(1)*( 3443.d0*fms(1) - 5966.d0*fms(0) + 1602.d0*fms(-1) ) &
            + fms(0)*( 2843.d0*fms(0) - 1642.d0*fms(-1) ) + 267.d0*fms(-1)**2
    beta(4) = fms(1)*( 2107.d0*fms(1) - 9402.d0*fms(0) + 7042.d0*fms(-1) - 1854.d0*fms(-2) ) &
            + fms(0)*( 11003.d0*fms(0) - 17246.d0*fms(-1) + 4642.d0*fms(-2) ) &
            + fms(-1)*( 7043.d0*fms(-1) - 3882.d0*fms(-2) ) + 547.d0*fms(-2)**2
    beta(:) = beta(:)/1000.d0
    
    tao7 = abs(beta(1) - beta(4))
    tao = tao7
    aa = 100.d0
    ac = aa*( ( beta(1) + beta(4) - tao + eps )/( tao + eps ) )**2
    tao12 = ( fms(-2) - 6.d0*fms(-1) + 15.d0*fms(0) - 20.d0*fms(1) &
            + 15.d0*fms(2) - 6.d0*fms(3) + fms(4) )**2
    
    alpha(:) = omega(:)*( ac + (tao12/(beta(:)+eps))**2 )
    !alpha(:) = omega(:)*( 1.d0 + (tao12/(beta(:)+eps))**2 )
    totalpha = alpha(1) + alpha(2) + alpha(3)+ alpha(4)
    fomega(:) = alpha(:)/totalpha
    fmx = fu(1)*fomega(1) + fu(2)*fomega(2) + fu(3)*fomega(3)+ fu(4)*fomega(4)
    
    return
    end subroutine wenozn7_m