subroutine incoming_branch_cold_circulator(me,mdot,dPres,VarIn,deriv_wr_both_br,ssIn,FluxIn)
type(circulator_t), intent(inout) :: me
real(dp), intent(in) :: mdot,dPres
real(dp), intent(out) :: VarIn,ssIn
type(derivatives_involving_both_branches_t), intent(inout) :: deriv_wr_both_br
type(flux_and_derivatives_for_cold_circulator_t), intent(out) :: FluxIn
real(dp) :: Rol,uL,pL,rL,cL,ML,Zl,Area,Sl,g
real(dp) :: Rostar,pstar,Ustar,R_Correction_Star,eintStar,Tstar,Cstar
real(dp) :: BB,HHH,EtotStar
real(dp), dimension(Nb_VarC) :: U,Cons_star,DerUstardUl,DerRostardUl,DerRcorrstardUl,DerPstardUl,DerBBdUl
real(dp), dimension(Nb_VarC) :: DerUstardUoth,DerRostardUoth,DerRcorrstardUoth,DerPstardUoth
real(dp), dimension(Nb_VarC) :: dMdotdUl,dMdotdUoth
real(dp), dimension(Nb_VarC,Nb_VarC) :: Jacob_star,DerConstardUl,DerConstardUoth
real(dp) :: kk, KKK, dedT, dPdT, HH, dPde
real(dp), dimension(Nb_VarP) :: PPP
real(dp), dimension(Nb_VarP) :: PrimStar
real(dp) :: dTdp_Ro,dTdRo_p,dRstardRo,dRstardT,Rstar
Rol=me%br(1)%p%Prim(Pri_ro)
uL=me%br(1)%p%Prim(Pri_u)
pL=me%br(1)%p%Prim(Pri_P)
rL=me%br(1)%p%Prim(Pri_R)
cL=me%br(1)%p%Prim(Pri_c)
ML=uL/cL
Zl=Rol*cL
Area=me%br(1)%p%Area
U(:)=me%br(1)%p%Cons(:)
! Approximate RH relation
Sl=-(abs(uL)+cL)
Ustar=Sl*mdot/(mdot-Rol*(uL-Sl)*Area)
Rostar=Rol*(uL-Sl)/(Ustar-Sl)
pstar=pL+Rol*(uL-Sl)*(uL-Ustar)
call state_roP(Rostar, pstar, R_Correction_Star, eintStar, Tstar, Cstar)
if(R_Correction) R_Correction_Star=rL*(uL-Sl)/(Ustar-Sl)
eintStar=(R_Correction_Star-pstar)/Rostar
EtotStar=eintStar+0.5_dp*Ustar**2
FluxIn%flx(Con_Mas)=Rostar*Ustar
FluxIn%flx(Con_Qdm)=Rostar*Ustar*Ustar+pstar
FluxIn%flx(Con_Ene)=(Rostar*EtotStar+pstar)*Ustar
FluxIn%flx(Con_R)=0.0_dp
if(R_Correction) FluxIn%flx(Con_R)=R_Correction_Star*Ustar
FluxIn%vit=Ustar
!------------------------------------------------------------------!
! Derivatives !
dMdotdUl(:)=deriv_wr_both_br%DerMdotDerUin(:)
DerUstardUl(1)=(Sl*dMdotdUl(1)*(mdot-(U(2)-U(1)*Sl)*Area)-Sl*mdot*(dMdotdUl(1)+Sl*Area))/((mdot-(U(2)-U(1)*Sl)*Area)**2)
DerUstardUl(2)=(Sl*dMdotdUl(2)*(mdot-(U(2)-U(1)*Sl)*Area)-Sl*mdot*(dMdotdUl(2)-Area))/((mdot-(U(2)-U(1)*Sl)*Area)**2)
DerUstardUl(3)=(Sl*dMdotdUl(3)*(mdot-(U(2)-U(1)*Sl)*Area)-Sl*mdot*dMdotdUl(3))/((mdot-(U(2)-U(1)*Sl)*Area)**2)
DerUstardUl(4)=(Sl*dMdotdUl(4)*(mdot-(U(2)-U(1)*Sl)*Area)-Sl*mdot*dMdotdUl(4))/((mdot-(U(2)-U(1)*Sl)*Area)**2)
DerRostardUl(1)=(-Sl*(Ustar-Sl)-(U(2)-U(1)*Sl)*DerUstardUl(1))/((Ustar-Sl)**2)
DerRostardUl(2)=((Ustar-Sl)-(U(2)-U(1)*Sl)*DerUstardUl(2))/((Ustar-Sl)**2)
DerRostardUl(3)=-(U(2)-U(1)*Sl)*DerUstardUl(3)/((Ustar-Sl)**2)
DerRostardUl(4)=-(U(2)-U(1)*Sl)*DerUstardUl(4)/((Ustar-Sl)**2)
if(R_Correction) then
Cons_star(Con_Mas)=Rostar
Cons_star(Con_Qdm)=Rostar*Ustar
Cons_star(Con_Ene)=Rostar*EtotStar
Cons_star(Con_R)=R_Correction_Star
call jacobian_star(Cons_star,Jacob_star)
BB=U(4)-U(3)+0.5_dp*(U(2)**2)/U(1)
DerBBdUl(1)=-0.5_dp*(U(2)**2)/(U(1)**2)
DerBBdUl(2)=U(2)/U(1)
DerBBdUl(3)=-1.0_dp
DerBBdUl(4)=1.0_dp
DerPstardUl(1)=DerBBdUl(1)-Sl*(U(2)/U(1)-Ustar)+(U(2)-U(1)*Sl)*(-U(2)/(U(1)**2)-DerUstardUl(1))
DerPstardUl(2)=DerBBdUl(2)+(U(2)/U(1)-Ustar)+(U(2)-U(1)*Sl)*(1.0_dp/U(1)-DerUstardUl(2))
DerPstardUl(3)=DerBBdUl(3)-(U(2)-U(1)*Sl)*DerUstardUl(3)
DerPstardUl(4)=DerBBdUl(4)-(U(2)-U(1)*Sl)*DerUstardUl(4)
DerRcorrstardUl(1)=-U(4)*(U(2)*(Ustar-Sl)/(U(1)**2)+(U(2)/U(1)-Sl)*DerUstardUl(1))/((Ustar-Sl)**2)
DerRcorrstardUl(2)=U(4)*((Ustar-Sl)/U(1)-(U(2)/U(1)-Sl)*DerUstardUl(2))/((Ustar-Sl)**2)
DerRcorrstardUl(3)=-U(4)*(U(2)/U(1)-Sl)*DerUstardUl(3)/((Ustar-Sl)**2)
DerRcorrstardUl(4)=(U(2)/U(1)-Sl)*(Ustar-Sl-U(4)*DerUstardUl(4))/((Ustar-Sl)**2)
DerConstardUl(:,1)=DerRostardUl(:)
DerConstardUl(:,2)=Rostar*DerUstardUl(:)+DerRostardUl(:)*Ustar
DerConstardUl(:,3)=DerRcorrstardUl(:)+0.5_dp*(DerRostardUl(:)*Ustar*Ustar+2.0_dp*Rostar*Ustar*DerUstardUl(:))-DerPstardUl(:)
DerConstardUl(:,4)=DerRcorrstardUl(:)
FluxIn%DerFlxdBr_a(:,:)=matmul(DerConstardUl(:,:),Jacob_star(:,:))
FluxIn%DerVitdBr_a(:)=DerUstardUl(:)
else
PrimStar(Pri_ro)=Rostar
PrimStar(Pri_u)=Ustar
PrimStar(Pri_p)=pstar
PrimStar(Pri_e)=eintStar
PrimStar(Pri_T)=Tstar
PrimStar(Pri_c)=Cstar
call jacobian_star_Req_rem(PrimStar,Jacob_star)
PPP(:)=me%br(1)%p%Prim(:)
HH=PPP(Pri_e)+0.5_dp*PPP(Pri_u)**2+PPP(Pri_p)/PPP(Pri_ro)
call jacobian_roT(PPP(Pri_ro), PPP(Pri_T), dedT, dPdT)
dPde=dPdT/dedT
kk=dPde/PPP(Pri_ro)
KKK=PPP(Pri_c)**2+kk*(PPP(Pri_u)**2-HH)
BB=pL
DerBBdUl(1)=KKK
DerBBdUl(2)=-kk*PPP(Pri_u)
DerBBdUl(3)=kk
DerBBdUl(4)=0.0_dp
DerPstardUl(1)=DerBBdUl(1)-Sl*(U(2)/U(1)-Ustar)+(U(2)-U(1)*Sl)*(-U(2)/(U(1)**2)-DerUstardUl(1))
DerPstardUl(2)=DerBBdUl(2)+(U(2)/U(1)-Ustar)+(U(2)-U(1)*Sl)*(1.0_dp/U(1)-DerUstardUl(2))
DerPstardUl(3)=DerBBdUl(3)-(U(2)-U(1)*Sl)*DerUstardUl(3)
DerPstardUl(4)=DerBBdUl(4)-(U(2)-U(1)*Sl)*DerUstardUl(4)
call jacobian_roT(Rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
DerRcorrstardUl(1)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUl(1)+dRstardT*dTdp_Ro*DerPstardUl(1)
DerRcorrstardUl(2)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUl(2)+dRstardT*dTdp_Ro*DerPstardUl(2)
DerRcorrstardUl(3)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUl(3)+dRstardT*dTdp_Ro*DerPstardUl(3)
DerRcorrstardUl(4)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUl(4)+dRstardT*dTdp_Ro*DerPstardUl(4)
DerConstardUl(:,1)=DerRostardUl(:)
DerConstardUl(:,2)=Rostar*DerUstardUl(:)+DerRostardUl(:)*Ustar
DerConstardUl(:,3)=DerRcorrstardUl(:)+0.5_dp*(DerRostardUl(:)*Ustar*Ustar+2.0_dp*Rostar*Ustar*DerUstardUl(:))-DerPstardUl(:)
DerConstardUl(:,4)=0.0_dp
FluxIn%DerFlxdBr_a(:,:)=matmul(DerConstardUl(:,:),Jacob_star(:,:))
FluxIn%DerFlxdBr_a(:,4)=0.0_dp
FluxIn%DerVitdBr_a(:)=DerUstardUl(:)
endif
!----------------------------------------------------------------------------------------------------------
! Contributions of the CC outlet
dMdotdUoth(:)=deriv_wr_both_br%DerMdotDerUout(:)
DerUstardUoth(:)=(Sl*dMdotdUoth(:)*(mdot-(U(2)-U(1)*Sl)*Area)-Sl*mdot*dMdotdUoth(:))/((mdot-(U(2)-U(1)*Sl)*Area)**2)
DerRostardUoth(:)=-(U(2)-U(1)*Sl)*DerUstardUoth(:)/((Ustar-Sl)**2)
DerPstardUoth(:)=-(U(2)-U(1)*Sl)*DerUstardUoth(:)
if(R_Correction) then
DerRcorrstardUoth(:)=-U(4)*(U(2)/U(1)-Sl)*DerUstardUoth(:)/((Ustar-Sl)**2)
else
call jacobian_roT(Rostar, Tstar, dedT, dPdT, dTdp_Ro, dTdRo_p, dRstardRo, dRstardT, Rstar)
DerRcorrstardUoth(:)=(dRstardRo+dRstardT*dTdRo_p)*DerRostardUoth(:)+dRstardT*dTdp_Ro*DerPstardUoth(:) ! DerRcorrstardUoth(4) already equals to 0
endif
DerConstardUoth(:,1)=DerRostardUoth(:)
DerConstardUoth(:,2)=Rostar*DerUstardUoth(:)+DerRostardUoth(:)*Ustar
DerConstardUoth(:,3)=DerRcorrstardUoth(:)+0.5_dp*(DerRostardUoth(:)*Ustar*Ustar+2.0_dp*Rostar*Ustar*DerUstardUoth(:))-&
DerPstardUoth(:)
DerConstardUoth(:,4)=DerRcorrstardUoth(:)
FluxIn%DerFlxdBr_b(:,:)=matmul(DerConstardUoth(:,:),Jacob_star(:,:))
if(.not. R_Correction) FluxIn%DerFlxdBr_b(:,4)=0.0_dp
FluxIn%DerVitdBr_b(:)=DerUstardUoth(:)
!------------------------------------------------------------------------------------------------------
g=9.81_dp
HHH=imposed_compression_energy(me,Rostar,g,dPres,me%dp0)
if(me%SubType=="pump") then
VarIn=R_Correction_Star/Rostar+0.5_dp*Ustar**2+HHH
call HHH_derivatives(me,rostar,g,DerRostardUl,deriv_wr_both_br)
deriv_wr_both_br%DerVarIndUin(:)=(DerRcorrstardUl(:)*Rostar-R_Correction_Star*DerRostardUl(:))/&
(Rostar**2)+Ustar*DerUstardUl(:)+deriv_wr_both_br%DerHHHdUin(:)
deriv_wr_both_br%DerVarIndUout(:)=(DerRcorrstardUoth(:)*Rostar-R_Correction_Star*DerRostardUoth(:))/&
(Rostar**2)+Ustar*DerUstardUoth(:)+deriv_wr_both_br%DerHHHdUout(:)
else if(me%SubType=="compressor") then
VarIn=R_Correction_Star/Rostar+0.5_dp*Ustar**2+g*HHH
call HHH_derivatives(me,rostar,g,DerRostardUl,deriv_wr_both_br)
deriv_wr_both_br%DerVarIndUin(:)=(DerRcorrstardUl(:)*Rostar-R_Correction_Star*DerRostardUl(:))/&
(Rostar**2)+Ustar*DerUstardUl(:)+g*deriv_wr_both_br%DerHHHdUin(:)
deriv_wr_both_br%DerVarIndUout(:)=(DerRcorrstardUoth(:)*Rostar-R_Correction_Star*DerRostardUoth(:))/&
(Rostar**2)+Ustar*DerUstardUoth(:)+g*deriv_wr_both_br%DerHHHdUout(:)
endif
ssIn=abs(uL)+cL
end subroutine incoming_branch_cold_circulator