subroutine wave_speed_calculation_channel(Cell_L,Cell_R,dUldUl,dUrdUr,sm,sl,sr,dSmdUl,dSmdUr,sim)
type(cell_state_t), intent(in) :: Cell_L,Cell_R
real(dp), dimension(Nb_VarC,Nb_VarC) :: dUldUl,dUrdUr
real(dp), intent(out) :: sm,sl,sr
real(dp), dimension(:), intent(out) :: dSmdUl,dSmdUr
type(simulation_t), intent(in) :: sim
real(dp), dimension(Nb_VarP) :: StatL,StatR
real(dp), dimension(Nb_VarC) :: Ul,Ur,SrcL,SrcR
real(dp), dimension(Nb_VarC,Nb_VarC) :: JacobL,JacobR,DerSrcL,DerSrcR
real(dp) :: SS,num,den,dxL,dxR
real(dp), dimension(2) :: Fl,Fr
SrcL=0.0_dp;SrcR=0.0_dp;DerSrcL=0.0_dp;DerSrcR=0.0_dp
StatL=Cell_L%Prim
StatR=Cell_R%Prim
Ul=Cell_L%Cons
Ur=Cell_R%Cons
JacobL=Cell_L%Jcb
JacobR=Cell_R%Jcb
SrcL(2)=Cell_L%SuFrct
SrcR(2)=Cell_R%SuFrct
DerSrcL(:,2)=Cell_L%DerSuFrct(:)
DerSrcR(:,2)=Cell_R%DerSuFrct(:)
dxL=Cell_L%dx
dxR=Cell_R%dx
SS=max(abs(StatL(Pri_u))+StatL(Pri_c),abs(StatR(Pri_u))+StatR(Pri_c))*sim%scheme_diff_fact
sr=SS
sl=-SS
Fl(1)=Ul(2)
Fl(2)=StatL(Pri_ro)*StatL(Pri_u)**2+StatL(Pri_P)
if(R_Correction) Fl(2)=(3.0_dp/2.0_dp)*(Ul(2)**2)/Ul(1)+Ul(4)-Ul(3)
Fr(1)=Ur(2)
Fr(2)=StatR(Pri_ro)*StatR(Pri_u)**2+StatR(Pri_P)
if(R_Correction) Fr(2)=(3.0_dp/2.0_dp)*(Ur(2)**2)/Ur(1)+Ur(4)-Ur(3)
num=Fr(2)-Fl(2)-sr*Ur(2)+sl*Ul(2)-(dxL*SrcL(2)+dxR*SrcR(2))/2.0_dp
den=Fr(1)-Fl(1)-sr*Ur(1)+sl*Ul(1)
sm=num/den
dSmdUl(:)=((-JacobL(:,2)+sl*dUldUl(:,2)-0.5_dp*dxL*DerSrcL(:,2))*den-num*(-JacobL(:,1)+sl*dUldUl(:,1)))/(den**2)
dSmdUr(:)=((JacobR(:,2)-sr*dUrdUr(:,2)-0.5_dp*dxR*DerSrcR(:,2))*den-num*(JacobR(:,1)-sr*dUrdUr(:,1)))/(den**2)
end subroutine wave_speed_calculation_channel