program propep_2ph

use variaveis
implicit none


call dados

call forca_empuxo



contains
!------------------------------------------------------------------------------------------
subroutine dados

use portlib

allocate( frac_molar_gas (100), n_frac_molar_gas (100), frac_molar_cond (100), n_frac_molar_cond (100)  )
frac_molar_gas=0;frac_molar_cond=0

n_frac_molar_gas = " "
n_frac_molar_cond = " "

inquire(file='entrada.in', exist=existe)
if(.not. existe) then
    open(1,file='entrada.in')
        write(1,'(A)') " --------------- Selecione o tipo da simulação ---------------"
        write(1,'(A)') " 0     :: 0=Com fluxo de massa conhecido, calcule pressão de estagnação, " // &
         "1=Com características do motor conhecido, calcule propriedades deste motor sem divergente" // &
         "2=Com base nos resultados experimentais de velocidade de ejeção, etc.., calcule pressao de estagnação"

    close(1)
endif

    dos = system("Notepad entrada.in")


open(1,file='entrada.in')
    read(1,*)
    read(1,*) escolha 

close(1)

select case (escolha)

case (0)
    inquire(file='entrada0.txt', exist=existe)
    if(.not. existe) then
        open(1,file='entrada0.txt')
            write(1,'(A)') " --------------- Arquivo de entrada do programa de calculo da pressao ---------------"
            write(1,'(A)') " 'RD-170'     :: Título da simulação"
            write(1,'(A)') " 235.5d-3        :: Diâmetro da garganta (m) "
            write(1,'(A)') " 1d-3   :: Tolerância para o sistema iterativo da diferenca de fluxos de massa" 
            write(1,'(A)') " 3  :: Propelente ( 0=KNSU (65/35), 1=Polvora do Roberto (150/50/20) KNO3, C, S,  2=KNSU (60/40), 3=Oxigenio Liquido e Kerosene (RP-1)) "
            write(1,'(A)') " 598.25d0   :: Fluxo de massa de propelente (kg/s)"
            write(1,'(A)') " 101.325d3    :: Pressão atmosferica (Pa)"
            write(1,'(A)') " 21.39d0    :: Tempo de queima (s)"
        close(1)
    endif

    dos = system("Notepad entrada0.txt")

    open(1,file='entrada0.txt')
        read(1,*) 
        read(1,*) Titulo_sim
        read(1,*) Dg
        read(1,*) tol
        read(1,*) propelente
        read(1,*) fmp_especificado
        read(1,*) Pa
        read(1,*) tq
        
    close(1)

case (1)
    inquire(file='entrada1.txt', exist=existe)
    if(.not. existe) then
        open(1,file='entrada1.txt')
            write(1,'(A)') " --------------- Arquivo de entrada do programa de calculo da pressao ---------------"
            write(1,'(A)') " 'Dg_10mm'     :: Título da simulação"
            write(1,'(A)') " 1517.668d0   :: Massa específica do propelente (kg/m3)"
            write(1,'(A)') " 10.0d-3        :: Diâmetro da garganta (m) "
            write(1,'(A)') " 0.d0          :: Diâmetro da Alma (mm) [Se a queima for cilíndrica deixar 0.d0 !!]"
            write(1,'(A)') " 60.d-3        :: Diâmetro Externo do grão (m) "
            write(1,'(A)') " 50.d-3        :: Comprimento do grão (m) "
            write(1,'(A)') " 9.06d4    :: Pressão atmosferica (Pa)"
            write(1,'(A)') " 2.6234d-3    :: Parâmetro a da taxa de queima (m/s/Pa)"
            write(1,'(A)') " 0.d0    :: Parâmetro n da taxa de queima (em função de Pa, mas este valor é adimensional)"
            write(1,'(A)') " 1d-15   :: Tolerância para o sistema iterativo, se blocado diff Fm, se nao Diferenca do Numero de Mach" 
            write(1,'(A)') " 0  :: Propelente ( 0=KNSU (65/35), 1=Polvora do Roberto (150/50/20) KNO3, C, S,  2=KNSU (60/40)) "
        close(1)
    endif

    dos = system("Notepad entrada1.txt")

    open(1,file='entrada1.txt')
        read(1,*) 
        read(1,*) Titulo_sim
        read(1,*) ro
        read(1,*) Dg
        read(1,*) Da
        read(1,*) De
        read(1,*) Lp
        read(1,*) Pa
        read(1,*) a
        read(1,*) n
        read(1,*) tol
        read(1,*) propelente
        
    close(1)
    
