subroutine flux_splitting_LLF(wavolcx, wavolcy)
    use commondata
    implicit none
    
    integer :: j0
    real*8 :: rho, u, v, p, c
    real*8 :: lamda, wavolcx, wavolcy
    
    wavolcx = 0.d0;    wavolcy = 0.d0
    j0 = -4
    if(FLG_prob==3) then    !  Forward-step boundary condition
        call corner_x
        call corner_y
    end if
    do i = -4, im+4
        if(FLG_prob==3) then
            if(i<=ifs+3) then
                j0 = -4
            else
                j0 = jfs-3
            end if
        end if
        do j = j0, jn+4
            u = var(2,i,j)
            p = var(4,i,j);    c = var(5,i,j)
            lamda = abs(u) + c
            wavolcx = max( wavolcx, lamda )
            
            fxp(1,i,j) = 0.5d0*( (u+lamda)*v1(1,i,j) )
            fxp(2,i,j) = 0.5d0*( (u+lamda)*v1(2,i,j) + p )
            fxp(3,i,j) = 0.5d0*( (u+lamda)*v1(3,i,j) )
            fxp(4,i,j) = 0.5d0*( (u+lamda)*v1(4,i,j) + u*p )
            
            fxm(1,i,j) = 0.5d0*( (u-lamda)*v1(1,i,j) )
            fxm(2,i,j) = 0.5d0*( (u-lamda)*v1(2,i,j) + p )
            fxm(3,i,j) = 0.5d0*( (u-lamda)*v1(3,i,j) )
            fxm(4,i,j) = 0.5d0*( (u-lamda)*v1(4,i,j) + u*p )
        end do
    end do
    do i = -4, im+4
        if(FLG_prob==3) then
            if(i<=ifs+3) then
                j0 = -4
            else
                j0 = jfs-3
            end if
        end if
        do j = j0, jn+4
            v = var(3,i,j)
            p = var(4,i,j);    c = var(5,i,j)
            lamda = abs(v) + c
            wavolcy = max( wavolcy, lamda )
            
            fyp(1,i,j) = 0.5d0*( (v+lamda)*v1(1,i,j) )
            fyp(2,i,j) = 0.5d0*( (v+lamda)*v1(2,i,j) )
            fyp(3,i,j) = 0.5d0*( (v+lamda)*v1(3,i,j) + p )
            fyp(4,i,j) = 0.5d0*( (v+lamda)*v1(4,i,j) + v*p )
            
            fym(1,i,j) = 0.5d0*( (v-lamda)*v1(1,i,j) )
            fym(2,i,j) = 0.5d0*( (v-lamda)*v1(2,i,j) )
            fym(3,i,j) = 0.5d0*( (v-lamda)*v1(3,i,j) + p )
            fym(4,i,j) = 0.5d0*( (v-lamda)*v1(4,i,j) + v*p )
        end do
    end do
    
end subroutine flux_splitting_LLF
    
    
subroutine flux_splitting_GLF(wavolcx, wavolcy)
    use commondata
    implicit none
    
    integer :: j0
    real*8 :: rho, u, v, p, c
    real*8 :: lamdax, lamday, wavolcx, wavolcy
    
    wavolcx = 0.d0;    wavolcy = 0.d0
    j0 = -4
    if(FLG_prob==3) then    !  Forward-step boundary condition
        call corner_x
        call corner_y
    end if
    do i = -4, im+4
        if(FLG_prob==3) then
            if(i<=ifs+3) then
                j0 = -4
            else
                j0 = jfs-3
            end if
        end if
        do j = j0, jn+4
            wavolcx = max( wavolcx, abs(var(2,i,j)) + var(5,i,j) )
            wavolcy = max( wavolcy, abs(var(3,i,j)) + var(5,i,j) )
        end do
    end do
    lamdax = wavolcx;    lamday = wavolcy
    do i = -4, im+4
        if(FLG_prob==3) then
            if(i<=ifs+3) then
                j0 = -4
            else
                j0 = jfs-3
            end if
        end if
        do j = j0, jn+4
            u = var(2,i,j)
            p = var(4,i,j);    c = var(5,i,j)
            
            fxp(1,i,j) = 0.5d0*( (u+lamdax)*v1(1,i,j) )
            fxp(2,i,j) = 0.5d0*( (u+lamdax)*v1(2,i,j) + p )
            fxp(3,i,j) = 0.5d0*( (u+lamdax)*v1(3,i,j) )
            fxp(4,i,j) = 0.5d0*( (u+lamdax)*v1(4,i,j) + u*p )
            
            fxm(1,i,j) = 0.5d0*( (u-lamdax)*v1(1,i,j) )
            fxm(2,i,j) = 0.5d0*( (u-lamdax)*v1(2,i,j) + p )
            fxm(3,i,j) = 0.5d0*( (u-lamdax)*v1(3,i,j) )
            fxm(4,i,j) = 0.5d0*( (u-lamdax)*v1(4,i,j) + u*p )
        end do
    end do
    do i = -4, im+4
        if(FLG_prob==3) then
            if(i<=ifs+3) then
                j0 = -4
            else
                j0 = jfs-3
            end if
        end if
        do j = j0, jn+4
            v = var(3,i,j)
            p = var(4,i,j);    c = var(5,i,j)
            
            fyp(1,i,j) = 0.5d0*( (v+lamday)*v1(1,i,j) )
            fyp(2,i,j) = 0.5d0*( (v+lamday)*v1(2,i,j) )
            fyp(3,i,j) = 0.5d0*( (v+lamday)*v1(3,i,j) + p )
            fyp(4,i,j) = 0.5d0*( (v+lamday)*v1(4,i,j) + v*p )
            
            fym(1,i,j) = 0.5d0*( (v-lamday)*v1(1,i,j) )
            fym(2,i,j) = 0.5d0*( (v-lamday)*v1(2,i,j) )
            fym(3,i,j) = 0.5d0*( (v-lamday)*v1(3,i,j) + p )
            fym(4,i,j) = 0.5d0*( (v-lamday)*v1(4,i,j) + v*p )
        end do
    end do
    
