Program empuxo1p4
    
! Diego
! versão 1.4:  10 Ago 15
! Colocado a eficiência da temperatura adiabática de chama, tabela dos parâmetros da taxa de queima
! Revisado o calculo no período de despressurização
! Tipo de propelente pode ser escolhido via arquivo de entrada    
    
! Diego
! versão 1.3:  17 Jun 15
! Nova equação é resolvida, baseada na massa de gases gerada apenas, novo método de resolução
! é colocado (bisseção modificado)

! Diego
! versão 1.2.2.1:  25 Jun 15
! Parâmetros da taxa de queima foram alterados, propriedades termodinâmicas são agora calculadas com base na pressão 
! de estagnação, grão tipo 6 foi colocado (tubular com coroas desinibidas) e resolvido erro no calculo do volume do
! grão queimado na câmara 
    
! Diego
! versão 1.2.2:  12 Fev 15
! Conservação da massa foi avaliada nos 4 tipos de grão e há a conservação da massa nos 3 tipos de grão

! Diego
! versão 1.2.1:  8 Dez 14
! Corrigido erro no cálculo do volume da câmara de combustão
! Arquivos de saída são colocados na pasta output, via título que é uma variável de entrada
! do programa

! Diego
! versão 1.2:  5 Dez 14
! Adição do tipo de grão 4 (grão de ativação constante) com animação em gif 
! , melhor interação com o usuário (dados via arquivo de texto)
! , adicionado caso onde a pressão gerada na câmara é inferior a pressão ambiente
! assim, a pressão se torna a pressão ambiente, não gerando portanto força alguma

! Marchi
! versão 1.1:  4 Jun 05
! versão 1.0: 30 Mai 05
! última alteração: 1 Nov 2013

USE PORTLIB

implicit none

real*8  :: z1, z2, z3
real*8  :: Vol, Rg, Tc, Ab, At, a, n, rop, gama, ac, cte, Vola, Vol0
real*8  :: Rga, Tca, gamaa, qui
real*8  :: L, t, dt, x, r, tol
real*8  :: Ri, W, fmt, Ve, Fd, Fmax, Cf, Fmed, Itot, Is, Ismed
real*8  :: pmed, pmax
real*8  :: p, pa, pao, paux, dp, Vtc, tq, z, g
real*8  :: Ai, Ae, Au, Voli, Vole, Volu, Re

real*8 :: p1, p2, p3
real*8 :: p1l, p2l, p3l

real*8 :: Dg, Di, Da, Dgrao, Mcri
real*8 :: r_ativ, t_inicio, mp

real*8, dimension(:), allocatable :: Rativ

real*8, parameter :: pi = dacos(-1.d0)

real*8 :: pedaco_grao, ra, temp
real*8 :: f, df
real*8 :: y1, y2, x1, x2, y, h
integer :: indice_grao, indice_grao1

integer :: passos, i, it, itmax, grao
integer :: ver, j

integer :: freq, passo
character*128 :: texto, arquivo
integer :: ativacao

real*8 :: relax
real*8 :: ho 

real*8 :: Cd, Et


integer :: prop, bates

integer :: tab_pq
integer :: n_tab, flag, itab
real*8, allocatable, dimension(:) :: P0tab, ntab, atab


character(128) :: entrada
t_inicio = 0.d0

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

!i= system("Notepad entrada.txt")
! dados
open(1,file=trim(adjustl(entrada)))
read(1,*) Arquivo
read(1,*) Dg
read(1,*) Di
read(1,*) Da
read(1,*) L
read(1,*) grao
read(1,*) freq
read(1,*) r_ativ
read(1,*) rop
read(1,*) a
read(1,*) n
read(1,*) pao
read(1,*) dt
read(1,*) passos ! Numero máximo de passos no tempo
read(1,*) itmax ! Numero máximo de iterações para o ponto do calculado
read(1,*) tol ! Tolerância para cada ponto calculado
read(1,*) relax ! Fator de relaxamento
read(1,*) ho ! Passso inicial para o processo iterativo da bisseção
read(1,*) cd ! Coeficiente de descarga
read(1,*) Et ! Eficiencia da To
read(1,*) tab_pq ! Usar ou não a tabela de parâmetros da taxa de queima
read(1,*) prop

! Se for usar os parâmetros da taxa de queima via tabela:
if(tab_pq == 1)then
    open(1,file='tabela_pq.txt')
    read(1,*)
    do i=1, 10000
        read(1,*, iostat=flag)
    
        if(flag /= 0) exit
    enddo
    
    close(1)
    n_tab = i-1
    
    allocate(P0tab(n_tab), ntab(n_tab), atab(n_tab))
    open(1,file='tabela_pq.txt')
    read(1,*)
    do i=1, n_tab
        read(1,*) P0tab(i), ntab(i), atab(i)
    enddo
    
    close(1)
    
endif


! Volume inicial da câmara de combustão (m3)
Vol0 =    0.25d0*pi*Da**2*L

! Área de queima para o grão 1 (m2)
Ab   =    0.25d0*pi*Di**2  


! Área da garganta (m2)
At   =    0.25d0*pi*Dg**2  

! Volume total da câmara de combustão
Vtc  =    pi*0.25d0*Di**2*L

! Raio interno do grão (m)
Ri   =    0.5d0*Da

! Raio externo do grão (m)
Re   =    0.5d0*Di

x      = 0.0d0
pa     = pao

! Aceleração gravitacional (m/s2)
g   = 9.80665d0

! L
! 1-Comprimento do grão (m)
! 2-Espessura do grão (m)
! 3-Metade da espessura do grão (m)
select case (grao)
case (2)
    W = L
    L = Re-Ri
    Ab = 2*pi*W*(Ri + x)
case(3)
    W = L
    L = (Di-Da)*0.250d0
