solve_junction Subroutine

subroutine solve_junction(dyn,FlJ)
    type(junction_dynamic_parameters_t), intent(inout) :: dyn
    type(flux_and_derivatives_t), intent(inout) :: FlJ(:)

    type(conservative_variable_derivatives_t),dimension(dyn%NbTotBr,dyn%NbTotBr) :: DerConStar
    type(primitive_variables_t),              dimension(dyn%NbTotBr) :: PrimStar
    type(conservative_variables_t),           dimension(dyn%NbTotBr) :: ConStar
    type(der_cons_vars_t), dimension(dyn%NbTotBr,dyn%NbTotBr+dyn%NbOut) :: DerFct,DerXStar
    type(der_cons_vars_t), dimension(dyn%NbTotBr,dyn%NbOut) :: DerCoef,DerAlphaGlob
    type(der_cons_vars_t), dimension(dyn%NbTotBr,dyn%NbTotBr) :: DerPStarGlob,DerVelStarGlob
    type(der_cons_vars_t), dimension(dyn%NbTotBr,dyn%NbTotBr) :: DerRPhysStarGlob,DerRCorStarGlob
    type(der_cons_vars_t), dimension(dyn%NbTotBr,dyn%NbTotBr) :: DerRhoStarGlob
    real(dp), dimension(dyn%NbTotBr) :: p_star,rhoStar,vitStar,ETotStar,rCorrStar,EnthalpyStar,eStar
    real(dp), dimension(dyn%NbTotBr) :: rhoStar_Int,QDMStar_Int,EnergyStar_Int,QDMStar,EnergyStar
    real(dp), dimension(dyn%NbTotBr) :: DerVitStar_dU1,DerVitStar_dU2,DerVitStar_dU3,DerVitStar_dU4
    real(dp), dimension(dyn%NbTotBr) :: DerRhoStar_dU1,DerRhoStar_dU2,DerRhoStar_dU3,DerRhoStar_dU4
    real(dp), dimension(dyn%NbTotBr) :: DerETotStar_Int_dU1,DerETotStar_Int_dU2,DerETotStar_Int_dU3
    real(dp), dimension(dyn%NbTotBr) :: DerHStar_dU1,DerHStar_dU2,DerHStar_dU3,DerHStar_dU4
    real(dp), dimension(dyn%NbTotBr) :: DerETotStar_dU1,DerETotStar_dU2,DerETotStar_dU3
    real(dp), dimension(dyn%NbTotBr) :: DerRCorrStar_dU1,DerRCorrStar_dU2,DerRCorrStar_dU3
    real(dp), dimension(dyn%NbTotBr) :: DerEntMixStar_dU1,DerEntMixStar_dU2,DerEntMixStar_dU3
    real(dp), dimension(dyn%NbTotBr) :: DerVitStar_dX,DerRhoStar_dX,DerETotStar_Int_dX
    real(dp), dimension(dyn%NbTotBr) :: DerRCorrStar_dX,DerEntMixStar_dX,DerHStar_dAlpha
    real(dp), dimension(dyn%NbTotBr) :: DerRhoETotStar_dAlpha,DerRhoETotStar_dX,DerHStar_dX
    real(dp), dimension(dyn%NbTotBr) :: DerETotStar_dX,DerETotStar_Int_dU4,DerEntMixStar_dU4
    real(dp), dimension(dyn%NbTotBr) :: DerRCorrStar_dU4,DerETotStar_dU4
    real(dp), dimension(Nb_VarC) :: DerBBdUk,DerDDdUk,DerRhoETotStar_dUk,Der_qj_dUk,Der_qRef_dUk
    real(dp), dimension(Nb_VarC,Nb_VarC) :: dFdFJ,dUdUjj
    real(dp), dimension(Nb_VarC,Nb_VarC,dyn%NbTotBr) :: Jacob_star
    real(dp), dimension(dyn%NbTotBr+dyn%NbOut) :: x, fVec
    real(dp), dimension(dyn%NbTotBr+dyn%NbOut,dyn%NbTotBr+dyn%NbOut) :: fJac, invJac
    real(dp), dimension(dyn%NbOut+1) :: DerCoef_dX
    real(dp), dimension(dyn%NbOut) :: Der_dFdRo_dP,ETotStar_Int
    real(dp), dimension(dyn%NbOut) :: alpha,dFdRo,CoefPLoss,Der_dFdRo_dRhoInt
    real(dp), dimension(dyn%NbTotBr) :: T_temp
    real(dp) :: BB,DD,dTdp_Ro,dTdRo_p,dRStar_dRo,dRStar_dT,RStar
    real(dp) :: denominator,EnthalpyMixStar,EnthalpyMixStarNum,qj,ps_ij
    real(dp) :: Der_qjdX,Der_qRef_dX,Der_qjdAlpha,DerCoef_dAlpha
    real(dp) :: kk, KKK, dedT, dPdT, HH, dPde
    integer :: info, NbOut, j, jj, NbBr, NbIn, i(dyn%NbIn), o(dyn%NbOut)
    logical, dimension(dyn%NbTotBr) :: om, im  
    real(dp), parameter :: gg = 9.81_dp

    ! Preparing short name for indexing, masking and sizes
    NbBr  = dyn%NbTotBr
    NbIn  = dyn%NbIn
    NbOut = dyn%NbOut
    om = .false.; im(1:NbOut)      = .true.
    im = .false.; im(NbOut+1:NbBr) = .true.
    i = [(j, j = NbOut+1, NbBr)]
    o = [(j, j = 1, dyn%NbOut)]

    ! Calling fSolve to provide solution of `p_star` and `alpha`
    x(1:NbBr) = dyn%br%p0
    x(NbBr+1:NbBr+NbOut) = 1.0e-1_dp
    call fSolve(dyn,NbBr+NbOut,x,fVec,1.0e-8_dp,info)
    if (info /= 1) then
        x(1:NbBr) = dyn%p_star_prev
        x(NbBr+1:NbBr+NbOut) = dyn%alpha_prev(1:NbOut)
        call fSolve(dyn,NbBr+NbOut,x,fVec,1.0e-8_dp,info)
        if (info /= 1) then
          call set_error('fSolve did not converge in solve_junction')
          return
        endif          
    end if
    p_star = x(1:NbBr)
    alpha  = x(NbBr+1:NbBr+NbOut)
    dyn%p_star_prev = p_star
    dyn%alpha_prev(1:NbOut) = alpha

    ! Incoming quantities as a function of p_star(j) 

    ! Approximate jump relations through the rarefaction wave
    vitStar(i)=dyn%br(i)%vit0-(p_star(i)-dyn%br(i)%p0-dyn%br(i)%rho0*(dyn%br(i)%q0 &
        -dyn%br(i)%qBar) + dyn%br(i)%dx0*dyn%br(i)%Frc/2.0_dp )/(dyn%br(i)%rho0 &
        *(dyn%br(i)%vit0-dyn%br(i)%SpeedS0))
    rhoStar(i)=dyn%br(i)%rho0*(dyn%br(i)%vit0-dyn%br(i)%SpeedS0)/(vitStar(i)-dyn%br(i)%SpeedS0)

    if(R_Correction) then
        rCorrStar(i)=(dyn%br(i)%rCorr0*(dyn%br(i)%vit0-dyn%br(i)%SpeedS0) &
            -dyn%br(i)%dx0*dyn%br(i)%FrcR/2.0_dp )/(vitStar(i)-dyn%br(i)%SpeedS0)
        ETotStar(i)=(rCorrStar(i)-p_star(i))/rhoStar(i)+0.5_dp*vitStar(i)**2 &
            +(dyn%br(i)%q0-dyn%br(i)%qBar)
    else
        ETotStar(i)=dyn%br(i)%ETot0+(dyn%br(i)%p0*dyn%br(i)%vit0-p_star(i) &
            *vitStar(i))/(dyn%br(i)%rho0*(dyn%br(i)%vit0-dyn%br(i)%SpeedS0)) &
            +(dyn%br(i)%q0-dyn%br(i)%qBar)
    endif
    eStar(i)=ETotStar(i)-0.5_dp*vitStar(i)**2
    EnthalpyStar(i)=ETotStar(i)+p_star(i)/rhoStar(i)

    EnthalpyMixStarNum = sum(dyn%br%AreaJGB * rhoStar * vitStar * EnthalpyStar, im)
    denominator = sum(dyn%br%AreaJGB * rhoStar * vitStar, im)

    if(abs(denominator)<1.0e-10_dp) denominator = 1.0e-10_dp
    EnthalpyMixStar = EnthalpyMixStarNum / denominator

    ! Outgoing quantities as a function of p_star(j) and correction term alpha(j)
    ! through the contact discontinuity

    ! Approximate jump relations through the shock wave
    vitStar(o)=dyn%br(o)%vit0-(p_star(o)-dyn%br(o)%p0-dyn%br(o)%rho0*(dyn%br(o)%q0 &
        -dyn%br(o)%qBar) + dyn%br(o)%dx0*dyn%br(o)%Frc/2.0_dp)/(dyn%br(o)%rho0 &
        *(dyn%br(o)%vit0-dyn%br(o)%SpeedS0))
    rhoStar_Int(o)=dyn%br(o)%rho0*(dyn%br(o)%vit0-dyn%br(o)%SpeedS0)/(vitStar(o) &
        -dyn%br(o)%SpeedS0)
    ETotStar_Int(o)=dyn%br(o)%ETot0+(dyn%br(o)%p0*dyn%br(o)%vit0-p_star(o)*vitStar(o)) &
        /(dyn%br(o)%rho0*(dyn%br(o)%vit0-dyn%br(o)%SpeedS0))+(dyn%br(o)%q0-dyn%br(o)%qBar)

    QDMStar_Int(o)=rhoStar_Int(o)*vitStar(o)
    EnergyStar_Int(o)=rhoStar_Int(o)*ETotStar_Int(o)

    ! Correction through the contact discontinuity
    ! Additional term related to the real gas EOS
    do jj = 1, NbOut
        dFdRo(o(jj)) = droeint_droP(rhoStar_Int(o(jj)), p_star(o(jj)))
    enddo
    rhoStar(o)=rhoStar_Int(o)+alpha(o)
    QDMStar(o)=QDMStar_Int(o)+alpha(o)*vitStar(o)
    EnergyStar(o)=EnergyStar_Int(o)+alpha(o)*(dFdRo(o)+0.5_dp*vitStar(o)*vitStar(o))
    ETotStar(o)=EnergyStar(o)/rhoStar(o)

    EnthalpyStar(o)=(EnergyStar(o)+p_star(o))/rhoStar(o)

    eStar(o)=ETotStar(o)-0.5_dp*vitStar(o)**2
    !do j=1,NbOut
        if(R_Correction) then
            do jj = 1, NbOut
                T_temp(o(jj)) = T_roE(rhoStar(o(jj)), eStar(o(jj)))
                rCorrStar(o(jj)) = r_roT(rhoStar(o(jj)), T_temp(o(jj)))
            enddo
        endif
    !enddo

    !            Derivative computation for implicit scheme
    !..........................................................................!
    ! Jacobian matrix evaluation for the computation of derivatives of p_star
    ! and alpha with respect to 0-variables (all branches involved)
    do j=1,NbBr
        DerVitStar_dX(j)=-1.0_dp/(dyn%br(j)%rho0*(dyn%br(j)%vit0-dyn%br(j)%SpeedS0))
        DerRhoStar_dX(j)=-(dyn%br(j)%rho0*(dyn%br(j)%vit0-dyn%br(j)%SpeedS0)) &
            *DerVitStar_dX(j)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)

        if(j<=NbOut) then
            DerETotStar_Int_dX(j)=-(vitStar(j)+p_star(j)*DerVitStar_dX(j)) &
                /(dyn%br(j)%rho0*(dyn%br(j)%vit0-dyn%br(j)%SpeedS0))
            Der_dFdRo_dRhoInt(j)=(droeint_droP(rhoStar_Int(j)+1.0e-2_dp,p_star(j)) &
                -droeint_droP(rhoStar_Int(j),p_star(j)))/1.0e-2_dp         
            Der_dFdRo_dP(j)=(droeint_droP(rhoStar_Int(j),p_star(j)+1.0e-1_dp) &
                -droeint_droP(rhoStar_Int(j),p_star(j)))/1.0e-1_dp
            DerRhoETotStar_dX(j)=DerRhoStar_dX(j)*ETotStar_Int(j)+rhoStar_Int(j) &
                *DerETotStar_Int_dX(j)+alpha(j)*(Der_dFdRo_dRhoInt(j)*DerRhoStar_dX(j) &
                +Der_dFdRo_dP(j)+vitStar(j)*DerVitStar_dX(j))
            DerHStar_dX(j)=((DerRhoETotStar_dX(j)+1.0_dp)*rhoStar(j)-(rhoStar(j)*ETotStar(j) &
                +p_star(j))*DerRhoStar_dX(j))/(rhoStar(j)**2)
            DerRhoETotStar_dAlpha(j)=(dFdRo(j)+0.5_dp*vitStar(j)*vitStar(j))
            DerHStar_dAlpha(j)=(DerRhoETotStar_dAlpha(j)*rhoStar(j)-(rhoStar(j)*ETotStar(j) &
                +p_star(j)))/(rhoStar(j)**2)
        else
            if(R_Correction) then
                DerRCorrStar_dX(j)=-(dyn%br(j)%rCorr0*(dyn%br(j)%vit0-dyn%br(j)%SpeedS0) &
                    -dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp )*DerVitStar_dX(j)/((vitStar(j) &
                    -dyn%br(j)%SpeedS0)**2)
                DerETotStar_dX(j)=((DerRCorrStar_dX(j)-1.0_dp)*rhoStar(j)-(rCorrStar(j) &
                    -p_star(j))*DerRhoStar_dX(j))/(rhoStar(j)**2)+vitStar(j)*DerVitStar_dX(j)
            else
                DerETotStar_dX(j)=-(vitStar(j)+p_star(j)*DerVitStar_dX(j))/(dyn%br(j)%rho0 &
                    *(dyn%br(j)%vit0-dyn%br(j)%SpeedS0))
            endif
            DerHStar_dX(j)=DerETotStar_dX(j)+(rhoStar(j)-p_star(j)*DerRhoStar_dX(j))/(rhoStar(j)**2)
        endif
    enddo
    
    DerEntMixStar_dX(1:NbOut)=0.0_dp
    do j=NbOut+1,NbBr
        DerEntMixStar_dX(j)=((dyn%br(j)%AreaJGB*(DerRhoStar_dX(j)*vitStar(j)+rhoStar(j) &
            *DerVitStar_dX(j))*EnthalpyStar(j)+dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j) &
            *DerHStar_dX(j))*denominator-EnthalpyMixStarNum*dyn%br(j)%AreaJGB*(DerRhoStar_dX(j) &
            *vitStar(j)+rhoStar(j)*DerVitStar_dX(j)))/(denominator**2)
    enddo
  
    fJac = 0.0_dp
    do j=1,NbBr+NbOut 
        if(j<=NbBr) then
            fJac(j,1)=dyn%br(j)%AreaJGB*(DerRhoStar_dX(j)*vitStar(j)+rhoStar(j)*DerVitStar_dX(j))
        else
            fJac(j,1)=dyn%br(j-NbBr)%AreaJGB*vitStar(j-NbBr)
        endif
    enddo

    do j=1,NbOut  ! Outgoing pipes
        if(abs(rhoStar(j)*vitStar(j))<epsCoef.or.abs(rhoStar(NbOut+1)*vitStar(NbOut+1))<epsCoef)then
            fJac(j,j+1)=-1.0_dp
            fJac(NbOut+1,j+1)=1.0_dp
            fJac(NbBr+j,j+1)=0.0_dp
        else
            qj=-dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j)/(dyn%br(NbOut+1)%AreaJGB &
                *rhoStar(NbOut+1)*vitStar(NbOut+1))
            ps_ij=dyn%br(NbOut+1)%AreaJGB/dyn%br(j)%AreaJGB

            Der_qjdX=-(dyn%br(j)%AreaJGB/(dyn%br(NbOut+1)%AreaJGB*rhoStar(NbOut+1) &
                *vitStar(NbOut+1)))*(DerRhoStar_dX(j)*vitStar(j)+rhoStar(j)*DerVitStar_dX(j))        
            Der_qRef_dX=(dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j)/dyn%br(NbOut+1)%AreaJGB) &
                *(DerRhoStar_dX(NbOut+1)*vitStar(NbOut+1)+rhoStar(NbOut+1)*DerVitStar_dX(NbOut+1)) &
                /((rhoStar(NbOut+1)*vitStar(NbOut+1))**2)                   
            CoefPLoss(j)=1.0_dp-cos((3.0_dp/4.0_dp)*(3.14159_dp-dyn%br(j)%ThetaJGB))/(qj*ps_ij)
            DerCoef_dX(j)=Der_qjdX*cos((3.0_dp/4.0_dp)*(3.14159_dp-dyn%br(j)%ThetaJGB)) &
                /(ps_ij*qj**2)
            DerCoef_dX(NbOut+1)=Der_qRef_dX*cos((3.0_dp/4.0_dp)*(3.14159_dp-dyn%br(j)%ThetaJGB)) &
                /(ps_ij*qj**2)
            Der_qjdAlpha=-dyn%br(j)%AreaJGB*vitStar(j)/(dyn%br(NbOut+1)%AreaJGB*rhoStar(NbOut+1) &
                *vitStar(NbOut+1))
            DerCoef_dAlpha=Der_qjdAlpha*cos((3.0_dp/4.0_dp)*(3.14159_dp-dyn%br(j)%ThetaJGB)) &
                /(ps_ij*qj**2)
            if(NbOut+1==dyn%idxKappaRef) then
                CoefPLoss(j)=CoefPLoss(j)*dyn%kappa
                DerCoef_dX(j)=DerCoef_dX(j)*dyn%kappa
                DerCoef_dX(NbOut+1)=DerCoef_dX(NbOut+1)*dyn%kappa
                DerCoef_dAlpha=DerCoef_dAlpha*dyn%kappa
            else
                if(j==dyn%idxKappaRef) then
                    CoefPLoss(j)=CoefPLoss(j)*dyn%kappa
                    DerCoef_dX(j)=DerCoef_dX(j)*dyn%kappa
                    DerCoef_dX(NbOut+1)=DerCoef_dX(NbOut+1)*dyn%kappa
                    DerCoef_dAlpha=DerCoef_dAlpha*dyn%kappa     
                endif
            endif

            fJac(j,j+1)=-1.0_dp-DerCoef_dX(j)*rhoStar(j)*vitStar(j)**2-CoefPLoss(j) &
                *(DerRhoStar_dX(j)*vitStar(j)**2+2.0_dp*rhoStar(j)*vitStar(j)*DerVitStar_dX(j))
            fJac(NbOut+1,j+1)=1.0_dp-rhoStar(j)*(vitStar(j)**2)*DerCoef_dX(NbOut+1)
            fJac(NbBr+j,j+1)=-vitStar(j)**2*(CoefPLoss(j)+rhoStar(j)*DerCoef_dAlpha)
        endif
    enddo

    do j=NbOut+2,NbBr
        fJac(j,j)=-1.0_dp
        fJac(NbOut+1,j)=1.0_dp
    enddo

    do j=1,NbOut
        fJac(j,NbBr+j)=DerHStar_dX(j)

        fJac(NbBr+j,NbBr+j)=DerHStar_dAlpha(j)

        do jj=NbOut+1,NbBr
            fJac(jj,NbBr+j)=-DerEntMixStar_dX(jj) 
        enddo
    enddo    
   
    !..........................................................................!

    ! Jacobian inverse
    invJac(:,:)=inv(fJac)
    if (sim_error > 0) return

    do j=1,NbOut
        DerCoef(:,j)%dU1=0.0_dp
        DerCoef(:,j)%dU2=0.0_dp
        DerCoef(:,j)%dU3=0.0_dp
        DerCoef(:,j)%dU4=0.0_dp
    enddo

    do j=1,NbBr
        if(R_Correction) then
            BB=p_star(j)-(dyn%br(j)%U4-dyn%br(j)%U3+0.5_dp*(dyn%br(j)%U2**2)/dyn%br(j)%U1)
            BB=BB-dyn%br(j)%U1*(dyn%br(j)%q0-dyn%br(j)%qBar)    + dyn%br(j)%dx0*dyn%br(j)%Frc/2.0_dp 
            DerBBdUk(1)=0.5_dp*(dyn%br(j)%U2**2)/(dyn%br(j)%U1**2)-(dyn%br(j)%q0-dyn%br(j)%qBar) &
                                                           + dyn%br(j)%dx0*dyn%br(j)%DerFric1/2.0_dp 
            DerBBdUk(2)=-dyn%br(j)%U2/dyn%br(j)%U1         + dyn%br(j)%dx0*dyn%br(j)%DerFric2/2.0_dp
            DerBBdUk(3)=1.0_dp                             + dyn%br(j)%dx0*dyn%br(j)%DerFric3/2.0_dp
            DerBBdUk(4)=-1.0_dp                            + dyn%br(j)%dx0*dyn%br(j)%DerFric4/2.0_dp
        else
            HH=dyn%br(j)%ETot0+dyn%br(j)%p0/dyn%br(j)%rho0
            call jacobian_roT(dyn%br(j)%rho0, dyn%br(j)%T0, dedT, dPdT)
            dPde=dPdT/dedT
            kk=dPde/dyn%br(j)%rho0
            KKK=dyn%br(j)%CSound0**2+kk*(dyn%br(j)%vit0**2-HH)

            BB=p_star(j)-dyn%br(j)%p0
            BB=BB-dyn%br(j)%U1*(dyn%br(j)%q0-dyn%br(j)%qBar)    + dyn%br(j)%dx0*dyn%br(j)%Frc/2.0_dp 
            DerBBdUk(1)=-KKK-(dyn%br(j)%q0-dyn%br(j)%qBar) + dyn%br(j)%dx0*dyn%br(j)%DerFric1/2.0_dp 
            DerBBdUk(2)=kk*dyn%br(j)%vit0                  + dyn%br(j)%dx0*dyn%br(j)%DerFric2/2.0_dp  
            DerBBdUk(3)=-kk                                + dyn%br(j)%dx0*dyn%br(j)%DerFric3/2.0_dp  
            DerBBdUk(4)=0.0_dp                             + dyn%br(j)%dx0*dyn%br(j)%DerFric4/2.0_dp  
        endif

        DerVitStar_dU1(j)=(-dyn%br(j)%U2/(dyn%br(j)%U1**2)-(DerBBdUk(1)*(dyn%br(j)%U2 &
            -dyn%br(j)%U1*dyn%br(j)%SpeedS0)+BB*dyn%br(j)%SpeedS0)/((dyn%br(j)%U2 &
            -dyn%br(j)%U1*dyn%br(j)%SpeedS0)**2))
        DerVitStar_dU2(j)=(1.0_dp/dyn%br(j)%U1-(DerBBdUk(2)*(dyn%br(j)%U2-dyn%br(j)%U1&
            *dyn%br(j)%SpeedS0)-BB)/((dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)**2))
        DerVitStar_dU3(j)=(-DerBBdUk(3)/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0))
        DerVitStar_dU4(j)=(-DerBBdUk(4)/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0))

        DerRhoStar_dU1(j)=(-dyn%br(j)%SpeedS0*(vitStar(j)-dyn%br(j)%SpeedS0)-&
            (dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)*DerVitStar_dU1(j)) &
            /((vitStar(j)-dyn%br(j)%SpeedS0)**2)
        DerRhoStar_dU2(j)=((vitStar(j)-dyn%br(j)%SpeedS0)-(dyn%br(j)%U2-dyn%br(j)%U1*&
            dyn%br(j)%SpeedS0)*DerVitStar_dU2(j))/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
        DerRhoStar_dU3(j)=-(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0) &
            *DerVitStar_dU3(j)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
        DerRhoStar_dU4(j)=-(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0) &
            *DerVitStar_dU4(j)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)

        if(j<=NbOut) then
            if(R_Correction) then
                DD=dyn%br(j)%U4-dyn%br(j)%U3+0.5_dp*(dyn%br(j)%U2**2)/dyn%br(j)%U1
                DerDDdUk(1)=-0.5_dp*(dyn%br(j)%U2**2)/(dyn%br(j)%U1**2)
                DerDDdUk(2)=dyn%br(j)%U2/dyn%br(j)%U1
                DerDDdUk(3)=-1.0_dp
                DerDDdUk(4)=1.0_dp
            else
                HH=dyn%br(j)%ETot0+dyn%br(j)%p0/dyn%br(j)%rho0
                call jacobian_roT(dyn%br(j)%rho0, dyn%br(j)%T0, dedT, dPdT)
                dPde=dPdT/dedT
                kk=dPde/dyn%br(j)%rho0
                KKK=dyn%br(j)%CSound0**2+kk*(dyn%br(j)%vit0**2-HH)
                
                DD=dyn%br(j)%p0
                DerDDdUk(1)=KKK
                DerDDdUk(2)=-kk*dyn%br(j)%vit0
                DerDDdUk(3)=kk
                DerDDdUk(4)=0.0_dp
            endif
    
            DerETotStar_Int_dU1(j)=-dyn%br(j)%U3/(dyn%br(j)%U1**2)+(((DerDDdUk(1) &
                *dyn%br(j)%U2/dyn%br(j)%U1-DD*dyn%br(j)%U2/(dyn%br(j)%U1**2)) &
                -p_star(j)*DerVitStar_dU1(j))*(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0) &
                +(DD*dyn%br(j)%U2/dyn%br(j)%U1-p_star(j)*vitStar(j))*dyn%br(j)%SpeedS0) &
                /((dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)**2)
            DerETotStar_Int_dU2(j)=(((DerDDdUk(2)*dyn%br(j)%U2/dyn%br(j)%U1+DD/dyn%br(j)%U1) &
                -p_star(j)*DerVitStar_dU2(j))*(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0) &
                -(DD*dyn%br(j)%U2/dyn%br(j)%U1-p_star(j)*vitStar(j)))/((dyn%br(j)%U2 &
                -dyn%br(j)%U1*dyn%br(j)%SpeedS0)**2)
            DerETotStar_Int_dU3(j)=(1.0_dp/dyn%br(j)%U1)+(DerDDdUk(3)*dyn%br(j)%U2/dyn%br(j)%U1 &
                -p_star(j)*DerVitStar_dU3(j))/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)
            DerETotStar_Int_dU4(j)=(DerDDdUk(4)*dyn%br(j)%U2/dyn%br(j)%U1-p_star(j) &
                *DerVitStar_dU4(j))/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)

            Der_dFdRo_dRhoInt(j)=(droeint_droP(rhoStar_Int(j)+1.0e-2_dp,p_star(j))-&
                                droeint_droP(rhoStar_Int(j),p_star(j)))/1.0e-2_dp                        

            DerRhoETotStar_dUk(1)=DerRhoStar_dU1(j)*ETotStar_Int(j)+rhoStar_Int(j) &
                *DerETotStar_Int_dU1(j)+alpha(j)*(Der_dFdRo_dRhoInt(j)*DerRhoStar_dU1(j) &
                +vitStar(j)*DerVitStar_dU1(j))
            DerRhoETotStar_dUk(2)=DerRhoStar_dU2(j)*ETotStar_Int(j)+rhoStar_Int(j) &
                *DerETotStar_Int_dU2(j)+alpha(j)*(Der_dFdRo_dRhoInt(j)*DerRhoStar_dU2(j) &
                +vitStar(j)*DerVitStar_dU2(j))
            DerRhoETotStar_dUk(3)=DerRhoStar_dU3(j)*ETotStar_Int(j)+rhoStar_Int(j) &
                *DerETotStar_Int_dU3(j)+alpha(j)*(Der_dFdRo_dRhoInt(j)*DerRhoStar_dU3(j) &
                +vitStar(j)*DerVitStar_dU3(j))
            DerRhoETotStar_dUk(4)=DerRhoStar_dU4(j)*ETotStar_Int(j)+rhoStar_Int(j) &
                *DerETotStar_Int_dU4(j)+alpha(j)*(Der_dFdRo_dRhoInt(j)*DerRhoStar_dU4(j) &
                +vitStar(j)*DerVitStar_dU4(j))

            DerHStar_dU1(j)=(DerRhoETotStar_dUk(1)*rhoStar(j)-(rhoStar(j)*ETotStar(j) &
                +p_star(j))*DerRhoStar_dU1(j))/(rhoStar(j)**2)
            DerHStar_dU2(j)=(DerRhoETotStar_dUk(2)*rhoStar(j)-(rhoStar(j)*ETotStar(j) &
                +p_star(j))*DerRhoStar_dU2(j))/(rhoStar(j)**2)
            DerHStar_dU3(j)=(DerRhoETotStar_dUk(3)*rhoStar(j)-(rhoStar(j)*ETotStar(j) &
                +p_star(j))*DerRhoStar_dU3(j))/(rhoStar(j)**2)
            DerHStar_dU4(j)=(DerRhoETotStar_dUk(4)*rhoStar(j)-(rhoStar(j)*ETotStar(j) &
                +p_star(j))*DerRhoStar_dU4(j))/(rhoStar(j)**2)

        else
            if(R_Correction) then
                DerRCorrStar_dU1(j)=((-dyn%br(j)%U4*dyn%br(j)%U2/(dyn%br(j)%U1**2) &
                    -dyn%br(j)%dx0*dyn%br(j)%DerFricR1/2.0_dp)*(vitStar(j)-dyn%br(j)%SpeedS0) &
                    -(dyn%br(j)%U4*(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0) &
                    -dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp)*DerVitStar_dU1(j)) &
                    /((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                DerRCorrStar_dU2(j)=((dyn%br(j)%U4/dyn%br(j)%U1-dyn%br(j)%dx0 &
                    *dyn%br(j)%DerFricR2/2.0_dp)*(vitStar(j)-dyn%br(j)%SpeedS0) &
                    -(dyn%br(j)%U4*(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0) &
                    -dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp)*DerVitStar_dU2(j)) &
                    /((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                DerRCorrStar_dU3(j)=((-dyn%br(j)%dx0*dyn%br(j)%DerFricR3/2.0_dp)*(vitStar(j)&
                    -dyn%br(j)%SpeedS0)-(dyn%br(j)%U4*(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0) &
                    -dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp)*DerVitStar_dU3(j))/((vitStar(j) &
                    -dyn%br(j)%SpeedS0)**2)
                DerRCorrStar_dU4(j)=((dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0 &
                    -dyn%br(j)%dx0*dyn%br(j)%DerFricR4/2.0_dp)*(vitStar(j)-dyn%br(j)%SpeedS0) &
                    -(dyn%br(j)%U4*(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0) &
                    -dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp )*DerVitStar_dU4(j))/((vitStar(j) &
                    -dyn%br(j)%SpeedS0)**2)
            
                DerETotStar_dU1(j)=(DerRCorrStar_dU1(j)*rhoStar(j)-(rCorrStar(j)-p_star(j)) &
                    *DerRhoStar_dU1(j))/(rhoStar(j)**2)+vitStar(j)*DerVitStar_dU1(j)
                DerETotStar_dU2(j)=(DerRCorrStar_dU2(j)*rhoStar(j)-(rCorrStar(j)-p_star(j)) &
                    *DerRhoStar_dU2(j))/(rhoStar(j)**2)+vitStar(j)*DerVitStar_dU2(j)
                DerETotStar_dU3(j)=(DerRCorrStar_dU3(j)*rhoStar(j)-(rCorrStar(j)-p_star(j)) &
                    *DerRhoStar_dU3(j))/(rhoStar(j)**2)+vitStar(j)*DerVitStar_dU3(j)
                DerETotStar_dU4(j)=(DerRCorrStar_dU4(j)*rhoStar(j)-(rCorrStar(j)-p_star(j)) &
                    *DerRhoStar_dU4(j))/(rhoStar(j)**2)+vitStar(j)*DerVitStar_dU4(j)

                DerHStar_dU1(j)=DerETotStar_dU1(j)-p_star(j)*DerRhoStar_dU1(j)/(rhoStar(j)**2)
                DerHStar_dU2(j)=DerETotStar_dU2(j)-p_star(j)*DerRhoStar_dU2(j)/(rhoStar(j)**2)
                DerHStar_dU3(j)=DerETotStar_dU3(j)-p_star(j)*DerRhoStar_dU3(j)/(rhoStar(j)**2)
                DerHStar_dU4(j)=DerETotStar_dU4(j)-p_star(j)*DerRhoStar_dU4(j)/(rhoStar(j)**2)          
            else
                HH=dyn%br(j)%ETot0+dyn%br(j)%p0/dyn%br(j)%rho0
                call jacobian_roT(dyn%br(j)%rho0, dyn%br(j)%T0, dedT, dPdT)
                dPde=dPdT/dedT
                kk=dPde/dyn%br(j)%rho0
                KKK=dyn%br(j)%CSound0**2+kk*(dyn%br(j)%vit0**2-HH)
                
                DD=dyn%br(j)%p0
                DerDDdUk(1)=KKK
                DerDDdUk(2)=-kk*dyn%br(j)%vit0
                DerDDdUk(3)=kk
                DerDDdUk(4)=0.0_dp

                DerETotStar_dU1(j)=-dyn%br(j)%U3/(dyn%br(j)%U1**2)+(((DerDDdUk(1)*dyn%br(j)%U2 &
                    /dyn%br(j)%U1-DD*dyn%br(j)%U2/(dyn%br(j)%U1**2))-p_star(j)*DerVitStar_dU1(j)) &
                    *(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)+(DD*dyn%br(j)%U2/dyn%br(j)%U1 &
                    -p_star(j)*vitStar(j))*dyn%br(j)%SpeedS0)/((dyn%br(j)%U2-dyn%br(j)%U1 &
                    *dyn%br(j)%SpeedS0)**2)
                DerETotStar_dU2(j)=(((DerDDdUk(2)*dyn%br(j)%U2/dyn%br(j)%U1+DD/dyn%br(j)%U1) &
                    -p_star(j)*DerVitStar_dU2(j))*(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0) &
                    -(DD*dyn%br(j)%U2/dyn%br(j)%U1-p_star(j)*vitStar(j)))/((dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)**2)
                DerETotStar_dU3(j)=(1.0_dp/dyn%br(j)%U1)+(DerDDdUk(3)*dyn%br(j)%U2/dyn%br(j)%U1 &
                    -p_star(j)*DerVitStar_dU3(j))/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)
                DerETotStar_dU4(j)=(DerDDdUk(4)*dyn%br(j)%U2/dyn%br(j)%U1-p_star(j) &
                    *DerVitStar_dU4(j))/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)           

                DerHStar_dU1(j)=DerETotStar_dU1(j)-p_star(j)*DerRhoStar_dU1(j)/(rhoStar(j)**2)
                DerHStar_dU2(j)=DerETotStar_dU2(j)-p_star(j)*DerRhoStar_dU2(j)/(rhoStar(j)**2)
                DerHStar_dU3(j)=DerETotStar_dU3(j)-p_star(j)*DerRhoStar_dU3(j)/(rhoStar(j)**2)
                DerHStar_dU4(j)=DerETotStar_dU4(j)-p_star(j)*DerRhoStar_dU4(j)/(rhoStar(j)**2)
            endif
        endif
    enddo

    DerEntMixStar_dU1(1:NbOut)=0.0_dp
    DerEntMixStar_dU2(1:NbOut)=0.0_dp
    DerEntMixStar_dU3(1:NbOut)=0.0_dp
    DerEntMixStar_dU4(1:NbOut)=0.0_dp
    do j=NbOut+1,NbBr
      DerEntMixStar_dU1(j)=((dyn%br(j)%AreaJGB*(DerRhoStar_dU1(j)*vitStar(j) &
        +rhoStar(j)*DerVitStar_dU1(j))*EnthalpyStar(j)+dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j) &
        *DerHStar_dU1(j))*denominator-EnthalpyMixStarNum*dyn%br(j)%AreaJGB*(DerRhoStar_dU1(j) &
        *vitStar(j)+rhoStar(j)*DerVitStar_dU1(j)))/(denominator**2)
      DerEntMixStar_dU2(j)=((dyn%br(j)%AreaJGB*(DerRhoStar_dU2(j)*vitStar(j) &
        +rhoStar(j)*DerVitStar_dU2(j))*EnthalpyStar(j)+dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j) &
        *DerHStar_dU2(j))*denominator-EnthalpyMixStarNum*dyn%br(j)%AreaJGB*(DerRhoStar_dU2(j) &
        *vitStar(j)+rhoStar(j)*DerVitStar_dU2(j)))/(denominator**2)
      DerEntMixStar_dU3(j)=((dyn%br(j)%AreaJGB*(DerRhoStar_dU3(j)*vitStar(j) &
        +rhoStar(j)*DerVitStar_dU3(j))*EnthalpyStar(j)+dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j) &
        *DerHStar_dU3(j))*denominator-EnthalpyMixStarNum*dyn%br(j)%AreaJGB*(DerRhoStar_dU3(j) &
        *vitStar(j)+rhoStar(j)*DerVitStar_dU3(j)))/(denominator**2)
      DerEntMixStar_dU4(j)=((dyn%br(j)%AreaJGB*(DerRhoStar_dU4(j)*vitStar(j) &
        +rhoStar(j)*DerVitStar_dU4(j))*EnthalpyStar(j)+dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j) &
        *DerHStar_dU4(j))*denominator-EnthalpyMixStarNum*dyn%br(j)%AreaJGB*(DerRhoStar_dU4(j) &
        *vitStar(j)+rhoStar(j)*DerVitStar_dU4(j)))/(denominator**2)
    enddo

    do j=1,NbOut
      if(abs(rhoStar(j)*vitStar(j))<epsCoef.or.abs(rhoStar(NbOut+1)*vitStar(NbOut+1))<epsCoef) then
        CoefPLoss(j)=0.0_dp
        DerCoef(j,j)%dU1=0.0_dp
        DerCoef(j,j)%dU2=0.0_dp
        DerCoef(j,j)%dU3=0.0_dp
        DerCoef(j,j)%dU4=0.0_dp
        DerCoef(NbOut+1,j)%dU1=0.0_dp
        DerCoef(NbOut+1,j)%dU2=0.0_dp
        DerCoef(NbOut+1,j)%dU3=0.0_dp
        DerCoef(NbOut+1,j)%dU4=0.0_dp
      else
        qj=-dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j)/(dyn%br(NbOut+1)%AreaJGB*rhoStar(NbOut+1) &
            *vitStar(NbOut+1))
        ps_ij=dyn%br(NbOut+1)%AreaJGB/dyn%br(j)%AreaJGB

        Der_qj_dUk(1)=-(dyn%br(j)%AreaJGB/(dyn%br(NbOut+1)%AreaJGB*rhoStar(NbOut+1) &
            *vitStar(NbOut+1)))*(DerRhoStar_dU1(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU1(j)) 
        Der_qj_dUk(2)=-(dyn%br(j)%AreaJGB/(dyn%br(NbOut+1)%AreaJGB*rhoStar(NbOut+1) &
            *vitStar(NbOut+1)))*(DerRhoStar_dU2(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU2(j)) 
        Der_qj_dUk(3)=-(dyn%br(j)%AreaJGB/(dyn%br(NbOut+1)%AreaJGB*rhoStar(NbOut+1) &
            *vitStar(NbOut+1)))*(DerRhoStar_dU3(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU3(j)) 
        Der_qj_dUk(4)=-(dyn%br(j)%AreaJGB/(dyn%br(NbOut+1)%AreaJGB*rhoStar(NbOut+1) &
            *vitStar(NbOut+1)))*(DerRhoStar_dU4(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU4(j))           

        Der_qRef_dUk(1)=(dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j)/dyn%br(NbOut+1)%AreaJGB) &
            *(DerRhoStar_dU1(NbOut+1)*vitStar(NbOut+1)+rhoStar(NbOut+1)*DerVitStar_dU1(NbOut+1)) &
            /((rhoStar(NbOut+1)*vitStar(NbOut+1))**2)
        Der_qRef_dUk(2)=(dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j)/dyn%br(NbOut+1)%AreaJGB) &
            *(DerRhoStar_dU2(NbOut+1)*vitStar(NbOut+1)+rhoStar(NbOut+1)*DerVitStar_dU2(NbOut+1)) &
            /((rhoStar(NbOut+1)*vitStar(NbOut+1))**2)
        Der_qRef_dUk(3)=(dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j)/dyn%br(NbOut+1)%AreaJGB) &
            *(DerRhoStar_dU3(NbOut+1)*vitStar(NbOut+1)+rhoStar(NbOut+1)*DerVitStar_dU3(NbOut+1)) &
            /((rhoStar(NbOut+1)*vitStar(NbOut+1))**2)
        Der_qRef_dUk(4)=(dyn%br(j)%AreaJGB*rhoStar(j)*vitStar(j)/dyn%br(NbOut+1)%AreaJGB) &
            *(DerRhoStar_dU4(NbOut+1)*vitStar(NbOut+1)+rhoStar(NbOut+1)*DerVitStar_dU4(NbOut+1)) &
            /((rhoStar(NbOut+1)*vitStar(NbOut+1))**2)                      

        CoefPLoss(j)=1.0_dp-cos((3.0_dp/4.0_dp)*(3.14159_dp-dyn%br(j)%ThetaJGB))/(qj*ps_ij)
        DerCoef(j,j)%dU1=Der_qj_dUk(1)*cos((3.0_dp/4.0_dp)*(3.14159_dp &
            -dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)
        DerCoef(j,j)%dU2=Der_qj_dUk(2)*cos((3.0_dp/4.0_dp)*(3.14159_dp &
            -dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)
        DerCoef(j,j)%dU3=Der_qj_dUk(3)*cos((3.0_dp/4.0_dp)*(3.14159_dp &
            -dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)
        DerCoef(j,j)%dU4=Der_qj_dUk(4)*cos((3.0_dp/4.0_dp)*(3.14159_dp &
            -dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)
        DerCoef(NbOut+1,j)%dU1=Der_qRef_dUk(1)*cos((3.0_dp/4.0_dp) &
            *(3.14159_dp-dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)
        DerCoef(NbOut+1,j)%dU2=Der_qRef_dUk(2)*cos((3.0_dp/4.0_dp) &
            *(3.14159_dp-dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)
        DerCoef(NbOut+1,j)%dU3=Der_qRef_dUk(3)*cos((3.0_dp/4.0_dp) &
            *(3.14159_dp-dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)
        DerCoef(NbOut+1,j)%dU4=Der_qRef_dUk(4)*cos((3.0_dp/4.0_dp) &
            *(3.14159_dp-dyn%br(j)%ThetaJGB))/(ps_ij*qj**2)

        if(NbOut+1==dyn%idxKappaRef) then
          CoefPLoss(j)=CoefPLoss(j)*dyn%kappa
          DerCoef(j,j)%dU1=DerCoef(j,j)%dU1*dyn%kappa
          DerCoef(j,j)%dU2=DerCoef(j,j)%dU2*dyn%kappa
          DerCoef(j,j)%dU3=DerCoef(j,j)%dU3*dyn%kappa
          DerCoef(j,j)%dU4=DerCoef(j,j)%dU4*dyn%kappa
          DerCoef(NbOut+1,j)%dU1=DerCoef(NbOut+1,j)%dU1*dyn%kappa
          DerCoef(NbOut+1,j)%dU2=DerCoef(NbOut+1,j)%dU2*dyn%kappa
          DerCoef(NbOut+1,j)%dU3=DerCoef(NbOut+1,j)%dU3*dyn%kappa
          DerCoef(NbOut+1,j)%dU4=DerCoef(NbOut+1,j)%dU4*dyn%kappa
        else
          if(j==dyn%idxKappaRef) then
            CoefPLoss(j)=CoefPLoss(j)*dyn%kappa
            DerCoef(j,j)%dU1=DerCoef(j,j)%dU1*dyn%kappa
            DerCoef(j,j)%dU2=DerCoef(j,j)%dU2*dyn%kappa
            DerCoef(j,j)%dU3=DerCoef(j,j)%dU3*dyn%kappa
            DerCoef(j,j)%dU4=DerCoef(j,j)%dU4*dyn%kappa
            DerCoef(NbOut+1,j)%dU1=DerCoef(NbOut+1,j)%dU1*dyn%kappa
            DerCoef(NbOut+1,j)%dU2=DerCoef(NbOut+1,j)%dU2*dyn%kappa
            DerCoef(NbOut+1,j)%dU3=DerCoef(NbOut+1,j)%dU3*dyn%kappa
            DerCoef(NbOut+1,j)%dU4=DerCoef(NbOut+1,j)%dU4*dyn%kappa              
          endif
        endif

      endif
    enddo

    do j=1,NbBr  
        DerFct(j,1)%dU1=dyn%br(j)%AreaJGB*(DerRhoStar_dU1(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU1(j))   
        DerFct(j,1)%dU2=dyn%br(j)%AreaJGB*(DerRhoStar_dU2(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU2(j))
        DerFct(j,1)%dU3=dyn%br(j)%AreaJGB*(DerRhoStar_dU3(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU3(j))
        DerFct(j,1)%dU4=dyn%br(j)%AreaJGB*(DerRhoStar_dU4(j)*vitStar(j)+rhoStar(j)*DerVitStar_dU4(j))      
    enddo

    do j=1,NbOut  ! Outgoing pipes
        DerFct(j,j+1)%dU1=-DerCoef(j,j)%dU1*rhoStar(j)*vitStar(j)**2-CoefPLoss(j)*&
                            (DerRhoStar_dU1(j)*vitStar(j)**2+2.0_dp*rhoStar(j)*vitStar(j)*&
                            DerVitStar_dU1(j))
        DerFct(j,j+1)%dU2=-DerCoef(j,j)%dU2*rhoStar(j)*vitStar(j)**2-CoefPLoss(j)*&
                            (DerRhoStar_dU2(j)*vitStar(j)**2+2.0_dp*rhoStar(j)*vitStar(j)*&
                            DerVitStar_dU2(j))
        DerFct(j,j+1)%dU3=-DerCoef(j,j)%dU3*rhoStar(j)*vitStar(j)**2-CoefPLoss(j)*&
                            (DerRhoStar_dU3(j)*vitStar(j)**2+2.0_dp*rhoStar(j)*vitStar(j)*&
                            DerVitStar_dU3(j))
        DerFct(j,j+1)%dU4=-DerCoef(j,j)%dU4*rhoStar(j)*vitStar(j)**2-CoefPLoss(j)*&
                            (DerRhoStar_dU4(j)*vitStar(j)**2+2.0_dp*rhoStar(j)*vitStar(j)*&
                            DerVitStar_dU4(j))

        DerFct(NbOut+1,j+1)%dU1=-rhoStar(j)*(vitStar(j)**2)*DerCoef(NbOut+1,j)%dU1
        DerFct(NbOut+1,j+1)%dU2=-rhoStar(j)*(vitStar(j)**2)*DerCoef(NbOut+1,j)%dU2
        DerFct(NbOut+1,j+1)%dU3=-rhoStar(j)*(vitStar(j)**2)*DerCoef(NbOut+1,j)%dU3
        DerFct(NbOut+1,j+1)%dU4=-rhoStar(j)*(vitStar(j)**2)*DerCoef(NbOut+1,j)%dU4
    enddo

    do j=NbOut+2,NbBr  ! Incoming pipes
        DerFct(j,j)%dU1=0.0_dp
        DerFct(j,j)%dU2=0.0_dp
        DerFct(j,j)%dU3=0.0_dp
        DerFct(j,j)%dU4=0.0_dp

        DerFct(NbOut+1,j)%dU1=0.0_dp
        DerFct(NbOut+1,j)%dU2=0.0_dp
        DerFct(NbOut+1,j)%dU3=0.0_dp
        DerFct(NbOut+1,j)%dU4=0.0_dp
    enddo

    do j=1,NbOut
        DerFct(j,NbBr+j)%dU1=DerHStar_dU1(j)
        DerFct(j,NbBr+j)%dU2=DerHStar_dU2(j)
        DerFct(j,NbBr+j)%dU3=DerHStar_dU3(j)
        DerFct(j,NbBr+j)%dU4=DerHStar_dU4(j)

        do jj=NbOut+1,NbBr
            DerFct(jj,NbBr+j)%dU1=-DerEntMixStar_dU1(jj)
            DerFct(jj,NbBr+j)%dU2=-DerEntMixStar_dU2(jj)
            DerFct(jj,NbBr+j)%dU3=-DerEntMixStar_dU3(jj)
            DerFct(jj,NbBr+j)%dU4=-DerEntMixStar_dU4(jj)
        enddo
    enddo

    do j=1,NbBr
        DerXStar(j,:)%dU1=-matmul(DerFct(j,:)%dU1,invJac(:,:))
        DerXStar(j,:)%dU2=-matmul(DerFct(j,:)%dU2,invJac(:,:))
        DerXStar(j,:)%dU3=-matmul(DerFct(j,:)%dU3,invJac(:,:))
        DerXStar(j,:)%dU4=-matmul(DerFct(j,:)%dU4,invJac(:,:))
    enddo

    do j=1,NbBr    !! Variables written in the junction referential
        DerPStarGlob(:,j)%dU1=DerXStar(:,j)%dU1
        DerPStarGlob(:,j)%dU2=DerXStar(:,j)%dU2
        DerPStarGlob(:,j)%dU3=DerXStar(:,j)%dU3
        DerPStarGlob(:,j)%dU4=DerXStar(:,j)%dU4
    enddo

    do j=1,NbOut    !! Variables written in the junction referential
        DerAlphaGlob(:,j)%dU1=DerXStar(:,NbBr+j)%dU1
        DerAlphaGlob(:,j)%dU2=DerXStar(:,NbBr+j)%dU2
        DerAlphaGlob(:,j)%dU3=DerXStar(:,NbBr+j)%dU3
        DerAlphaGlob(:,j)%dU4=DerXStar(:,NbBr+j)%dU4
    enddo

    do j=1,NbBr
        do jj=1,NbBr
            if(jj==j) then
                if(R_Correction) then
                    BB=p_star(j)-(dyn%br(j)%U4-dyn%br(j)%U3+0.5_dp*(dyn%br(j)%U2**2)/dyn%br(j)%U1)
                    BB=BB-dyn%br(j)%U1*(dyn%br(j)%q0-dyn%br(j)%qBar) &
                                                                + dyn%br(j)%dx0*dyn%br(j)%Frc/2.0_dp
                    DerBBdUk(1)=DerPStarGlob(j,j)%dU1+0.5_dp*(dyn%br(j)%U2**2)/(dyn%br(j)%U1**2) &
                        -(dyn%br(j)%q0-dyn%br(j)%qBar) +     dyn%br(j)%dx0*dyn%br(j)%DerFric1/2.0_dp
                    DerBBdUk(2)=DerPStarGlob(j,j)%dU2-dyn%br(j)%U2/dyn%br(j)%U1 &
                                                           + dyn%br(j)%dx0*dyn%br(j)%DerFric2/2.0_dp
                    DerBBdUk(3)=DerPStarGlob(j,j)%dU3+1.0_dp &
                                                           + dyn%br(j)%dx0*dyn%br(j)%DerFric3/2.0_dp
                    DerBBdUk(4)=DerPStarGlob(j,j)%dU4-1.0_dp &
                                                           + dyn%br(j)%dx0*dyn%br(j)%DerFric4/2.0_dp
                else
                    HH=dyn%br(j)%ETot0+dyn%br(j)%p0/dyn%br(j)%rho0
                    call jacobian_roT(dyn%br(j)%rho0, dyn%br(j)%T0, dedT, dPdT)
                    dPde=dPdT/dedT
                    kk=dPde/dyn%br(j)%rho0
                    KKK=dyn%br(j)%CSound0**2+kk*(dyn%br(j)%vit0**2-HH)

                    BB=p_star(j)-dyn%br(j)%p0
                    BB=BB-dyn%br(j)%U1*(dyn%br(j)%q0-dyn%br(j)%qBar) &
                                                                + dyn%br(j)%dx0*dyn%br(j)%Frc/2.0_dp
                    DerBBdUk(1)=DerPStarGlob(j,j)%dU1-KKK-(dyn%br(j)%q0-dyn%br(j)%qBar) &
                                                           + dyn%br(j)%dx0*dyn%br(j)%DerFric1/2.0_dp
                    DerBBdUk(2)=DerPStarGlob(j,j)%dU2+kk*dyn%br(j)%vit0 &
                                                           + dyn%br(j)%dx0*dyn%br(j)%DerFric2/2.0_dp
                    DerBBdUk(3)=DerPStarGlob(j,j)%dU3-kk   + dyn%br(j)%dx0*dyn%br(j)%DerFric3/2.0_dp
                    DerBBdUk(4)=DerPStarGlob(j,j)%dU4      + dyn%br(j)%dx0*dyn%br(j)%DerFric4/2.0_dp
                endif
            
                DerVelStarGlob(jj,j)%dU1=(-dyn%br(j)%U2/(dyn%br(j)%U1**2)-(DerBBdUk(1) &
                    *(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)+BB*dyn%br(j)%SpeedS0) &
                    /((dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)**2))
                DerVelStarGlob(jj,j)%dU2=(1.0_dp/dyn%br(j)%U1-(DerBBdUk(2)*(dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)-BB)/((dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)**2))
                DerVelStarGlob(jj,j)%dU3=(-DerBBdUk(3)/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0))
                DerVelStarGlob(jj,j)%dU4=(-DerBBdUk(4)/(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0))

                DerRhoStarGlob(jj,j)%dU1=(-dyn%br(j)%SpeedS0*(vitStar(j)-dyn%br(j)%SpeedS0) &
                    -(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0)*DerVelStarGlob(jj,j)%dU1) &
                    /((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                DerRhoStarGlob(jj,j)%dU2=((vitStar(j)-dyn%br(j)%SpeedS0)-(dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)*DerVelStarGlob(jj,j)%dU2)/((vitStar(j) &
                    -dyn%br(j)%SpeedS0)**2)
                DerRhoStarGlob(jj,j)%dU3=-(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0) &
                    *DerVelStarGlob(jj,j)%dU3/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                DerRhoStarGlob(jj,j)%dU4=-(dyn%br(j)%U2-dyn%br(j)%U1*dyn%br(j)%SpeedS0) &
                    *DerVelStarGlob(jj,j)%dU4/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
            else
                DerVelStarGlob(jj,j)%dU1=-DerPStarGlob(jj,j)%dU1/(dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)
                DerVelStarGlob(jj,j)%dU2=-DerPStarGlob(jj,j)%dU2/(dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)
                DerVelStarGlob(jj,j)%dU3=-DerPStarGlob(jj,j)%dU3/(dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)
                DerVelStarGlob(jj,j)%dU4=-DerPStarGlob(jj,j)%dU4/(dyn%br(j)%U2 &
                    -dyn%br(j)%U1*dyn%br(j)%SpeedS0)

                DerRhoStarGlob(jj,j)%dU1=-DerVelStarGlob(jj,j)%dU1*(dyn%br(j)%U2-dyn%br(j)%U1 &
                    *dyn%br(j)%SpeedS0)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                DerRhoStarGlob(jj,j)%dU2=-DerVelStarGlob(jj,j)%dU2*(dyn%br(j)%U2-dyn%br(j)%U1 &
                    *dyn%br(j)%SpeedS0)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                DerRhoStarGlob(jj,j)%dU3=-DerVelStarGlob(jj,j)%dU3*(dyn%br(j)%U2-dyn%br(j)%U1 &
                    *dyn%br(j)%SpeedS0)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                DerRhoStarGlob(jj,j)%dU4=-DerVelStarGlob(jj,j)%dU4*(dyn%br(j)%U2-dyn%br(j)%U1 &
                    *dyn%br(j)%SpeedS0)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)           
                                                    
            endif
            if(j<=NbOut) then
                DerRhoStarGlob(jj,j)%dU1=DerRhoStarGlob(jj,j)%dU1+DerAlphaGlob(jj,j)%dU1
                DerRhoStarGlob(jj,j)%dU2=DerRhoStarGlob(jj,j)%dU2+DerAlphaGlob(jj,j)%dU2
                DerRhoStarGlob(jj,j)%dU3=DerRhoStarGlob(jj,j)%dU3+DerAlphaGlob(jj,j)%dU3
                DerRhoStarGlob(jj,j)%dU4=DerRhoStarGlob(jj,j)%dU4+DerAlphaGlob(jj,j)%dU4
            endif   
        enddo
    enddo

    if(R_Correction) then
        do j=1,NbBr
            do jj=1,NbBr
                if(j<=NbOut) then
                    T_temp(j)=T_roE(rhoStar(j),eStar(j))
                    call jacobian_roT(rhoStar(j), T_temp(j), dedT, dPdT, dTdp_Ro, dTdRo_p, dRStar_dRo, dRStar_dT, RStar)
                    DerRCorStarGlob(jj,j)%dU1=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                        *DerRhoStarGlob(jj,j)%dU1+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU1
                    DerRCorStarGlob(jj,j)%dU2=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                        *DerRhoStarGlob(jj,j)%dU2+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU2  
                    DerRCorStarGlob(jj,j)%dU3=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                        *DerRhoStarGlob(jj,j)%dU3+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU3
                    DerRCorStarGlob(jj,j)%dU4=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                        *DerRhoStarGlob(jj,j)%dU4+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU4
                else
                    if(j==jj) then
                        DerRCorStarGlob(jj,j)%dU1=((-dyn%br(j)%U4*dyn%br(j)%U2/(dyn%br(j)%U1**2) &
                            -dyn%br(j)%dx0*dyn%br(j)%DerFricR1/2.0_dp)*(vitStar(j) &
                            -dyn%br(j)%SpeedS0)-(dyn%br(j)%U4*(dyn%br(j)%U2/dyn%br(j)%U1 &
                            -dyn%br(j)%SpeedS0)-dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp) &
                            *DerVelStarGlob(jj,j)%dU1)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                        DerRCorStarGlob(jj,j)%dU2=((dyn%br(j)%U4/dyn%br(j)%U1-dyn%br(j)%dx0 &
                            *dyn%br(j)%DerFricR2/2.0_dp)*(vitStar(j)-dyn%br(j)%SpeedS0) &
                            -(dyn%br(j)%U4*(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0) &
                            -dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp)*DerVelStarGlob(jj,j)%dU2) &
                            /((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                        DerRCorStarGlob(jj,j)%dU3=((-dyn%br(j)%dx0*dyn%br(j)%DerFricR3/2.0_dp) &
                            *(vitStar(j)-dyn%br(j)%SpeedS0)-(dyn%br(j)%U4*(dyn%br(j)%U2 &
                            /dyn%br(j)%U1-dyn%br(j)%SpeedS0)-dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp) &
                            *DerVelStarGlob(jj,j)%dU3)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                        DerRCorStarGlob(jj,j)%dU4=((dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0 &
                            -dyn%br(j)%dx0*dyn%br(j)%DerFricR4/2.0_dp)*(vitStar(j) &
                            -dyn%br(j)%SpeedS0)-(dyn%br(j)%U4*(dyn%br(j)%U2/dyn%br(j)%U1 &
                            -dyn%br(j)%SpeedS0) - dyn%br(j)%dx0*dyn%br(j)%FrcR/2.0_dp) &
                            *DerVelStarGlob(jj,j)%dU4)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                    else
                        DerRCorStarGlob(jj,j)%dU1=-DerVelStarGlob(jj,j)%dU1*(dyn%br(j)%U4 &
                            *(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0)-dyn%br(j)%dx0 &
                            *dyn%br(j)%FrcR/2.0_dp)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                        DerRCorStarGlob(jj,j)%dU2=-DerVelStarGlob(jj,j)%dU2*(dyn%br(j)%U4 &
                            *(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0)-dyn%br(j)%dx0 &
                            *dyn%br(j)%FrcR/2.0_dp)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                        DerRCorStarGlob(jj,j)%dU3=-DerVelStarGlob(jj,j)%dU3*(dyn%br(j)%U4 &
                            *(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0)-dyn%br(j)%dx0 &
                            *dyn%br(j)%FrcR/2.0_dp)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)
                        DerRCorStarGlob(jj,j)%dU4=-DerVelStarGlob(jj,j)%dU4*(dyn%br(j)%U4 &
                            *(dyn%br(j)%U2/dyn%br(j)%U1-dyn%br(j)%SpeedS0)-dyn%br(j)%dx0 &
                            *dyn%br(j)%FrcR/2.0_dp)/((vitStar(j)-dyn%br(j)%SpeedS0)**2)  
                    endif
                endif        
            enddo
        enddo
    else
        do j=1,NbBr
            do jj=1,NbBr
                T_temp(j)=T_roE(rhoStar(j),eStar(j))
                call jacobian_roT(rhoStar(j), T_temp(j), dedT, dPdT, dTdp_Ro, dTdRo_p, dRStar_dRo, dRStar_dT, RStar)
                DerRPhysStarGlob(jj,j)%dU1=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                    *DerRhoStarGlob(jj,j)%dU1+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU1
                DerRPhysStarGlob(jj,j)%dU2=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                    *DerRhoStarGlob(jj,j)%dU2+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU2  
                DerRPhysStarGlob(jj,j)%dU3=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                    *DerRhoStarGlob(jj,j)%dU3+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU3
                DerRPhysStarGlob(jj,j)%dU4=(dRStar_dRo+dRStar_dT*dTdRo_p) &
                    *DerRhoStarGlob(jj,j)%dU4+dRStar_dT*dTdp_Ro*DerPStarGlob(jj,j)%dU4
            enddo
        enddo
    endif

    do j=1,NbBr
        ConStar(j)%VarC(Con_Mas)=rhoStar(j)
        ConStar(j)%VarC(Con_Qdm)=rhoStar(j)*vitStar(j)
        ConStar(j)%VarC(Con_Ene)=rhoStar(j)*ETotStar(j)
        ConStar(j)%VarC(Con_R)=1.0_dp
        if(R_Correction) ConStar(j)%VarC(Con_R)=rCorrStar(j)

        if(R_Correction) then
            call jacobian_star_loc(ConStar(j)%VarC,Jacob_star(:,:,j))
        else
            PrimStar(j)%VarP(Pri_ro)=rhoStar(j)
            PrimStar(j)%VarP(Pri_u)=vitStar(j)
            PrimStar(j)%VarP(Pri_p)=p_star(j)
            call state_roP(rhoStar(j), p_star(j), PrimStar(j)%VarP(Pri_R), PrimStar(j)%VarP(Pri_e), &
                PrimStar(j)%VarP(Pri_T), PrimStar(j)%VarP(Pri_c))
            call jacobian_star_Req_rem_loc(PrimStar(j)%VarP,Jacob_star(:,:,j))
        endif

        dFdFJ(:,:)=0.0_dp
        dFdFJ(1,1)=dyn%br(j)%SgnJGB
        dFdFJ(2,2)=1.0_dp
        dFdFJ(3,3)=dyn%br(j)%SgnJGB
        dFdFJ(4,4)=dyn%br(j)%SgnJGB
      
        do jj=1,NbBr
            DerConStar(jj,j)%DerVarC(:,1)=[DerRhoStarGlob(jj,j)%dU1,&
                                    DerRhoStarGlob(jj,j)%dU2,&
                                    DerRhoStarGlob(jj,j)%dU3,&
                                    DerRhoStarGlob(jj,j)%dU4]

            DerConStar(jj,j)%DerVarC(:,2)=[&
                rhoStar(j)*DerVelStarGlob(jj,j)%dU1+DerRhoStarGlob(jj,j)%dU1*vitStar(j),&
                rhoStar(j)*DerVelStarGlob(jj,j)%dU2+DerRhoStarGlob(jj,j)%dU2*vitStar(j),&
                rhoStar(j)*DerVelStarGlob(jj,j)%dU3+DerRhoStarGlob(jj,j)%dU3*vitStar(j),&
                rhoStar(j)*DerVelStarGlob(jj,j)%dU4+DerRhoStarGlob(jj,j)%dU4*vitStar(j)]
        
            if(R_Correction) then
                DerConStar(jj,j)%DerVarC(:,3)=[&
                    DerRCorStarGlob(jj,j)%dU1+0.5_dp*(DerRhoStarGlob(jj,j)%dU1*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU1)-DerPStarGlob(jj,j)%dU1,&
                    DerRCorStarGlob(jj,j)%dU2+0.5_dp*(DerRhoStarGlob(jj,j)%dU2*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU2)-DerPStarGlob(jj,j)%dU2,&
                    DerRCorStarGlob(jj,j)%dU3+0.5_dp*(DerRhoStarGlob(jj,j)%dU3*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU3)-DerPStarGlob(jj,j)%dU3,&
                    DerRCorStarGlob(jj,j)%dU4+0.5_dp*(DerRhoStarGlob(jj,j)%dU4*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU4)-DerPStarGlob(jj,j)%dU4]                        
            else
                DerConStar(jj,j)%DerVarC(:,3)=[&
                    DerRPhysStarGlob(jj,j)%dU1+0.5_dp*(DerRhoStarGlob(jj,j)%dU1*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU1)-DerPStarGlob(jj,j)%dU1,&
                    DerRPhysStarGlob(jj,j)%dU2+0.5_dp*(DerRhoStarGlob(jj,j)%dU2*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU2)-DerPStarGlob(jj,j)%dU2,&
                    DerRPhysStarGlob(jj,j)%dU3+0.5_dp*(DerRhoStarGlob(jj,j)%dU3*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU3)-DerPStarGlob(jj,j)%dU3,&
                    DerRPhysStarGlob(jj,j)%dU4+0.5_dp*(DerRhoStarGlob(jj,j)%dU4*vitStar(j)*vitStar(j)+&
                        2.0_dp*rhoStar(j)*vitStar(j)*DerVelStarGlob(jj,j)%dU4)-DerPStarGlob(jj,j)%dU4]
            endif

            DerConStar(jj,j)%DerVarC(:,4)=0.0_dp
            if(R_Correction) DerConStar(jj,j)%DerVarC(:,4) = [&
                    DerRCorStarGlob(jj,j)%dU1,DerRCorStarGlob(jj,j)%dU2,&
                    DerRCorStarGlob(jj,j)%dU3,DerRCorStarGlob(jj,j)%dU4]          
        enddo

        ! Velocity star rewritten in the referential intrinsic to each branch
        vitStar(j)=dyn%br(j)%SgnJGB*vitStar(j)
        FlJ(j)%Vit=vitStar(j)
        FlJ(j)%Flx(Con_Mas)=rhoStar(j)*vitStar(j)
        FlJ(j)%Flx(Con_Qdm)=rhoStar(j)*vitStar(j)*vitStar(j)+p_star(j)
        FlJ(j)%Flx(Con_Ene)=(rhoStar(j)*ETotStar(j)+p_star(j))*vitStar(j)
        FlJ(j)%Flx(Con_R)=0.0_dp
        if(R_Correction) FlJ(j)%Flx(Con_R)=rCorrStar(j)*vitStar(j)

        do jj=1,NbBr
            dUdUjj(:,:)=0.0_dp
            dUdUjj(1,1)=1.0_dp
            dUdUjj(2,2)=dyn%br(jj)%SgnJGB
            dUdUjj(3,3)=1.0_dp
            dUdUjj(4,4)=1.0_dp

            FlJ(j)%DerFlx(:,:,jj) = matmul(dUdUjj(:,:), matmul(DerConStar(jj,j)%DerVarC(:,:),&
                    matmul(Jacob_star(:,:,j),dFdFJ(:,:))))
            if(.not. R_Correction) FlJ(j)%DerFlx(:,4,jj)=0.0_dp

            FlJ(j)%DerVit(1,jj)=dyn%br(j)%SgnJGB*DerVelStarGlob(jj,j)%dU1*dUdUjj(1,1)
            FlJ(j)%DerVit(2,jj)=dyn%br(j)%SgnJGB*DerVelStarGlob(jj,j)%dU2*dUdUjj(2,2)
            FlJ(j)%DerVit(3,jj)=dyn%br(j)%SgnJGB*DerVelStarGlob(jj,j)%dU3*dUdUjj(3,3)
            FlJ(j)%DerVit(4,jj)=dyn%br(j)%SgnJGB*DerVelStarGlob(jj,j)%dU4*dUdUjj(4,4)       
        enddo
    enddo

end subroutine solve_junction