end subroutine flux_splitting_GLF
    
    
subroutine flux_splitting_SW(nd, fp, fm, wavolc)
    use commondata
    implicit none
    
    integer :: nd, j0
    real*8 :: fp(kn,-4:im+4,-4:jn+4), fm(kn,-4:im+4,-4:jn+4), wavolc
    real*8 :: alpha, beta, zeta, epsi
    real*8 :: rho, ux, uy, p, c
    real*8 :: c1, cha(5), ch(5), wii
    
    epsi = 1.d0
    wavolc = 0.d0
    j0 = -4
    if(nd==1) then
        alpha = 1.d0
        beta = 0.d0
        if(FLG_prob==3) call corner_x
    else if(nd==2) then
        alpha = 0.d0
        beta = 1.d0
        if(FLG_prob==3) call corner_y
    end if
    zeta = sqrt( alpha**2 + beta**2 )
    do i = -4, im+4
        if(FLG_prob==3) then
            if(i<=ifs+3) then
                j0 = -4
            else
                j0 = jfs-3
            end if
        end if
        do j = j0, jn+4
            rho = var(1,i,j)
            ux = var(2,i,j);    uy = var(3,i,j)
            p  = var(4,i,j);    c  = var(5,i,j)
            c1 = rho/(2.d0*gamma)
            
            cha(1) = alpha*ux + beta*uy
            cha(2) = cha(1)
            cha(3) = cha(1) - zeta*c
            cha(4) = cha(1) + zeta*c
            cha(5) = 0.d0
            wavolc = max( wavolc, abs(cha(3)), abs(cha(4)) )
            
            !ch(:) = 0.5d0*( cha(:) + sqrt( cha(:)**2 + epsi**2 ) )
            do k = 1, 5
                if( k<3 .or. k>4 ) then
                    ch(k) = 0.5d0*( cha(k) + abs(cha(k)) )
                else
                    ch(k) = 0.5d0*( cha(k) + sqrt( cha(k)**2 + epsi**2 ) )
                end if
            end do
            wii = (3.d0-gamma)*( ch(3) + ch(4) )*c**2/(2.d0*(gamma-1.d0))
            fp(1,i,j) = 2.d0*(gamma-1.d0)*ch(1) + ch(3) + ch(4)
            fp(2,i,j) = 2.d0*(gamma-1.d0)*ch(1)*ux &
                + ch(3)*( ux-alpha*c ) + ch(4)*( ux+alpha*c )
            fp(3,i,j) = 2.d0*(gamma-1.d0)*ch(1)*uy &
                + ch(3)*( uy-beta*c ) + ch(4)*( uy+beta*c )
            fp(4,i,j) = (gamma-1.d0)*ch(1)*( ux**2+uy**2 ) + wii &
                + 0.5*ch(3)*( (ux-alpha*c)**2 + (uy-beta*c)**2 ) &
                + 0.5*ch(4)*( (ux+alpha*c)**2 + (uy+beta*c)**2 )
            
            !ch(:) = 0.5d0*( cha(:) - sqrt( cha(:)**2 + epsi**2 ) )
            do k = 1, 5
                if( k<3 .or. k>4 ) then
                    ch(k) = 0.5d0*( cha(k) - abs(cha(k)) )
                else
                    ch(k) = 0.5d0*( cha(k) - sqrt( cha(k)**2 + epsi**2 ) )
                end if
            end do
            wii = (3.d0-gamma)*( ch(3) + ch(4) )*c**2/(2.d0*(gamma-1.d0))
            fm(1,i,j) = 2.d0*(gamma-1.d0)*ch(1) + ch(3) + ch(4)
            fm(2,i,j) = 2.d0*(gamma-1.d0)*ch(1)*ux &
                + ch(3)*( ux-alpha*c ) + ch(4)*( ux+alpha*c )
            fm(3,i,j) = 2.d0*(gamma-1.d0)*ch(1)*uy &
                + ch(3)*( uy-beta*c ) + ch(4)*( uy+beta*c )
            fm(4,i,j) = (gamma-1.d0)*ch(1)*( ux**2+uy**2 ) + wii &
                + 0.5*ch(3)*( (ux-alpha*c)**2 + (uy-beta*c)**2 ) &
                + 0.5*ch(4)*( (ux+alpha*c)**2 + (uy+beta*c)**2 )
            
            fp(:,i,j) = fp(:,i,j)*c1
            fm(:,i,j) = fm(:,i,j)*c1
        end do
    end do
    