case (2)
    inquire(file='entrada2.txt', exist=existe)
    if(.not. existe) then
        open(1,file='entrada2.txt')
            write(1,'(A)') " --------------- Arquivo de entrada do programa de calculo da pressao ---------------"
            write(1,'(A)') " 'C6-0_Bandeirante'     :: Título da simulação"
            write(1,'(A)') " 3.0d-3        :: Diâmetro da garganta (m) "
            write(1,'(A)') " 1d-3   :: Tolerância para o sistema iterativo da diferenca de fluxos de massa" 
            write(1,'(A)') " 1  :: Propelente ( 0=KNSU (65/35), 1=Polvora do Roberto (150/50/20) KNO3, C, S,  2=KNSU (60/40), 3=Oxigenio Liquido e Kerosene (RP-1)) "
            write(1,'(A)') " 3.18d0      :: Forca media experimental (N) "
            write(1,'(A)') " 701.d0      :: Velocidade de ejecao efetiva media dos gases, experimental (m/s) "
            write(1,'(A)') " 91.745d3    :: Pressão atmosferica (Pa)"
        close(1)
    endif

    dos = system("Notepad entrada2.txt")

    open(1,file='entrada2.txt')
        read(1,*) 
        read(1,*) Titulo_sim
        read(1,*) Dg
        read(1,*) tol
        read(1,*) propelente
        read(1,*) Fe
        read(1,*) Ve
        read(1,*) Pa
        
    close(1)
        
end select


end subroutine dados
!------------------------------------------------------------------------------------------
subroutine forca_empuxo

use portlib

real*8 :: Pi=dacos(-1.d0)
real*8 :: Ab, Ag
real*8 :: Fm_p, Fm_t

if(Da == 0.d0) then
    Ab = Pi/4*(De)**2
endif

Ag = Pi/4*(Dg)**2

! Starts cpu time
call DATE_AND_TIME(VALUES=tarray1)
!VALUES(1):  The year
!VALUES(2):  The month
!VALUES(3):  The day of the month
!VALUES(4):  Time difference with UTC in minutes
!VALUES(5):  The hour of the day
!VALUES(6):  The minutes of the hour
!VALUES(7):  The seconds of the minute
!VALUES(8):  The milliseconds of the second
! Valor em segundos absoluto, contando dias, horas, minutos, segundos e milisegundos
tcpu1 =  tarray1(3)*24.d0*60.d0*60.d0 + 60.d0*60.d0*tarray1(5)+ &
 60.d0*tarray1(6) + 1.d0*tarray1(7) + 1.d-3*tarray1(8)


