subroutine flux_internal_calculation_channel(i,HepropPP,Cell_L,Cell_R,dtMaxLoc,flux,sim)
type(cell_state_t), intent(in) :: Cell_L,Cell_R
real(dp), intent(out) :: dtMaxLoc
type(flux_He_channel_t), intent(inout) :: flux
integer, intent(in) :: i
type(He_prop_t), intent(in) :: HePropPP
type(simulation_t), intent(in) :: sim
real(dp), dimension(Nb_VarC,Nb_VarC) :: dUldUl,dUrdUr
real(dp), dimension(Nb_VarP) :: StatLstar,StatRstar
real(dp), dimension(Nb_VarC) :: Fl,Fr,dSmdUl,dSmdUr,dpstardUl,dpstardUr,ConsLstar,ConsRstar
real(dp), dimension(Nb_VarC) :: Fstar,LastTermL,LastTermR
real(dp), dimension(Nb_VarC,Nb_VarC) :: derUstarL_dUL,derUstarL_dUR,derUstarR_dUL,derUstarR_dUR
real(dp), dimension(Nb_VarC,Nb_VarC) :: dFstardUl,dFstardUr,LastMatrix,AddThermL,AddThermR
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) :: Etstar,pstar,dxL,dxR
real(dp) :: sm,sl,sr
integer :: n,imPP
StatLstar=0.0_dp;StatRstar=0.0_dp;Fl=0.0_dp;Fr=0.0_dp
imPP=HepropPP%NbCells
dUldUl(:,:)=id_4x4(:,:)
dUrdUr(:,:)=id_4x4(:,:)
SrcL=0.0_dp;SrcR=0.0_dp;DerSrcL=0.0_dp;DerSrcR=0.0_dp
call wave_speed_calculation_channel(Cell_L,Cell_R,dUldUl,dUrdUr,sm,sl,sr,dSmdUl(:),dSmdUr(:),sim)
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(:)
if(R_Correction) then
SrcL(4)=Cell_L%SrFrct
SrcR(4)=Cell_R%SrFrct
DerSrcL(:,4)=Cell_L%DerSrFrct(:)
DerSrcR(:,4)=Cell_R%DerSrFrct(:)
endif
dxL=Cell_L%dx
dxR=Cell_R%dx
if(max(abs(sl),abs(sr))<=1.0e-8_dp) then
call set_error('flux denominator zero - CFL criterion for non-uniform mesh channels')
return
endif
if(i==0) then
dtMaxLoc=Cell_R%dx/max(abs(sl),abs(sr))
else if(i==imPP) then
dtMaxLoc=Cell_L%dx/max(abs(sl),abs(sr))
else
dtMaxLoc=min(Cell_L%dx/max(abs(sl),abs(sr)),Cell_R%dx/max(abs(sl),abs(sr))) ! Facteur 1/2 ?
endif
Fl(Con_Mas)=StatL(Pri_ro)*StatL(Pri_u)
Fl(Con_Qdm)=StatL(Pri_ro)*StatL(Pri_u)*StatL(Pri_u)+StatL(Pri_p)
Fl(Con_Ene)=(StatL(Pri_ro)*(StatL(Pri_e)+0.5_dp*StatL(Pri_u)*StatL(Pri_u))+StatL(Pri_p))*StatL(Pri_u)
Fl(Con_R)=0.0_dp
if(R_Correction) Fl(Con_R)=StatL(Pri_R)*StatL(Pri_u)
Fr(Con_Mas)=StatR(Pri_ro)*StatR(Pri_u)
Fr(Con_Qdm)=StatR(Pri_ro)*StatR(Pri_u)*StatR(Pri_u)+StatR(Pri_p)
Fr(Con_Ene)=(StatR(Pri_ro)*(StatR(Pri_e)+0.5_dp*StatR(Pri_u)*StatR(Pri_u))+StatR(Pri_p))*StatR(Pri_u)
Fr(Con_R)=0.0_dp
if(R_Correction) Fr(Con_R)=StatR(Pri_R)*StatR(Pri_u)
! Wave sampling
pstar=0.5_dp*(Fr(2)-sr*Ur(2)-sm*(Fr(1)-sr*Ur(1))+Fl(2)-sl*Ul(2)-sm*(Fl(1)-sl*Ul(1)))
if(sm>=0.0_dp) then
! Left star state
StatLstar(Pri_ro)=StatL(Pri_ro)*(StatL(Pri_u)-sl)/(sm-sl)
Etstar=(StatL(Pri_e)+0.5_dp*StatL(Pri_u)*StatL(Pri_u))+(StatL(Pri_p)*&
StatL(Pri_u)-pstar*sm)/(StatL(Pri_ro)*(StatL(Pri_u)-sl))
StatLstar(Pri_R)=(StatL(Pri_R)*(StatL(Pri_u)-sl)+dxL*SrcL(4)/2.0_dp)/(sm-sl)
flux%Cons(Con_Mas,i)=Fl(Con_Mas)+sl*(StatLstar(Pri_ro)-StatL(Pri_ro))
flux%Cons(Con_Qdm,i)=Fl(Con_Qdm)+sl*(StatLstar(Pri_ro)*sm-StatL(Pri_ro)*StatL(Pri_u))+dxL*SrcL(2)/2.0_dp
flux%Cons(Con_Ene,i)=Fl(Con_Ene)+sl*(StatLstar(Pri_ro)*Etstar-StatL(Pri_ro)*(StatL(Pri_e)+&
0.5_dp*StatL(Pri_u)*StatL(Pri_u)))
flux%Cons(Con_R,i)=0.0_dp
if(R_Correction) flux%Cons(Con_R,i)=Fl(Con_R)+sl*(StatLstar(Pri_R)-StatL(Pri_R))+dxL*SrcL(4)/2.0_dp
flux%VitTNC(i)=sm
else
! Right star state
StatRstar(Pri_ro)=StatR(Pri_ro)*(StatR(Pri_u)-sr)/(sm-sr)
Etstar=(StatR(Pri_e)+0.5_dp*StatR(Pri_u)*StatR(Pri_u))+(StatR(Pri_p)*&
StatR(Pri_u)-pstar*sm)/(StatR(Pri_ro)*(StatR(Pri_u)-sr))
StatRstar(Pri_R)=(StatR(Pri_R)*(StatR(Pri_u)-sr)-dxR*SrcR(4)/2.0_dp)/(sm-sr)
flux%Cons(Con_Mas,i)=Fr(Con_Mas)+sr*(StatRstar(Pri_ro)-StatR(Pri_ro))
flux%Cons(Con_Qdm,i)=Fr(Con_Qdm)+sr*(StatRstar(Pri_ro)*sm-StatR(Pri_ro)*StatR(Pri_u))-dxR*SrcR(2)/2.0_dp
flux%Cons(Con_Ene,i)=Fr(Con_Ene)+sr*(StatRstar(Pri_ro)*Etstar-StatR(Pri_ro)*(StatR(Pri_e)+&
0.5_dp*StatR(Pri_u)*StatR(Pri_u)))
flux%Cons(Con_R,i)=0.0_dp
if(R_Correction) flux%Cons(Con_R,i)=Fr(Con_R)+sr*(StatRstar(Pri_R)-StatR(Pri_R))-dxR*SrcR(4)/2.0_dp
flux%VitTNC(i)=sm
endif
! Flux derivative calculation
Fstar(1)=0.0_dp
Fstar(2)=pstar
Fstar(3)=pstar*sm
Fstar(4)=0.0_dp
LastTermL(:)=Fl(:)-sl*Ul(:)-Fstar(:)+dxL*SrcL(:)/2.0_dp
LastTermR(:)=Fr(:)-sr*Ur(:)-Fstar(:)-dxR*SrcR(:)/2.0_dp
if(.not. R_Correction) LastTermL(4)=0.0_dp; LastTermR(4)=0.0_dp
ConsLstar(:)=LastTermL(:)/(sm-sl)
ConsRstar(:)=LastTermR(:)/(sm-sr)
do n=1,Nb_VarC
AddThermL(:,n)=(ConsRstar(n)-ConsLstar(n))*dSmdUl(:)
AddThermR(:,n)=(ConsRstar(n)-ConsLstar(n))*dSmdUr(:)
enddo
! DerFlux_dUL calculation
dpstardUl(:)=0.5_dp*(JacobL(:,2)-sl*dUldUl(:,2)-sm*(JacobL(:,1)-sl*dUldUl(:,1))-&
(Fl(1)-sl*Ul(1))*dSmdUl(:)-(Fr(1)-sr*Ur(1))*dSmdUl(:))
dFstardUl(:,:)=0.0_dp
dFstardUl(:,2)=dpstardUl(:)
dFstardUl(:,3)=sm*dpstardUl(:)+pstar*dSmdUl(:)
do n=1,4
LastMatrix(:,n)=LastTermL(n)*dSmdUl(:)
enddo
derUstarL_dUL(:,:)=((JacobL(:,:)-sl*id_4x4(:,:)-dFstardUl(:,:)+dxL*DerSrcL(:,:)/2.0_dp)*(sm-sl)-LastMatrix(:,:))/((sm-sl)**2)
do n=1,4
LastMatrix(:,n)=LastTermR(n)*dSmdUl(:)
enddo
derUstarR_dUL(:,:)=((-dFstardUl(:,:))*(sm-sr)-LastMatrix(:,:))/((sm-sr)**2)
if(i/=0) THEN
flux%Cons_DerdUL(:,:,i)=0.5_dp*JacobL(:,:)-&
sign(1.0_dp,sl)*(sl/2.0_dp)*(derUstarL_dUL(:,:)-id_4x4(:,:))-&
sign(1.0_dp,sm)*(sm/2.0_dp)*(derUstarR_dUL(:,:)-derUstarL_dUL(:,:))+&
sign(1.0_dp,sr)*(sr/2.0_dp)*derUstarR_dUL(:,:)-&
0.5_dp*sign(1.0_dp,sm)*AddThermL(:,:)
flux%VitTNC_DerdUL(:,i)=dSmdUl(:)
endif
! DerFlux_dUR calculation
dpstardUr(:)=0.5_dp*(JacobR(:,2)-sr*dUrdUr(:,2)-sm*(JacobR(:,1)-sr*dUrdUr(:,1))-&
(Fr(1)-sr*Ur(1))*dSmdUr(:)-(Fl(1)-sl*Ul(1))*dSmdUr(:))
dFstardUr(:,:)=0.0_dp
dFstardUr(:,2)=dpstardUr(:)
dFstardUr(:,3)=sm*dpstardUr(:)+pstar*dSmdUr(:)
do n=1,4
LastMatrix(:,n)=LastTermL(n)*dSmdUr(:)
enddo
derUstarL_dUR(:,:)=((-dFstardUr(:,:))*(sm-sl)-LastMatrix(:,:))/((sm-sl)**2)
do n=1,4
LastMatrix(:,n)=LastTermR(n)*dSmdUr(:)
enddo
derUstarR_dUR(:,:)=((JacobR(:,:)-sr*id_4x4(:,:)-dFstardUr(:,:)-dxR*DerSrcR(:,:)/2.0_dp)*(sm-sr)-LastMatrix(:,:))/((sm-sr)**2)
if(i/=imPP) then
flux%Cons_DerdUR(:,:,i)=0.5_dp*JacobR(:,:)-&
sign(1.0_dp,sl)*(sl/2.0_dp)*derUstarL_dUR(:,:)-&
sign(1.0_dp,sm)*(sm/2.0_dp)*(derUstarR_dUR(:,:)-derUstarL_dUR(:,:))-&
sign(1.0_dp,sr)*(sr/2.0_dp)*(id_4x4(:,:)-derUstarR_dUR(:,:))-&
0.5_dp*sign(1.0_dp,sm)*AddThermR(:,:)
flux%VitTNC_DerdUR(:,i)=dSmdUr(:)
endif
end subroutine flux_internal_calculation_channel