case(4)

    ! Calculo do numero de discretizações para a ativação do grão-propelente
    ativacao = int(L/r_ativ/dt)+1

    ! Discretização do raio do grão-propelente em 'ativacao' partes
    allocate(Rativ(ativacao))
    
    Rativ = 0.d0
    
    ! Dgrao = comprimento de cada pedaço de grão discretizado (m)
    Dgrao = L/ativacao
    ! Diferencas finitas!!
    ! i =1, x = 0
    ! i = ativacao+1, x= L
    write(*,*) "delta_grao = ", Dgrao
    write(*,*) "ativacao = ", ativacao
    
    Rativ = Ri


endselect
    

open(8,file="./output/" // trim(adjustl(arquivo)) // ".txt")

write(8,10)
10 format("#", t8,"i", t15,"tempo (s)", t30,"Fd (N)",   t45,"r (mm/s)", t60,"p (bar)", &
                  t75,"K (adim.)", t90,"fmt (kg/s)", t105,"Ve (m/s)",              &
                 t120,"x (mm)",   t135,"Cf (adim.)", t150,"Is (s)",   &
                 t165,"dp (Pa)",   t175,"iterações",   t190,"a",   t205,"n")

! instante inicial

p   = pao
Vol = Vol0

if(tab_pq == 1) then
    ! Escolhe os parametros baseado na pressão
    do itab=1, n_tab
        if(p <= P0tab(itab) ) then 
            exit   
        endif
    enddo   
    if(itab>n_tab) itab=n_tab
    a = atab(itab)
    n = ntab(itab)
endif
r   = a*(p**n)
t   = 0.0d0
fmt = 0.0d0
Ve  = 0.0d0
Fd  = 0.0d0
Fmax=0
Cf  = 0.0d0
Is  = 0.0d0
dp  = 0.0d0
Itot= 0.0d0
Fmed= 0.0d0
pmed= 0.0d0 
pmax= 0.0d0
it  = 0
i   = 0
temp= 0
passo=0
t_inicio= 0.0d0 
Ismed = 0
!write(8,20) i, t, Fd, r*1000, p/101325.0d0, Ab/At, fmt, Ve, x*1000, Cf, Is, dp, it
20 format (i8, 11(1pe15.4), i8, 10(1pE15.4))



indice_grao = 0

write(*,*) " ::: INCIANDO CALCULOS :::"

do i = 1, passos



   Vola = Vol

   p  = pa
   ra = r
   
    if(tab_pq == 1) then
        ! Escolhe os parametros baseado na pressão
        do itab=1, n_tab
            if(p <= P0tab(itab) ) then 
                exit   
            endif
        enddo   
        if(itab>n_tab) itab=n_tab
        a = atab(itab)
        n = ntab(itab)
    endif
   r  = a*(p**n)

   x  = x + r*dt

   t  = i * dt
   
   
   select case (grao)
         
   case(1)
       
      Vol = Vol0 + Ab*x

    case(2)
      Ab  = 2*pi*W*(Ri + x)
      Vol = Vol0 + pi*W*( ((Ri+x)**2) - (Ri**2) )
      
    case(3)
      Ai  = 2 * pi * (Ri+x) * (W-2*x)
      Ae  = 2 * pi * (Re-x) * (W-2*x)
      Au  = 2 * pi * ( ((Re-x)**2) - ((Ri+x)**2) )
      Ab  = Ai + Ae + Au
      Voli= pi * W * ( (Ri+x)**2 )
      Vole= pi * W * ( (Re**2) - ((Re-x)**2) )
      Volu= pi * 2 * x * ( ((Re-x)**2) - ((Ri+x)**2) )
      Vol = Vol0 + Voli + Vole + Volu

    case (4)
        indice_grao1 = indice_grao
        if(i <= ativacao) then
            Rativ(1:i) = Rativ(1:i) + r*dt
        else
            Rativ(indice_grao+1:ativacao) = Rativ(indice_grao+1:ativacao) + r*dt
        endif
        do j=1, ativacao           
            ! Se for além do máximo permitido (Diâmetro externo do grão) colocar o valor de Re no Rativ
            if( Rativ(j) > Re) then 

                indice_grao = j
                Rativ(j) = Re
            endif
        enddo
        
        ! Calculo das áreas de queima e do volume atual
        Ab = 0.d0
        
        if(i > ativacao) then
            do j=indice_grao+2, ativacao
                    Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j))
            enddo
        else
            if(i==1) Ab = Ab + 0.5d0*Rativ(j-1)
            do j=indice_grao+2, i
                    Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j))
            enddo
        endif
        
        if(indice_grao - indice_grao1 > 0) Ab = Ab + (indice_grao - indice_grao1)*Re
        
        Vol= 0.d0
        do j=2, ativacao 
                Vol= Vol+ ( 0.5d0*(Rativ(j-1)+Rativ(j)) )**2 
        enddo
!        write(4,'(i5,i5, 19(1PE15.5))') indice_grao1, indice_grao - indice_grao1,&
!        Rativ(indice_grao), Ab*2*pi*dgrao,  Vol*pi*dgrao
        
        Ab  = Ab*2*pi*dgrao 
        Vol = Vol*pi*dgrao
        