end subroutine flux_splitting_SW
    
    
subroutine flux_splitting_VL(nd, fp, fm, wavolc)
    use commondata
    implicit none
    
    integer :: nd, j0
    real*8 :: fp(kn,-4:im+4,-4:jn+4), fm(kn,-4:im+4,-4:jn+4), wavolc
    real*8 :: alpha, beta, zeta
    real*8 :: rho, ux, uy, p, c, uv, xm
    real*8 :: gamma1, f1p, f1m, phi
    
    wavolc = 0.d0
    j0 = -4
    if(nd==1) then
        alpha = 1.d0
        beta = 0.d0
        if(FLG_prob==3) call corner_x
    else if(nd==2) then
        alpha = 0.d0
        beta = 1.d0
        if(FLG_prob==3) call corner_y
    end if
    zeta = sqrt( alpha**2 + beta**2 )
    do i = -4, im+4
        if(FLG_prob==3) then
            if(i<=ifs+3) then
                j0 = -4
            else
                j0 = jfs-3
            end if
        end if
        do j = j0, jn+4
            rho = var(1,i,j)
            ux = var(2,i,j);    uy = var(3,i,j)
            p  = var(4,i,j);    c  = var(5,i,j)
            uv = alpha*ux + beta*uy
            wavolc = max( wavolc, abs(uv-c), abs(uv+c) )
            xm = uv/c
            if(xm>=1.d0) then
                fp(1,i,j) = rho*uv
                fp(2,i,j) = rho*uv*ux + alpha*p
                fp(3,i,j) = rho*uv*uy + beta*p
                fp(4,i,j) = uv*( p/(gamma-1.d0) + 0.5d0*rho*( ux**2+uy**2 ) + p )
                fm(:,i,j) = 0.d0
            else if(xm<=-1.d0) then
                fm(1,i,j) = rho*uv
                fm(2,i,j) = rho*uv*ux + alpha*p
                fm(3,i,j) = rho*uv*uy + beta*p
                fm(4,i,j) = uv*( p/(gamma-1.d0) + 0.5d0*rho*( ux**2+uy**2 ) + p )
                fp(:,i,j) = 0.d0
            else
                gamma1 = gamma - 1.d0
                f1p =  zeta*rho*c*( xm+1.d0 )**2/4.d0
                f1m = -zeta*rho*c*( xm-1.d0 )**2/4.d0
                phi = gamma1/2.d0*( ux**2+uy**2 )
                fp(1,i,j) = f1p
                fp(2,i,j) = f1p*( alpha*(-uv+2.d0*c)/gamma + ux )
                fp(3,i,j) = f1p*(  beta*(-uv+2.d0*c)/gamma + uy )
                fp(4,i,j) = f1p*( (-gamma1*uv**2 + 2.d0*gamma1*uv*c + 2.d0*c**2)/(gamma**2-1.d0) + phi/gamma1 )
                fm(1,i,j) = f1m
                fm(2,i,j) = f1m*( alpha*(-uv-2.d0*c)/gamma + ux )
                fm(3,i,j) = f1m*(  beta*(-uv-2.d0*c)/gamma + uy )
                fm(4,i,j) = f1m*( (-gamma1*uv**2 - 2.d0*gamma1*uv*c + 2.d0*c**2)/(gamma**2-1.d0) + phi/gamma1 )                
            end if
        end do
    end do
    
end subroutine flux_splitting_VL