open(4,file=trim(adjustl(Titulo_sim)) // ".txt")

    write(4,'(1x, A/)') ":: Arquivo de saida de dados do programa forca_empuxo_propep_2ph ::"

    write(4,'(1x, A/)') ":: DADOS INICIAIS ::"

    5 format(1PE25.14,A)
    
    select case (escolha)
    case (0)
        write(4,5)  Dg, " :: Diâmetro da garganta (m) "
        write(4,5)  tol, " :: Tolerância para o sistema iterativo" 
        write(4,'(I25,A)') propelente, ":: Propelente ( 0=KNSU (65/35), 1=Polvora do Roberto (150/50/20) KNO3, C, S,  2=KNSU (60/40))"
        write(4,5) fmp_especificado, " :: Se diferente de zero, especifica o fluxo de massa de propelente (usos diversos) (kg/s)"
        write(4,5)  Pa, " :: Pressão atmosferica (Pa)"
    case (1)
        texto = " :: Massa específica do propelente (kg/m3) "; call imprime_e3(ro, 4, texto )
        write(4,5)  Dg, " :: Diâmetro da garganta (m) "
        write(4,5)  Da, " :: Diâmetro da Alma (mm) "
        write(4,5)  De, " :: Diâmetro Externo do grão (m) "
        write(4,5)  Lp, " :: Comprimento do grão (m) "     
        write(4,5)  Pa, " :: Pressão atmosferica (Pa)"
        write(4,5)  a, " :: Parâmetro a da taxa de queima (m/s/Pa)" 
        write(4,5)  n, " :: Parâmetro n da taxa de queima (em função de Pa, mas este valor é adimensional)"   
        write(4,5)  tol, " :: Tolerância para o sistema iterativo" 
        write(4,'(I25,A)') propelente, ":: Propelente ( 0=KNSU (65/35), 1=Polvora do Roberto (150/50/20) KNO3, C, S,  2=KNSU (60/40))"
    case (2)
        write(4,5)  Dg, " :: Diâmetro da garganta (m) "
        write(4,5)  tol, " :: Tolerância para o sistema iterativo" 
        write(4,'(I25,A)') propelente, ":: Propelente ( 0=KNSU (65/35), 1=Polvora do Roberto (150/50/20) KNO3, C, S,  2=KNSU (60/40))"
        write(4,5) Fe, " :: Forca media experimental (N) "
        write(4,5) ve, " :: Velocidade de ejecao efetiva media dos gases, experimental (m/s) "
        write(4,5)  Pa, " :: Pressão atmosferica (Pa)"

    end select
    write(4,*) 
    
    write(4,'(1x, A/)') ":: Iterações ::"

! Utilizando a pressão ambiente como Chute Inicial
P0 = Pa

call calcula_pressao
! SAIDAS= R,T0,k2ph

if(escolha /= 2) then
        ! Fluxo de massa disponivel do propelente
        Fm_p = ro*Ab*a*P0**n
        if(Fmp_especificado /= 0.d0) Fm_p = Fmp_especificado


        ! Fluxo de massa necessário pela teoria para escoamento blocado
        ! Regime permanente
        Fm_t = P0*Ag*dsqrt(k2ph/(R*T0)*(2/(k2ph+1))**((k2ph+1)/(k2ph-1)))

        write(4,'(10x,A, 3X, A, 21x, A)') "It", "P0", "Variacao(Fm_p - Fm_t)"
        write(4,*) 1,P0, dabs(Fm_p-Fm_t)

        write(*,'(10x,A, 3X, A, 21x, A)') "It", "P0", "Variacao(Fm_p - Fm_t)"
        write(*,*) 1,P0, dabs(Fm_p-Fm_t)

        if(Fm_p >= Fm_t) then ! ESCOAMENTO BLOCADO!!!!!!!!!!!!!!!!!

        do ka=1, 1000000
            P0 = Fm_p/(Ag*dsqrt(k2ph/(R*T0)*(2/(k2ph+1))**((k2ph+1)/(k2ph-1))))
            
            call calcula_pressao

            ! Fluxo de massa disponivel do propelente
            Fm_p = ro*Ab*a*P0**n
            if(Fmp_especificado /= 0.d0) Fm_p = Fmp_especificado

            
            ! Fluxo de massa necessário pela teoria para escoamento blocado
            ! Regime permanente
            Fm_t = P0*Ag*dsqrt(k2ph/(R*T0)*(2/(k2ph+1))**((k2ph+1)/(k2ph-1)))


            write(*,*) ka+1,P0, dabs(Fm_p-Fm_t)
            write(4,*) ka+1,P0, dabs(Fm_p-Fm_t)
            if( dabs(Fm_p-Fm_t)< tol ) exit    

        enddo

        ! PROPRIEDADES NA GARGANTA    
        Tg = T0/(1.d0+(k2ph-1.d0)/2.d0)
        pg = P0*(1.d0+(k2ph-1.d0)/2.d0)**(-k2ph/(k2ph-1.d0))
        Vg = dsqrt(k2ph*R*Tg)
            


        elseif(Fm_p < Fm_t) then ! PRESSAO INSUFICIENTE PARA ESCOAMENTO BLOCADO NA GARGANTA!!!!

        ! Chute inicial para o numero de mach na garganta
        Mach = 0.5d0

        do k=1, 10000
            ! Propriedades na garganta:
            Tg = T0/(1.d0+(k2ph-1.d0)*0.5d0*Mach**2)
            pg = P0*(1.d0+(k2ph-1.d0)*0.5d0*Mach**2)**(-k2ph/(k2ph-1.d0))
            Vg = Fm_p*R*Tg/(Ag*pg)

            Mach2 = Vg/dsqrt(k2ph*R*Tg)

            write(*,*) k, abs(Mach2-Mach)
            if(abs(Mach2-Mach)<= tol) exit
                
            Mach = Mach2
        enddo


        !do k=1, 10000
        !    ! X = X - f/df ! NEWTON!!!!
        !    Mach = Mach2 - ( Fm_p/Ag/P0*dsqrt(R*T0/k2ph)*(1+(k2ph-1)/2*Mach**2)**((k2ph+1)/(2*(k2ph-1)))-Mach)/&
        !    (Fm_p*dsqrt(R*T0)*(1+(k2ph-1)/2*Mach**2)**((k2ph+1)/(2*k2ph-2))*(k2ph+1)*(k2ph-1)*Mach/(Ag*P0*dsqrt(k2ph)*(2*k2ph-2)*(1+(k2ph-1)/2*Mach**2))-1)
        !
        !    write(*,*) k, abs(Mach2-Mach)
        !    if(abs(Mach2-Mach)<= tol) exit
        !
        !    Mach2 = Mach
        !
        !enddo

        ! Propriedades na garganta:
        Tg = T0/(1.d0+(k2ph-1.d0)*0.5d0*Mach**2)
        pg = P0*(1.d0+(k2ph-1.d0)*0.5d0*Mach**2)**(-k2ph/(k2ph-1.d0))
        Vg = Mach*dsqrt(k2ph*R*Tg)


        endif
else

write(*,'(10x,A, 3X, A, 21x, A)') "It", "P0", "Variacao(P0_a - P0)"
write(4,'(10x,A, 3X, A, 21x, A)') "It", "P0", "Variacao(P0_a - P0)"

do k=1, 1000000      
    P0_ant = P0
    
    call calcula_pressao
    
    
    
    gamma = k2ph
    
    Mach= 2.d0*Ve/dsqrt(-2.d0*Ve**2*gamma+2.d0*Ve**2+4.d0*gamma*R*T0)


    Ps = P0*(1.d0+(gamma-1.d0)/2.d0*Mach**2)**(-gamma/(gamma-1.d0)) 


    Epsilon = 1.d0/mach*(2.d0/(gamma+1.d0)*(1.d0+(gamma-1.d0)/2.d0*mach**2))**((gamma+1.d0)/(2.d0*(gamma-1.d0)))


    Cf_vac = dsqrt((2.d0*gamma**2/(gamma-1.d0)*(2.d0/(gamma+1.d0))**((gamma+1.d0)/(gamma-1.d0))*(1.d0-(Ps/P0)**((gamma-1.d0)/gamma))))+(Ps/P0)*Epsilon;

    Cf = Cf_vac-Pa/P0*Epsilon
    
    if(Cf < 0.0d0) then
        write(*,*) "ERRO -- Valor de Cf eh menor que 0.0...., estimando o absoluto"
        Cf = abs(Cf)
    endif
    
    
    P0 = Fe*4.d0/(Cf*(pi*(Dg)**2))
    
    P0 = P0_ant + (P0 - P0_ant)*1D-1

    
        write(4,*) k,P0, dabs(P0-P0_ant)

        write(*,*) k,P0, dabs(P0-P0_ant)
    
    
    if( dabs(P0-P0_ant) <= tol ) exit

enddo

! Pos - processamento:

! Massa específica de estagnacao
ro0 = P0/R/T0

! Propriedades na saida
TS  = T0*(1.d0 + (gamma-1.d0)/2.d0*Mach**2)**(-1.d0)
roS = ro0 * (1.d0 + (gamma-1.d0)/2.d0*Mach**2)**(-1.d0/(gamma-1.d0))
Ds = dsqrt(Epsilon*Dg**2)

endif

! Ends cpu time
call DATE_AND_TIME(VALUES=tarray2)
!VALUES(1):  The year
!VALUES(2):  The month
!VALUES(3):  The day of the month
!VALUES(4):  Time difference with UTC in minutes
!VALUES(5):  The hour of the day
!VALUES(6):  The minutes of the hour
!VALUES(7):  The seconds of the minute
!VALUES(8):  The milliseconds of the second
! Valor em segundos absoluto, contando dias, horas, minutos, segundos e milisegundos
tcpu2 = tarray2(3)*24.d0*60.d0*60.d0 + 60.d0*60.d0*tarray2(5) + &
 60.d0*tarray2(6) + 1.d0*tarray2(7) + 1.d-3*tarray2(8)
 
 
tc = tcpu2 - tcpu1


! Pos - processamento:


    
    write(4,'(/1x, A/)') ":: RESULTADOS PRINCIPAIS ::"
    
    
    texto = " :: Pressao de estagnacao (Pa) "; call imprime_e3(P0, 4, texto )
    texto = " :: Pressao (Psi) "; call imprime_e0(P0/6894.8d0, 4, texto )
    texto = " :: Pressao (Atm) "; call imprime_e0(P0/101325d0, 4, texto )
    write(4,*)
    texto = " :: Temperatura adiabatica de chama (K) "; call imprime_e0(T0, 4, texto )
    write(4,*)
    texto = " :: Constante do gas (apenas gas) (J/(kg.K)) " ; call imprime_e0(R, 4, texto )
    texto = " :: Razao de calores especificos utilizada "; call imprime_e0(k2ph, 4, texto )
    write(4,*)
    
    if(escolha /=2) then
        write(4,'(/1x, A/)') ":: PROPRIEDADES NA SAIDA (GARGANTA) ::"

        texto = " :: Pressao na garganta (Pa) "; call imprime_e3(Pg, 4, texto )
        texto = " :: Temperatura na garganta(K) "; call imprime_e0(Tg, 4, texto )

        write(4,*)

        texto = " :: Velocidade dos gases na garganta (m/s) "; call imprime_e0(Vg, 4, texto )
        
        write(4,*)
        
        if(escolha == 1) then
            if(Da == 0.d0) then
                tq = Lp / (a * P0 **n )
            else
                tq = (De-Da) / (a * P0 **n )/2
            endif
            
            write(4,5)  tq, " :: Tempo de queima medio (s) "
        endif
        
        Fdin = fm_p*Vg; write(4,5)  Fdin, " :: Forca dinamica (N)"   
        
        Fp = Ag*(Pg-Pa); write(4,5)  Fp, " :: Forca de pressao (N)"  
        
        Ftotal = Fdin + Fp; write(4,5)  Ftotal, " :: Forca total (N)"  
        
        write(4,*)
        
        It = Ftotal * tq ; write(4,5)  It, " :: Impulso total (N.s)"  
        
        Is = Ftotal /fm_p/g ; write(4,5)  Is, " :: Impulso especifico (s)"
        
    else
    
    write(4,'(/1x, A/)') ":: PROPRIEDADES NA SAIDA ::"

    texto = " :: Pressao (Pa) "; call imprime_e3(Ps, 4, texto )
    texto = " :: Pressao (Psi) "; call imprime_e0(Ps/6894.8d0, 4, texto )
    texto = " :: Pressao (Atm) "; call imprime_e0(Ps/101325d0, 4, texto )

    write(4,*)

    texto = " :: Temperatura (K) "; call imprime_e0(Ts, 4, texto )

    write(4,5)  roS, " :: Massa especifica  (kg/m3) "
    write(4,5)  Cf, " :: Coeficiente de empuxo (adimensional)"
    
    write(4,*)
    texto = " :: Numero de Mach (adimensional) "; call imprime_e0(Mach, 4, texto )
    texto = " :: Razão de expansao efetiva (adimensional) "; call imprime_e0(Epsilon, 4, texto )
    texto = " :: Diametro efetivo da saida (mm) "; call imprime_e0(Ds, 4, texto )
    
    endif
        
        
        
        write(4,'(/1x, A/)') ":: RESULTADOS SECUNDARIOS ::"
    
!    endif

    
    texto = " :: Massa molar do gas (apenas gas) (g/mol) " ; call imprime_e0(Mmg, 4, texto )

    
!    if(escolha ==1) then
        write(4,5)  ro0, " :: Massa especifica do gas(kg/m3) "

        write(4,*)   
        texto = " :: Razao de calores especificos do gas "; call imprime_e0(kgas, 4, texto )
        write(4,*)
        texto = " :: Razao de calores especificos da mistura "; call imprime_e0(kmix, 4, texto )
        
!    endif

    write(4,'(2/, A, F15.3)') " Tempo computacional (s) = ", tc

close(4)

dos = system("Notepad " // trim(adjustl(Titulo_sim)) // ".txt")


end subroutine forca_empuxo
!------------------------------------------------------------------------------------------
subroutine calcula_pressao

use portlib

implicit none

iteracao = 1
if(iteracao<10) then
    write(entrada_propep,'(A,I1,A)') "input_00", iteracao , ".propep" 
    write( saida_propep,'(A,I1,A)') "output_00", iteracao , ".propep" 
elseif(iteracao<100) then
    write(entrada_propep,'(A,I2,A)') "input_0", iteracao , ".propep" 
    write( saida_propep,'(A,I2,A)') "output_0", iteracao , ".propep" 
elseif(iteracao<1000) then 
    write(entrada_propep,'(A,I3,A)') "input_", iteracao , ".propep" 
    write( saida_propep,'(A,I3,A)') "output_", iteracao , ".propep" 
endif

select case (propelente)
    case (0) ! KNSU 65/35
        open( 2,file=".\entrada\" // trim(adjustl(entrada_propep)) )
            write(2,'(A  )') "Propellant MOTOR_EXPERIMENTAL_KNSU"
            write(2,'(A  )') "+821     65 g "
            write(2,'(A  /)') "+897    35 g "
            write(2, '(A)') "HP"
            write(2, '(A,F20.10, A)') "+chamber_pressure ",P0*1d-3, " kPa"
        close(2)
        
        n_comp_prop = 2
        MASSA = 100.d0
    case (1) ! POLVORA DO ROBERTO
        open( 2,file=".\entrada\" // trim(adjustl(entrada_propep)) )
            write(2,'(A  )') "Propellant POLVORA_ROBERTO"
            write(2,'(A  )') "+821     150 g "
            write(2,'(A  )') "+245     50 g "
            write(2,'(A  /)') "+898     20 g "
            write(2, '(A)') "HP"
            write(2, '(A,F20.10, A)') "+chamber_pressure ",P0*1d-3, " kPa"
        close(2)

        n_comp_prop = 3
        MASSA = 220.d0
    case (2) ! KNSU 60/40
        open( 2,file=".\entrada\" // trim(adjustl(entrada_propep)) )
            write(2,'(A  )') "Propellant MOTOR_EXPERIMENTAL_KNSU"
            write(2,'(A  )') "+821     60 g "
            write(2,'(A  /)') "+897    40 g "
            write(2, '(A)') "HP"
            write(2, '(A,F20.10, A)') "+chamber_pressure ",P0*1d-3, " kPa"
        close(2)
        
        n_comp_prop = 2
        MASSA = 100.d0
        
    case (3) ! Liquid Oxygen and Kerosene (RP1)
        open( 2,file=".\entrada\" // trim(adjustl(entrada_propep)) )
            write(2,'(A  )') "Propellant Oxigenio_LIQUIDO_Querosene_RP1"
            write(2,'(A  )') "+737     263 g "
            write(2,'(A  /)') "+855    100 g "
            write(2, '(A)') "HP"
            write(2, '(A,F20.10, A)') "+chamber_pressure ",P0*1d-3, " kPa"
        close(2)
        
        n_comp_prop = 2
        MASSA = 263 + 100
    case (4) ! Liquid Oxygen and Kerosene (RP1)
        open( 2,file=".\entrada\" // trim(adjustl(entrada_propep)) )
            write(2,'(A  )') "Propellant Oxigenio_LIQUIDO_Querosene_RP1_L15"
            write(2,'(A  )') "+737     158 g "
            write(2,'(A  /)') "+855    100 g "
            write(2, '(A)') "HP"
            write(2, '(A,F20.10, A)') "+chamber_pressure ",P0*1d-3, " kPa"
        close(2)
        
        n_comp_prop = 2
        MASSA = 158 + 100
end select


! Usa o propep para encontrar as caracteristicas do gas a dada temperatura
dos = system("cpropep.exec -f " // ".\entrada\"// trim(adjustl(entrada_propep)) // " > " // ".\saida\" // trim(adjustl(saida_propep)))

open( 3,file= ".\saida\" // trim(adjustl(saida_propep)) )

do pulo=1, 15+n_comp_prop
    read(3,*) temp1
enddo

read(3,'(A,F15.3)') temp1, T0

do pulo=1, 16
    read(3,'(A,F15.3)')
enddo
    
do i=1, 100
    read(3,*,iostat=flag) temp1
    if(trim(adjustl(temp1)) == 'Condensed') exit
    if(flag /=0) exit
enddo
n_f_molares_gas = i-1

close(3)

    open( 3,file= ".\saida\" // trim(adjustl(saida_propep)) )
    do pulo=1, 15+n_comp_prop
        read(3,*) temp1
    enddo

    do pulo=1, 14
        read(3,*) 
    enddo
    read(3,'(A,1PE15.4)')
    read(3,'(A,1PE15.4)')
    read(3,'(A,1PE15.4)')
       
    !------------------------------------------------------------------
    ! CALCULO DA FASE GASOSA
    !------------------------------------------------------------------

    do pulo=1, n_f_molares_gas
        read(3,'(A,1PE15.4)') n_frac_molar_gas(pulo), frac_molar_gas(pulo)
    !        write(*,*) n_frac_molar_gas(pulo), frac_molar_gas(pulo)

    enddo


    call get_molar_mass( n_f_molares_gas, n_frac_molar_gas, frac_molar_gas, Mmg )
    !-------------------------------------------------------------


!------------------------------------------------------------------
! SE HOUVER FASE CONDENSADA!!!
!------------------------------------------------------------------
if(flag ==0) then 

    read(3,*) temp1

    do i=1, 100
        read(3,*,iostat=flag) temp1
        if(flag /= 0) exit
    enddo
    n_f_molares_cond = i-1
    close(3)

    ! n_f_molares_cond é o número de fases condensadas


    open( 3,file= ".\saida\" // trim(adjustl(saida_propep)) )
    do pulo=1, 15+n_comp_prop
        read(3,*) temp1
    enddo

    do pulo=1, 14
        read(3,*) 
    enddo
    read(3,'(A,1PE15.4)')
    read(3,'(A,1PE15.4)')
    read(3,'(A,1PE15.4)')

    do pulo=1, n_f_molares_gas
        read(3,*) temp1

    enddo
    read(3,*)
    do pulo=1, n_f_molares_cond
        read(3,'(A,1PE15.4)') n_frac_molar_cond(pulo), frac_molar_cond(pulo)
    !        write(*,*) n_frac_molar_cond(pulo), frac_molar_cond(pulo)
    enddo




    call get_molar_mass( n_f_molares_cond, n_frac_molar_cond, frac_molar_cond, Mmc )

       


    close(3)
    
    n_mol =MASSA/(Mmg + Mmc)


    ! MASSA MOLAR DOS GASES
    Mmg = MASSA/(( sum(frac_molar_gas) ) * n_mol)


    ! CALCULO DO CALOR ESPECIFICO DOS COMPONENTES
     call specific_heat( n_f_molares_gas, n_frac_molar_gas, frac_molar_gas, T0, Cp_gas)
     
     call specific_heat( n_f_molares_cond, n_frac_molar_cond, frac_molar_cond, T0, Cp_condensado)


    Cp_mistura = ( Cp_gas + Cp_condensado) /( sum(frac_molar_gas) )


    Cp_gas = Cp_gas /( sum(frac_molar_gas) )


    ! DOS GASES
    R = 8314.472d0 / Mmg
    kmix=Cp_mistura/(Cp_mistura - 8.314472d0)
    kgas=Cp_gas/(Cp_gas - 8.314472d0)

    qui = n_mol*Mmc/MASSA

    psi = qui/(1.d0-qui)

    Cs = Cp_condensado / sum(frac_molar_cond) 

    k2ph= kgas*(1.d0+psi*Cs/Cp_gas)/(1.d0+kgas*psi*Cs/Cp_gas)


    ! Massa específica de estagnacao
    ro0 = P0/R/T0


else ! NAO HA FASE CONDENSADA!!!!!

    close(3)
    
    n_mol =MASSA/Mmg


    ! MASSA MOLAR DOS GASES
    Mmg = MASSA/(( sum(frac_molar_gas) ) * n_mol)


    ! CALCULO DO CALOR ESPECIFICO DOS COMPONENTES
     call specific_heat( n_f_molares_gas, n_frac_molar_gas, frac_molar_gas, T0, Cp_gas)
     
    Cp_gas = Cp_gas /( sum(frac_molar_gas) )


    ! DOS GASES
    R = 8314.472d0 / Mmg

    kgas=Cp_gas/(Cp_gas - 8.314472d0)

    ! Massa específica de estagnacao
    ro0 = P0/R/T0
    
    k2ph = kgas


endif



end subroutine calcula_pressao
!------------------------------------------------------------------------------------------
subroutine imprime_e3(valor, flag, TEXTO)
real*8, intent(in) :: valor
integer, intent(in) :: flag
character*90, intent(in) :: TEXTO
integer :: inteiro
inteiro = dlog10(Valor)
if(inteiro==8) then
    write(flag,'(6PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==7) then
    write(flag,'(5PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==6) then
    write(flag,'(4PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==5) then
    write(flag,'(3PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==4) then
    write(flag,'(2PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==3) then
    write(flag,'(1PE25.14,A)')  Valor,  trim(TEXTO)
endif

end subroutine imprime_e3
!------------------------------------------------------------------------------------------
subroutine imprime_e0(valor, flag, TEXTO)
real*8, intent(in) :: valor
integer, intent(in) :: flag
character*90, intent(in) :: TEXTO
integer :: inteiro
inteiro = dlog10(Valor)
if(inteiro==8) then
    write(flag,'(9PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==7) then
    write(flag,'(8PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==6) then
    write(flag,'(7PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==5) then
    write(flag,'(6PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==4) then
    write(flag,'(5PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==3) then
    write(flag,'(4PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==2) then
    write(flag,'(3PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==1) then
    write(flag,'(2PE25.14,A)')  Valor,  trim(TEXTO)
    elseif(inteiro==0) then
    write(flag,'(1PE25.14,A)')  Valor,  trim(TEXTO)
endif

end subroutine imprime_e0
!------------------------------------------------------------------------------------------
subroutine get_molar_mass( n_f_molares, n_frac_molar, frac_molar, Mm )
integer, intent(in) :: n_f_molares

character*20, intent(in), dimension(n_f_molares) :: n_frac_molar
real*8, intent(in), dimension(n_f_molares) :: frac_molar

real*8, intent(out) :: Mm

real*8  :: Mm_componente
integer :: i

Mm = 0
do i=1, n_f_molares

    call reconhecer_Mm_componente( n_frac_molar(i), Mm_componente )
    
    Mm = Mm + Mm_componente * frac_molar(i)
enddo


end subroutine get_molar_mass
!------------------------------------------------------------------------------------------
subroutine reconhecer_Mm_componente(n_frac_molar, Mm_componente)
character*20, intent(in) :: n_frac_molar
real*8, intent(out) :: Mm_componente

integer :: quantidade_m_molar

character*15, dimension(:), allocatable :: nome_frac_molar_componente
real*8      , dimension(:), allocatable :: frac_molar_componente

integer :: i, j, flag

character*15 :: temp1

Mm_componente = 0.d0

open(1,file='tabelas\componentes_massas_molares.txt')
    do i=1, 1000
        read(1,*,iostat=flag) temp1
        if(flag /= 0) exit
    enddo
quantidade_m_molar = i-1
close(1)


allocate( nome_frac_molar_componente(quantidade_m_molar), frac_molar_componente(quantidade_m_molar) )


open(1,file='tabelas\componentes_massas_molares.txt')

do i=1, quantidade_m_molar
    read(1,'(A,1x,1PE21.15)') nome_frac_molar_componente(i), frac_molar_componente(i)
enddo

close(1)

do i=1, quantidade_m_molar
if( trim(adjustl( n_frac_molar )) == trim(adjustl( nome_frac_molar_componente(i) )) ) then
    Mm_componente = frac_molar_componente(i)
    exit
endif

enddo


if(Mm_componente == 0.d0) then

    write(*,*) " ---------------------- ERRO ------------------------ "
    write(*,*)
    write(*,*) " VALOR DA MASSA MOLECULAR DA FRACAO NAO EH CONHECIDA "
    write(*,*) n_frac_molar
    write(*,*)
    write(*,*) " APERTE QUALQUER TECLA PARA SAIR "
    pause    
    stop
endif

deallocate( nome_frac_molar_componente, frac_molar_componente )

end subroutine reconhecer_Mm_componente
!------------------------------------------------------------------------------------------
subroutine specific_heat( n_f_molares, n_frac_molar, frac_molar, T0,  Cp_total)
 
! Entradas:
! Quantidade de frações molares, nome da fração molar, valor da fração molar, Temperatura de estagnação

! Saida:
! Calor específico

integer, intent(in) :: n_f_molares

character*20, intent(in), dimension(n_f_molares) :: n_frac_molar
real*8, intent(in), dimension(n_f_molares) :: frac_molar
real*8, intent(in) :: T0

real*8, intent(out) :: Cp_total

real*8  :: T, Tmenos1 
real*8  :: Cp, Cpmenos1, Cp_componente
integer :: i, j

integer :: n

character*40 :: arquivo

logical existe



Cp_total = 0
do i=1, n_f_molares
    ! Encontra o nome da tabela
    write(arquivo,'(A,A)') trim(adjustl(n_frac_molar(i))), ".txt"
    
    ! Pergunta se a tabela existe
    inquire(file="tabelas\" // arquivo, exist = existe)
    
    ! Se não existir mostra na tela, e desconsidera este
    if(.not. existe) then
        write(4,*) " TABELA DE CALORES ESPECIFICOS INEXISTENTE PARA ESTA FRACAO: ", trim(adjustl(n_frac_molar(i))), " )       "
        goto 199
    
    endif
    
    ! Abre a tabela
    open(1,file= "tabelas\" // arquivo)
    read(1,*)
    read(1,*)
    
    
    !-----------------------------------------------
    ! Conta no número de linhas do arquivo
    read(1,*)
    do j=1, 1000
        read(1,*, iostat=n)
        if(n /= 0) exit
    enddo
    n = j
    close(1)
    open(1,file= "tabelas\" // arquivo)
    read(1,*)
    read(1,*)
    !-----------------------------------------------
    
    
    ! Pega o primeiro valor de Temperatura e calor específico para comparação
    read(1,*) Tmenos1, Cpmenos1
    
    do j = 1, n
    
        read(1,*) T, Cp
        

        if( T0 > Tmenos1 .and. T0 < T ) then
            
            Cp_componente = frac_molar(i)*( Cp + (Cpmenos1-Cp)*(T0-T)/(Tmenos1-T) ) 
            exit        
        else
            Tmenos1 = T
            Cpmenos1 = Cp
        
        endif
              
    
    enddo

    close(1)
    Cp_total = Cp_total + Cp_componente
    
    199 continue
enddo



end subroutine specific_heat
!---------------------------------------------------------------------------------
end program propep_2ph