!        write(*,'(I10,3(1PE15.3))') i, Ab, vol, Vtc

    if (Ab <= 1E-12) exit
        
    if(freq > 0) then
        if(mod(i,freq) == 0) then
        
        

        passo = passo + 1
        if(passo<10) then
            write(texto,'("0000",I1)') passo
        elseif(passo<100) then
            write(texto,'("000",I2)') passo
        elseif(passo<1000) then
            write(texto,'("00",I3)') passo
        elseif(passo<10000) then
            write(texto,'("0",I4)') passo
        elseif(passo<100000) then
            write(texto,'(I5)') passo
        endif
        
        OPEN(2,FILE='./temp/'// trim(adjustl(texto)) //'.txt')
        DO J=1, ativacao
                WRITE(2,*) dgrao*(j-1), Rativ(J)
        ENDDO
        WRITE(2,*) dgrao*(j-1), Rativ(J-1)
        close(2)
        endif
    endif
    
    ! Queima cilíndrica com parte exterior (externa+coroas) desinibidas
    case(5)
      ! Área das coroas
      Au =  2 * pi * (Re-x) ** 2  
      ! Área externa
      Ae  = 2 * pi * (Re-x) * (L-2*x)
      Ab  = Ae + Au
      ! Volume das coroas 
      Volu= pi* 2 * x  * ((Re-x)**2)
      ! Volume externo
      Vole= pi * ( Re**2 - (Re-x)**2)*L
      Vol = Vol0 + Vole + Volu
      
      
    case(6) ! Queima tubular com  2 coroas desinibidas
     ! Área interna        
      Ai  = 2 * pi * (Ri+x) * (L-2*x)
      ! Área das coroas
      Au  = 2 * pi * ( ((Re)**2) - ((Ri+x)**2) )
      ! Área de queima total
      Ab  = Ai + Au
      ! Volume interno
      Voli= pi * (L-2*x) * ( (Ri+x)**2 - Ri**2 )
      ! Volume das coroas 
      Volu= pi * 2 * x * ( ((Ri+x)**2) - (Ri**2))
      
      Vol = Vol0 + Voli + Volu
      
    case(7) ! Bates com j subdivisão
        Ai = 0
        Au = 0
        Voli = 0
        Volu = 0
        bates = 4
       do j = 1, bates
        ! Área interna da parte j
          Ai  = Ai + 2 * pi * (Ri+x) * (L/bates-2*x)
      
          ! Área das coroas da parte j
          Au  = Au + 2 * pi * ( Re**2 - (Ri+x)**2 )

          ! Volume interno da parte j
          Voli= Voli + pi * (L/bates-2*x) * ( (Ri+x)**2 - Ri**2 )

          ! Volume das coroas da parte j
          Volu= Volu + pi * 2 * x * ( (Ri+x)**2 - Ri**2)
      
    enddo
    
      ! Área de queima total
      Ab  = Ai + Au
      Vol = Vol0 + Voli + Volu
      
           
      
    case(8) !  j subdivisão
        Ai = 0
        Au = 0
        Voli = 0
        Volu = 0
        bates = 2
       do j = 1, bates
        ! Área interna da parte j
          Ai  = Ai + 2 * pi * (Ri+x) * (L/j-2*x)
      
          ! Área das coroas da parte j
          Au  = Au + 2 * pi * ( ((Re)**2) - ((Ri+x)**2) )

          ! Volume interno da parte j
          Voli= Voli + pi * (L/j-2*x) * ( (Ri+x)**2 - Ri**2 )

          ! Volume das coroas da parte j
          Volu= Volu + pi * 2 * x * ( ((Ri+x)**2) - (Ri**2))
      
    enddo
    
      ! Área de queima total
      Ab  = Ai + Au
      Vol = Vol0 + Voli + Volu
      
           
    end select

    ! RESOLUCAO NUMERICA DA EQ DO TRANSIENTE DA PRESSURIZACAO DO MOTOR-FOGUETE
    ! VIA METODO DA SECANTE MODIFICADO
   h = ho
   
       
    ! ------------------------------------------------------------
    ! INICIO DO calculo para encontrar o ponto 1
    ! ------------------------------------------------------------
    call Fk(prop,p,gama)
    call FR(prop,p,Rg)
    call FT0(prop,p,Tc)
    Tc = Tc*Et
    call FQui(prop,p, Qui)
    

    ! Calcula Mcri
    call M_cri(p, pao, gama, Mcri)
    
    ! Escoamento totalmente subsônico
    if(Mcri< 1) then
        p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc !Ab*a*p**n*(rop*(1-qui)-p/Rg/Tc)
        p2 = Cd*Mcri*At*p*sqrt(gama/Rg/Tc*(1+0.5d0*(gama-1)*Mcri**2)**(-(gama+1)/(gama-1)))
        p3 = (Vol*p - Vola*pa)/Rg/Tc
        
        f = (p1-p2)*dt-p3
        
        
    ! Escoamento totalmente supersônico
    elseif(Mcri >= 1) then
        p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc
        p2 = Cd*At*p*sqrt(gama/Rg/Tc*((1+gama)/2)**(-(gama+1)/(gama-1)))
        p3 = (Vol*p - Vola*pa)/Rg/Tc
        
        f = (p1-p2)*dt-p3

        
    endif
    ! ------------------------------------------------------------
    ! Atribuição do ponto 1
    ! ------------------------------------------------------------
    y1 = f
    x1 = p
   do it = 1, itmax

        paux = p

        
        ! ------------------------------------------------------------
        ! Calcula o valor do próximo ponto de pressão a ser calculado
        ! ------------------------------------------------------------
        
        p = p + h
        
        
        ! ------------------------------------------------------------
        ! INICIO DO calculo para encontrar o ponto 2
        ! ------------------------------------------------------------
        call Fk(prop,p,gama)
        call FR(prop,p,Rg)
        call FT0(prop,p,Tc)
        Tc = Tc*Et
        call FQui(prop,p, Qui)
        
    
        ! Calcula Mcri
        call M_cri(p, pao, gama, Mcri)
        
        if(Mcri /= Mcri .or. Mcri == 0) Mcri = 1d-1
        ! Escoamento totalmente subsônico
        if(Mcri< 1) then
            p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc !Ab*a*p**n*(rop*(1-qui)-p/Rg/Tc)
            p2 = Cd* (Mcri*At*p*sqrt(gama/Rg/Tc*(1+0.5d0*(gama-1)*Mcri**2)**(-(gama+1)/(gama-1))))
            p3 = (Vol*p - Vola*pa)/Rg/Tc
            
            f = (p1-p2)*dt-p3

            
        ! Escoamento totalmente supersônico
        elseif(Mcri >= 1) then
            p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc
            p2 = Cd* ( At*p*sqrt(gama/Rg/Tc*((1+gama)/2)**(-(gama+1)/(gama-1))) )
            p3 = (Vol*p - Vola*pa)/Rg/Tc
            
            f = (p1-p2)*dt-p3

            
        endif
        ! ------------------------------------------------------------
        ! Atribuição do ponto 2
        ! ------------------------------------------------------------
        y2 = f
        x2 = p
        
        
        ! ------------------------------------------------------------
        ! Novo valor a ser calculado para chegar ao zero da função y com interpolação linear
        ! ------------------------------------------------------------
        h = (x2-x1)/(y2-y1)*(0.d0-y1)*relax
        
!        write(*,*) it, p, h, dp
!        pause
        y1 = y2
        x1 = x2

          dp = dabs(p-paux)
          
            

          ! Limites
          if(p <= pao) then
                p = (pao + abs(p))

          elseif(p>= 1e10) then
                p = 0.5d0*(1e10 + abs(p))
                
          endif
          

            if( p /= p) then
                p = pao
                exit
                endif
            if(dp < tol) exit
    

   end do
    ! SE o processo iterativo não convergir, coloca-se a pressão ambiente no  valor da pressão
      if(it == itmax+1) then
        p = pao
      endif 
   
    call Fk(prop,p,gama)
    call FR(prop,p,Rg)
    call FT0(prop,p,Tc)
    Tc = Tc*Et
    call FQui(prop,p, Qui)
        
    
    ! Calcula Mcri
    call M_cri(p, pao, gama, Mcri)
    
    ! Calcula fluxo de massa considerando as possibilidades
    ! Escoamento totalmente subsônico
    If(Mcri<1) then
        fmt = Cd*Mcri*At*p*sqrt(gama/Rg/Tc*(1+0.5d0*(gama-1)*Mcri**2)**(-(gama+1)/(gama-1)))
    ! Escoamento totalmente supersônico
    elseif(Mcri>=1) then
        fmt = Cd*At*p*sqrt(gama/Rg/Tc*((1+gama)/2)**(-(gama+1)/(gama-1)))
    endif
    
   ve = Mcri*sqrt(gama*Rg*Tc*(1+(gama-1)/2*Mcri**2)**(-1))
          
   Fd  = fmt*Ve

   Itot= Itot + Fd*dt
   Is  = Fd/(fmt*g)
   
   ! Se o fluxo de massa der zero é posível que dê problemas no impulso específico
   if(Is /= Is) Is = 0
   Cf  = Fd/(p*At)
   Ismed = Ismed + Is*dt
   
   temp = temp + fmt
   
   ! Esquema para encontrar a força de empuxo máxima
   if(Fd>=Fmax) Fmax = Fd
   
    ! Esquema para encontrar a pressão máxima
   if(p>=pmax) pmax = p
   
    if(p>pao) pmed = pmed + p*dt


   write(8,20) i, t, Fd, r*1000, p/1d5, Ab/At, fmt*1e3, Ve, x*1000, Cf, Is, dp, it, a, n

    select case (grao)
    case default
        if ( x > L ) exit
    case(5)
        if( Ab < tol) exit
    case(6)
        if( x >= (Re-Ri)) then
            write(*,*) "x > (Re-Ri)"
            exit
        elseif(x>(L/2)) then
            write(*,*) "x > L/2"
            exit
        endif
    case(7)        
        if(x >= (Re-Ri) .or. x > L/bates/2 ) exit
   end select

   pa = p

end do
! despressurização


write(*,*) " ::: DESPRESSURIZOU :::"

if(grao /= 4) then
    
    Ab = 0.0d0
    r  = 0.0d0
    dp = 0.0d0
    it = 0
    tq = t
    vol = vola

    do j = i+1, passos

       pa = p
       
        if(tab_pq == 1) then
            ! Escolhe os parametros baseado na pressão
            do itab=1, n_tab
                if(p <= P0tab(itab) ) then 
                    exit   
                endif
            enddo   
            if(itab>n_tab) itab=n_tab
            a = atab(itab)
            n = ntab(itab)
        endif
       t  = j * dt

    ! RESOLUCAO NUMERICA DA EQ DO TRANSIENTE DA PRESSURIZACAO DO MOTOR-FOGUETE
    ! VIA METODO DA BISSEÇÂO MODIFICADO
   h = ho
   
    ! ------------------------------------------------------------
    ! INICIO DO calculo para encontrar o ponto 1
    ! ------------------------------------------------------------
    call Fk(prop,p,gama)
    call FR(prop,p,Rg)
    call FT0(prop,p,Tc)
    Tc = Tc*Et
    call FQui(prop,p, Qui)
    

    ! Calcula Mcri
    call M_cri(p, pao, gama, Mcri)
    
    ! Escoamento totalmente subsônico
    if(Mcri< 1) then
        p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc !Ab*a*p**n*(rop*(1-qui)-p/Rg/Tc)
        p2 = Cd*Mcri*At*p*sqrt(gama/Rg/Tc*(1+0.5d0*(gama-1)*Mcri**2)**(-(gama+1)/(gama-1)))
        p3 = (Vol*p - Vola*pa)/Rg/Tc
        
        f = (p1-p2)*dt-p3
        
        
    ! Escoamento totalmente supersônico
    elseif(Mcri >= 1) then
        p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc
        p2 = Cd*At*p*sqrt(gama/Rg/Tc*((1+gama)/2)**(-(gama+1)/(gama-1)))
        p3 = (Vol*p - Vola*pa)/Rg/Tc
        
        f = (p1-p2)*dt-p3

        
    endif
    ! ------------------------------------------------------------
    ! Atribuição do ponto 1
    ! ------------------------------------------------------------
    y1 = f
    x1 = p
       
   do it = 1, itmax

        paux = p

        
        ! ------------------------------------------------------------
        ! Calcula o valor do próximo ponto de pressão a ser calculado
        ! ------------------------------------------------------------
        
        p = p + h
        
        
        ! ------------------------------------------------------------
        ! INICIO DO calculo para encontrar o ponto 2
        ! ------------------------------------------------------------
        call Fk(prop,p,gama)
        call FR(prop,p,Rg)
        call FT0(prop,p,Tc)
        Tc = Tc*Et
        call FQui(prop,p, Qui)
        
    
        ! Calcula Mcri
        call M_cri(p, pao, gama, Mcri)
        
        if(Mcri /= Mcri .or. Mcri == 0) Mcri = 1d-1
        ! Escoamento totalmente subsônico
        if(Mcri< 1) then
            p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc 
            p2 = Cd*Mcri*At*p*sqrt(gama/Rg/Tc*(1+0.5d0*(gama-1)*Mcri**2)**(-(gama+1)/(gama-1)))
            p3 = (Vol*p - Vola*pa)/Rg/Tc
            
            f = (p1-p2)*dt-p3

            
        ! Escoamento totalmente supersônico
        elseif(Mcri >= 1) then
            p1 = Ab*a*p**n*rop*(1-qui) - Ab*a*p**(n+1)/Rg/Tc
            p2 = Cd*At*p*sqrt(gama/Rg/Tc*((1+gama)/2)**(-(gama+1)/(gama-1)))
            p3 = (Vol*p - Vola*pa)/Rg/Tc
            
            f = (p1-p2)*dt-p3

            
        endif
        ! ------------------------------------------------------------
        ! Atribuição do ponto 2
        ! ------------------------------------------------------------
        y2 = f
        x2 = p
        
        
        ! ------------------------------------------------------------
        ! Novo valor a ser calculado para chegar ao zero da função y com interpolação linear
        ! ------------------------------------------------------------
        h = (x2-x1)/(y2-y1)*(0.d0-y1)*relax
        
!        write(*,*) it, p, h, dp
!        pause
        y1 = y2
        x1 = x2


          dp = dabs(p-paux)
          
            

          ! Limites
          if(p <= pao) then
                p = (pao + abs(p))

          elseif(p>= 1e10) then
                p = 0.5d0*(1e10 + abs(p))
                
          endif
          

            if( p /= p) then
                p = pao
                exit
                endif
            if(dp < tol) exit
    

       end do 
   
        call Fk(prop,p,gama)
        call FR(prop,p,Rg)
        call FT0(prop,p,Tc)
        Tc = Tc*Et
        call FQui(prop,p, Qui)
        
    
        ! Calcula Mcri
        call M_cri(p, pao, gama, Mcri)
    
        If(Mcri<1) then
            fmt = Cd*Mcri*At*p*sqrt(gama/Rg/Tc*(1+0.5d0*(gama-1)*Mcri**2)**(-(gama+1)/(gama-1)))
        elseif(Mcri>=1) then
            fmt = Cd*At*p*sqrt(gama/Rg/Tc*((1+gama)/2)**(-(gama+1)/(gama-1)))
        endif
    
       ve = Mcri*sqrt(gama*Rg*Tc*(1+(gama-1)/2*Mcri**2)**(-1))
          
       Fd  = fmt*Ve

       Itot= Itot + Fd*dt
       Is  = Fd/(fmt*g)
       Cf  = Fd/(p*At)
       Ismed = Ismed + Is*dt

       temp = temp + fmt


       write(8,20) j, t, Fd, r*1000, p/1.0d5, Ab/At, fmt*1e3, Ve, x*1000, Cf, Is, dp, it, a, n

        if ( Fd <= 0.3d0 ) exit
    end do

endif

close(8)

tq = t

pmed = pmed / tq

Fmed  = Itot / tq


Ismed = Ismed / tq


write(*,'(//a)') " ::: FIM DOS CALCULOS :::"
open(8,file="./output/" // trim(adjustl(arquivo)) // "_resultados.txt")

87 FORMAT (I10, A)
88 FORMAT (1PE10.3, A)
   write(8,'(/A/)') "::: PARAMETROS DE ENTRADA DA SIMULAÇÃO :::"
write(8,'(A)') "'" // trim(adjustl(Arquivo)) // "' :: Titulo (Até 128 caracteres)"
write(8,88) Dg, " :: Diâmetro da garganta (m)"
write(8,88) Di, " :: Diâmetro externo do grão (m)"
write(8,88) Da, " :: Diâmetro da alma do grão (m)"
write(8,88) L,  " :: Comprimento do grão (m)"
write(8,87) grao, " :: grao (1 = cilíndrico (queima cigarro), 2 = tubular (queima interna), 3 = tubular (queima em toda superfície), 4 = taxa de ativação constante (de fora pra dentro), 5 = Cilindrico com superfície externa (externa+coroas) desinibidas, 6 = Tubular com coroas desinibidas)"
write(8,87) freq, " :: P/ grao tipo 4: Frequencia para criacao da animacao gif (freq = 0 para desabilitar)"
write(8,88) r_ativ, " :: P/ grao tipo 4: Taxa de ativação do grão-propelente (m/s)"
write(8,88) rop, " :: Massa específica do propelente antes da queima(kg/m3)"
write(8,88) a, " :: a = m/s/Pa, Parâmetros da taxa de queima do propelente r=a*p^n" 
write(8,88) n, " :: n = adim, Parâmetros da taxa de queima do propelente r=a*p^n"
write(8,88) pao, " :: Pressão atmosférica (Pa)"
write(8,88) dt, " :: Passo de tempo (s)"
write(8,87) passos, " :: Numero máximo de passos no tempo"
write(8,87) itmax, " :: Numero máximo de iterações para o ponto do calculado"
write(8,88) tol, " :: Tolerância para cada ponto calculado"
write(8,88) relax, " :: relax  = Fator de relaxamento do metodo da bisseção"
write(8,88) ho, " :: ho  = Passo inicial do chute "
write(8,88) cd, " :: Cd = Coeficiente de Descarga"
write(8,88) Et, " :: Et = Eficiência da temperatura adiabática de chama"
write(8,87) tab_pq, " :: tab_pq = 0: Não usa tabela dos paràmetros da taxa de queima, 1: usa"
write(8,87) prop," :: prop: 0=KNSu 65/35, 1=Polvora Bandeirante (150/50/20), 2=KNSu 60/40, 3=KNSu 70/30"

write(8,'(/A/)') "::: RESULTADOS DA SIMULAÇÃO :::"

    select case (grao)
    case(1)
            mp = 0.25d0*pi*(Di**2)*L*rop
            write(8,'(10(//, 4x, 1PE22.15, A))') temp*dt, " = Massa de gás ejetada ", &
            mp, " = Massa total embarcada de propelente no motor "    
    case (2,3)
            mp = 0.25d0*pi*(Di**2-Da**2)*W*rop
            write(8,'(10(//, 4x, 1PE22.15, A))') temp*dt, " = Massa de gás ejetada ", &
            mp, " = Massa total embarcada de propelente no motor "
    case (4,6)
            mp = 0.25d0*pi*(Di**2-Da**2)*L*rop
            write(8,'(10(//, 4x, 1PE22.15, A))') temp*dt, " = Massa de gás ejetada ", &
            mp, " = Massa total embarcada de propelente no motor "    
    case (5)
            mp = 0.25d0*pi*Di**2*L*rop
            write(8,'(10(//, 4x, 1PE22.15, A))') temp*dt, " = Massa de gás ejetada ", &
            mp, " = Massa total embarcada de propelente no motor "    
    end select

call classe(Itot, Fmed, texto)

write(8,'(/A,A)') "CLASSE: ", trim(adjustl(texto))

write(8,30) tq, t-tq, Fmed, Fmax, Itot, Pmed/1d5, Pmax/1d5, Itot/temp/dt/g, Itot/temp/dt, Itot/mp/g, Itot/mp
30 format(/, 5x, "Tempo de queima total (s)     = ", F15.4, &
         /, 5x, "Tempo de despressurização (s) = ", F15.4, &
         /, 5x, "Empuxo médio (N)              = ", F15.4, &
         /, 5x, "Empuxo máximo (N)             = ", F15.4, &
         /, 5x, "Impulso total (N.s)           = ", F15.4, &
         /, 5x, "Pressão média (bar)           = ", F15.4, &
         /, 5x, "Pressão máxima (bar)          = ", F15.4, &
         /, 5x, "Impulso específico (s)        = ", F15.4, &
         /, 5x, "Velocidade de ejeção efetiva média (m/s)  = ", F15.4, &
         //, 5x,":: Parâmetros baseados na massa de propelente embarcada (real) :: ", &
         /, 5x, "Impulso específico (s)        = ", F15.4, &
         /, 5x, "Velocidade de ejeção efetiva média (m/s)  = ", F15.4  )


close(8)


! GRAFICO DA FORÇA
open(2,file="./output/" // trim(adjustl(arquivo)) // '_forca.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "_forca.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'F (N)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:3 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "_forca.gnu")

! GRAFICO DA PRESSAO
open(2,file="./output/" // trim(adjustl(arquivo)) // '_pressao.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "_pressao.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'pressão (bar)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:5 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "_pressao.gnu")

! GRAFICO DA RAZAO DE QUEIMA (Kn)
open(2,file="./output/" // trim(adjustl(arquivo)) // '_Kn.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "_Kn.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'Kn (adim)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:6 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // '_Kn.gnu')

! GRAFICO Do fluxo de massa
open(2,file="./output/" // trim(adjustl(arquivo)) // '_fm.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "fm.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'Fluxo de Massa (g/s)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:7 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // '_fm.gnu')

! GRAFICO Da velocidade de ejeção
open(2,file="./output/" // trim(adjustl(arquivo)) // '_ve.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "ve.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'Velocidade de Ejeção (m/s)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:8 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "_ve.gnu")


if(freq>0 .and. grao==4) then
    OPEN(2,FILE='./temp/grao.txt')
    write(2,*) 0.d0, Ri
    write(2,*) L, Ri
    write(2,*) L, Re
    write(2,*) 0.d0, Re
    write(2,*) 0.d0, Ri
    close(2)

    open(2,file='./output/' // trim(adjustl(arquivo)) // 'plot_gif.gnu')

    write(2,*) "set terminal gif animate delay 10 size 1280, 720"
    write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "progresso.gif'"
    write(2,*) "set style data lines"
    write(2,*) "set grid"
    write(2,*) "set xlabel 'x (m)'"
    write(2,*) "set ylabel 'r (m)'"
    write(2,'(A,F13.4,A,F13.4,A)') "set xrange [", -0.001d0, ":", L+.001d0, "]"
    write(2,'(A,F13.4,A,F13.4,A)') "set yrange [", -0.001d0, ":", Re+.001d0, "]"


    do i=1, passo
        if(i<10) then
            write(texto,'("0000",I1)') i
        elseif(i<100) then
            write(texto,'("000",I2)') i
        elseif(i<1000) then
            write(texto,'("00",I3)') i
        elseif(i<10000) then
            write(texto,'("0",I4)') i
        elseif(i<100000) then
            write(texto,'(I5)') i
        endif
        
        write(2,*) "plot './temp/"// trim(adjustl(texto)) //".txt' using 1:2 title 'Frente de chama' \"
        write(2,*) ", './temp/grao.txt' using 1:2 title 'Grão-propelente' lw 2 "

    enddo
    close(2)
    i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "plot_gif.gnu")
endif

!ver = system("Notepad saida2.txt")

if(grao == 4) ver = system("del .\temp\*.txt")

!ver = system("start Notepad ./output/" // trim(adjustl(arquivo)) // "_resultados.txt")

contains
!------------------------------------------------------------------------------------------
subroutine FQui(prop,P0, Qui)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: Qui
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0
    


select case (prop)
    ! KNSU 65/35
    case (0)
A0 =-2.04799653852393E+01
A1 = 3.62113343656072E+01
A2 =-2.77291460981005E+01
A3 = 1.21466806224666E+01
A4 =-3.34936868681469E+00
A5 = 6.02653493751515E-01
A6 =-7.07220352937977E-02
A7 = 5.21876383009207E-03
A8 =-2.19830376383906E-04
A9 = 4.03090411144927E-06
    ! Polvora Roberto 150/50/20
    case (1)
A0 = 7.26043568864702E+01
A1 =-1.44059368611776E+02
A2 = 1.26207419659459E+02
A3 =-6.37177913512304E+01
A4 = 2.04215720064625E+01
A5 =-4.30677147952935E+00
A6 = 5.97321868085242E-01
A7 =-5.25030031896640E-02
A8 = 2.65184174254155E-03
A9 =-5.85854750762796E-05

    ! Polvora KNSU 60/40
    case (2)
A0 = 1.62460275168875E+01
A1 =-2.76930932627632E+01
A2 = 2.07639961634161E+01
A3 =-8.85836008568855E+00
A4 = 2.37098159313768E+00
A5 =-4.12691242859834E-01
A6 = 4.67339057507249E-02
A7 =-3.32430078032562E-03
A8 = 1.35010149479413E-04
A9 =-2.38958104337635E-06
 
    ! Polvora KNSU 70/30
    case (3)
A0 = 1.99837170806653E+00
A1 =-4.70322436642678E+00
A2 = 4.66891945154650E+00
A3 =-2.54946842877398E+00
A4 = 8.56376845398039E-01
A5 =-1.84348060636415E-01
A6 = 2.55205817914470E-02
A7 =-2.19476961416862E-03
A8 = 1.06509629421422E-04
A9 =-2.22494962331610E-06

end select

logP0 = dlog10(P0)
! 0.01 bar a 11 kbar
Qui = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0

end subroutine Fqui
!------------------------------------------------------------------------------------------
subroutine FR(prop,P0, R)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: R
    
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0
    

select case (prop)
    ! KNSU 65/35
    case (0)
A0 = 3.93261975545070E+02
A1 =-2.78261444164441E+02
A2 = 2.21731823032611E+02
A3 =-1.00207260435184E+02
A4 = 2.74321375819001E+01
A5 =-4.67641633470184E+00
A6 = 4.86891326745650E-01
A7 =-2.87651336070012E-02
A8 = 8.00061878913493E-04
A9 =-5.32904201135432E-06
    ! Polvora Roberto 150/50/20
    case (1)
A0 =-4.61816631501257E+04
A1 = 9.26089352007268E+04
A2 =-8.12663522815713E+04
A3 = 4.11056870481495E+04
A4 =-1.32027764825664E+04
A5 = 2.79128507736463E+03
A6 =-3.88244760033038E+02
A7 = 3.42396307890670E+01
A8 =-1.73615214889782E+00
A9 = 3.85331201552525E-02
        
    ! Polvora KNSU 60/40
    case (2)
A0 = 6.71732634788817E+03
A1 =-1.13880372797209E+04
A2 = 8.76949555557240E+03
A3 =-3.87374827622120E+03
A4 = 1.08094863188746E+03
A5 =-1.97577175720993E+02
A6 = 2.36492984951256E+01
A7 =-1.78690783631224E+00
A8 = 7.73150298950068E-02
A9 =-1.45940036826555E-03

    ! Polvora KNSU 70/30
    case (3)
A0 = 2.47645890800131E+02
A1 = 1.16173534468038E+02
A2 =-2.08988520875701E+02
A3 = 1.46939752486607E+02
A4 =-5.73358292473525E+01
A5 = 1.36290650072157E+01
A6 =-2.02380940138111E+00
A7 = 1.83269113930935E-01
A8 =-9.24730277758449E-03
A9 = 1.99039233116305E-04

end select

logP0 = dlog10(P0)
! 0.01 bar a 11 kbar
R = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0

end subroutine FR
!------------------------------------------------------------------------------------------
subroutine FT0(prop,P0, T0)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: T0
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0
    
select case (prop)
    ! KNSU 65/35
    case (0)
A0 =-1.37888992619268E+04
A1 = 2.63269770368983E+04
A2 =-2.03742329478610E+04
A3 = 9.02448771642549E+03
A4 =-2.50843498438118E+03
A5 = 4.53324034236469E+02
A6 =-5.31925781292568E+01
A7 = 3.90421818333940E+00
A8 =-1.62672493431325E-01
A9 = 2.93507026655967E-03
    ! Polvora Roberto 150/50/20
    case (1)
A0 =-1.81135947499828E+05
A1 = 3.63202971159722E+05
A2 =-3.18593433119871E+05
A3 = 1.61107084078521E+05
A4 =-5.17303063313623E+04
A5 = 1.09329676609606E+04
A6 =-1.52011083262037E+03
A7 = 1.34001521365648E+02
A8 =-6.79117562963686E+00
A9 = 1.50634004319533E-01

    ! Polvora KNSU 60/40
    case (2)
A0 =-3.93833638581385E+04
A1 = 7.15176408651588E+04
A2 =-5.54248331358298E+04
A3 = 2.46941412204016E+04
A4 =-6.96071798591121E+03
A5 = 1.28700072060828E+03
A6 =-1.55984837441572E+02
A7 = 1.19385199904895E+01
A8 =-5.23114870016792E-01
A9 = 9.99299233223467E-03
        
    ! Polvora KNSU 70/30
    case (3)
A0 =-2.47644986635254E+03
A1 = 4.22534690529101E+03
A2 =-1.60315577906195E+03
A3 =-4.80970623690143E+01
A4 = 2.41319545777364E+02
A5 =-8.83748663032046E+01
A6 = 1.61431174239937E+01
A7 =-1.65434932347325E+00
A8 = 9.06303622630090E-02
A9 =-2.06798110867363E-03

end select

logP0 = dlog10(P0)
! 0.01 bar a 11 kbar
T0 = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0

end subroutine FT0
!------------------------------------------------------------------------------------------
subroutine Fk(prop,P0, k)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: k
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0


select case (prop)
    ! KNSU 65/35
    case (0)
! kgas
!A0=-4.15198741081744E+00
!A1= 9.43179477126738E+00
!A2=-7.17478856172397E+00
!A3= 3.12235584269792E+00
!A4=-8.57362675678530E-01
!A5= 1.54068099574775E-01
!A6=-1.81219194658271E-02
!A7= 1.34583650541077E-03
!A8=-5.72884438851859E-05
!A9= 1.06539054876278E-06

! kmix
!A0= 1.89225403506303E+00
!A1=-1.23412794184640E+00
!A2= 9.73970487442510E-01
!A3=-4.41653674918613E-01
!A4= 1.25118018390715E-01
!A5=-2.29286320619050E-02
!A6= 2.71256388130721E-03
!A7=-1.99431987464383E-04
!A8= 8.26545014672395E-06
!A9=-1.47293882113353E-07

! k2ph
A0= 3.09510840624045E+00
A1=-3.47788297095420E+00
A2= 2.68589018179646E+00
A3=-1.19233536051624E+00
A4= 3.33445002556658E-01
A5=-6.08847172457605E-02
A6= 7.25433618821334E-03
A7=-5.43620288790910E-04
A8= 2.32502840539917E-05
A9=-4.32683123306007E-07
    ! Polvora Roberto 150/50/20
    case (1)
A0= 5.78092194482495E+00
A1=-9.31175283694160E+00
A2= 8.13298643510669E+00
A3=-4.10099875062160E+00
A4= 1.31571165602474E+00
A5=-2.78400484688017E-01
A6= 3.88294967051815E-02
A7=-3.43954747961506E-03
A8= 1.75404902262076E-04
A9=-3.91867465686797E-06

    ! Polvora KNSU 60/40
    case (2)
A0=-6.11043165506052E-01
A1= 2.98058454366955E+00
A2=-2.25564722929034E+00
A3= 9.67756473233875E-01
A4=-2.60214357140101E-01
A5= 4.54419230940925E-02
A6=-5.15377674951999E-03
A7= 3.66450558281284E-04
A8=-1.48496162006620E-05
A9= 2.61844906051524E-07
        
    ! Polvora KNSU 70/30
    case (3)
A0= 2.11131140257238E+00
A1=-1.48872601067561E+00
A2= 1.02571757620582E+00
A3=-4.04537302415782E-01
A4= 9.82597170151715E-02
A5=-1.51035569743172E-02
A6= 1.44690561421264E-03
A7=-8.09074405121067E-05
A8= 2.24135022850270E-06
A9=-1.84145422651684E-08

end select

logP0 = dlog10(P0)

! 0.01 bar a 11 kbar
k = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0


end subroutine Fk
!------------------------------------------------------------------------------------------
subroutine classe(Itot, Fmed, texto)

real*8, intent(in) :: Itot, Fmed

character*128, intent(out) :: texto

write(texto,'(F15.2)') Fmed

40 format(A, " - ", A)
if(Itot<=0.3125d0) then
    write(texto,40) "Micro", trim(adjustl(texto))
elseif(Itot<=0.625d0) then
    write(texto,40) "1/4 A", trim(adjustl(texto))
elseif(Itot<=1.250d0) then
    write(texto,40) "1/2 A", trim(adjustl(texto))
elseif(Itot<=2.500d0) then
    write(texto,40) " A", trim(adjustl(texto))
elseif(Itot<=5.000d0) then
    write(texto,40) " B", trim(adjustl(texto))
elseif(Itot<=10.00d0) then
    write(texto,40) " C", trim(adjustl(texto))
elseif(Itot<=20.00d0) then
    write(texto,40) " D", trim(adjustl(texto))
elseif(Itot<=40.00d0) then
    write(texto,40) " E", trim(adjustl(texto))
elseif(Itot<=80.00d0) then
    write(texto,40) " F", trim(adjustl(texto))
elseif(Itot<=160.0d0) then
    write(texto,40) " G", trim(adjustl(texto))
elseif(Itot<=320.0d0) then
    write(texto,40) " H", trim(adjustl(texto))
elseif(Itot<=640.0d0) then
    write(texto,40) " I", trim(adjustl(texto))
elseif(Itot<=1280.d0) then
    write(texto,40) " J", trim(adjustl(texto))
elseif(Itot<=2560.d0) then
    write(texto,40) " K", trim(adjustl(texto))
elseif(Itot<=5120.d0) then
    write(texto,40) " L", trim(adjustl(texto))    
elseif(Itot<=10200d0) then
    write(texto,40) " M", trim(adjustl(texto))  
endif
    
    
end subroutine classe
!------------------------------------------------------------------------------------------
subroutine M_cri(p0, Pamb, k, Mcri)
implicit none
real*8,intent(in) :: p0, pamb, k
real*8,intent(out) :: Mcri

Mcri = sqrt(((p0/Pamb)**((k-1.d0)/k)-1.d0)*2.d0/(k-1.d0))

end subroutine M_cri
!------------------------------------------------------------------------------------------
end Program empuxo1p4
