program Gibbs_1p3_main

! Cálculo de propriedades para reações químicas em equilíbrio e não-equilíbrio
! para o propelente H2/O2

! Implementação da versão atual (1.3):
!      Programador: Carlos Henrique Marchi
!      Local:       LENA/DEMEC/UFPR, Curitiba, PR
!      Período:     14 de janeiro a 24 de fevereiro de 2005
!      Última alteração: 3 de maio de 2005

! Implementação da versão atual (1.2):
!      Programador: Carlos Henrique Marchi
!      Local:       LENA/DEMEC/UFPR, Curitiba, PR
!      Período:     10 a 13 de janeiro de 2005
!      Última alteração: 13 de janeiro de 2005

! Implementação da versão atual (1.1):
!      Programador: Carlos Henrique Marchi
!      Local:       LENA/DEMEC/UFPR, Curitiba, PR
!      Período:     7 a 7 de janeiro de 2005
!      Última alteração: 7 de janeiro de 2005

! Implementação da versão inicial (1.0):
!      Programador: Carlos Henrique Marchi
!      Local:       LENA/DEMEC/UFPR, Curitiba, PR
!      Período:     4 a 6 de janeiro de 2005

!-------------------------------------------------------------------------------

! MODELO NUMÉRICO (resumo):
!   Até 8 espécies químicas disponíveis:
!     1 = H2O   2 = O2   3 = H2    4 = OH   
!     5 = O     6 = H    7 = HO2   8 = H2O2
!   Propriedades termodinâmicas e de trasnporte são calculadas conforme
!     NASA TM-4513 (1993)
!   O equilíbrio químico é resolvido com o método das constantes de	equilíbrio, 
!     que necessita definir as reações e as espécies envolvidas, e o método de
!     Newton simples para resolver as taxas de dissociação
!   A temperatura de combustão é resolvida com o método da bisseção
!   Linguagem Fortran
!   Aplicativo usado: Microsoft Developer Studio, Fortran PowerStation 4.0
!   Tipo de projeto: Console Application
!   O programa permite resolver 5 tipos de problemas:
!     1 = TERMO: calcular propriedades termodinâmicas e de transporte de cada
!                espécie química, na base mássica ou molar
!     2 = CONGELADO: calcular propriedades termodinâmicas e de transporte de
!                    uma mistura de espécies químicas
!     3 = EQUILIBRIO: calcular a composição de uma mistura de gases;
!                     há 9 modelos disponíveis, com 0 a 18 reações de
!                     dissociação, envolvendo 3 a 8 espécies químicas
!     4 = CÂMARA: calcular a temperatura de combustão e a composição de uma
!                 mistura gases; há 9 modelos disponíveis, com 0 a 18 reações de
!                 dissociação, envolvendo 3 a 8 espécies químicas
!     5 = TAXA FINITA: calcular a taxa de geração de massa por unidade de
!                      volume de cada espécie química;
!                      há 6 modelos disponíveis, com 4 a 18 reações de
!                      dissociação, envolvendo 6 ou 8 espécies químicas

! ARQUIVOS envolvidos no programa:
!   Gibbs_1p3_modulo.f90   = sub-rotinas funcionais
!   Gibbs_1p3_main.f90     = programa principal
!   Gibbs_coeficientes_cp_h_g.dat = coeficientes para propriedades termodinâmicas
!   Gibbs_coeficientes_ki_mi.dat  = coeficientes para propriedades de transporte
!   Gibbs_1p3.in              = arquivo de dados principal
!   Gibbs_1p3_termo.in        = arquivo de dados do bloco TERMO
!   Gibbs_1p3_congelado.in    = arquivo de dados do bloco CONGELADO
!   Gibbs_1p3_equilibrio.in   = arquivo de dados do bloco EQUILIBRIO
!   Gibbs_1p3_camara.in       = arquivo de dados do bloco CAMARA
!   Gibbs_1p3_taxa_finita.in  = arquivo de dados 1 do bloco TAXA FINITA
!   Gibbs_1p3_taxa_finita2.in = arquivo de dados 2 do bloco TAXA FINITA
!   *.*         = arquivo de saída do programa (nome definido pelo usuário)
!   Notepad.exe = editor de arquivos
!   Gibbs_1p3_LEIA_ME.txt = informações básicas sobre o programa
!   Gibbs_1p3.exe         = programa executável

!-------------------------------------------------------------------------------

use GIBBS_1p3

use PORTLIB

!-------------------------------------------------------------------------------

implicit none

!-------------------------------------------------------------------------------

! *** constantes ***

real*8,parameter :: po = 101325.0d+0 ! pressão ao nível do mar (Pa)

!-------------------------------------------------------------------------------

integer :: calculo ! Tipo de cálculo:
                   ! 1=TERMO; 2=CONGELADO; 3=EQUILÍBRIO; 4=CÂMARA; 5=TAXA FINITA


integer :: N       ! número de espécies químicas

integer :: base    ! base das propriedades: 1=mássica; 2=molar

integer :: i       ! número da espécie química
                   ! 1 = H2O   2 = O2   3 = H2    4 = OH   
                   ! 5 = O     6 = H    7 = HO2   8 = H2O2

real*8  :: T       ! temperatura (K)

real*8  :: p       ! pressão (atm)

real*8  :: p_si    ! pressão (Pa)

real*8  :: OF      ! razão em massa oxidante/combustível (adim.)

real*8  :: tol_e   ! tolerância para resolver e(j)

real*8  :: tol_n   ! tolerância para resolver n

real*8  :: tol_Tc  ! tolerância para obter Tc

integer :: itimax  ! número máximo de iterações para resolver e(j)

integer :: itemax  ! número máximo de iterações para resolver n

integer :: itcmax  ! número máximo de iterações para resolver Tc

integer :: modelo  ! número do modelo de reações em equilíbrio químico

real*8  :: T_O     ! temperatura de injeção do oxidante na câmara (K)

real*8  :: T_F     ! temperatura de injeção do combustível na câmara (K)

real*8  :: h_O     ! entalpia de injeção do oxidante na câmara (J/mol)

real*8  :: h_F     ! entalpia de injeção do combustível na câmara (J/mol)

integer :: tipo_hr ! tipo de HR: 1=automático; 2=prescrito

integer :: sai     ! =1, escreve variáveis secundárias durante as iterações

real*8,dimension(:),allocatable :: ni   ! número de moles (mol)

integer, parameter :: in   =  7 ! número do arquivo de entrada (dados)
integer, parameter :: out  = 10 ! número do arquivo de saída principal

character*50  :: caso      ! nome do arquivo de saída
character*90  :: titulo    ! título do gráfico
character*90  :: note_caso ! notepad + caso
integer       :: dos       ! acessa prompt dos

integer :: k      ! auxiliar

integer :: repete ! número de repetições dos cálculos para medir tcpu

real*8  :: tcpu   ! tempo de processamento (s)

!-------------------------------------------------------------------------------

! leitura dos dados

  dos = system('notepad Gibbs_1p3.in') ! lista dados

  open(in,file='Gibbs_1p3.in')

    read(in,*) calculo
    read(in,*) repete

  close(in)

  select case ( calculo )

    case ( 1 )
               dos = system('notepad Gibbs_1p3_Termo.in') ! lista dados
               open(in,file='Gibbs_1p3_Termo.in')
               read(in,*) caso
               read(in,*) N
               read(in,*) base
               read(in,*) T

    case ( 2 )
               dos = system('notepad Gibbs_1p3_Congelado.in') ! lista dados
               open(in,file='Gibbs_1p3_Congelado.in')
               read(in,*) caso
               read(in,*) N
               read(in,*) T
               allocate ( ni(N) )
               do i = 1, N
                  read(in,*) ni(i)
               end do

    case ( 3 )
               dos = system('notepad Gibbs_1p3_Equilibrio.in') ! lista dados
               open(in,file='Gibbs_1p3_Equilibrio.in')
               read(in,*) caso
               read(in,*) T
               read(in,*) p_si
               read(in,*) OF
               read(in,*) sai
               read(in,*) itimax
               read(in,*) tol_e
               read(in,*) itemax
               read(in,*) tol_n
               read(in,*) modelo

    case ( 4 )
               dos = system('notepad Gibbs_1p3_Camara.in') ! lista dados
               open(in,file='Gibbs_1p3_Camara.in')
               read(in,*) caso
               read(in,*) p_si
               read(in,*) OF
               read(in,*) sai
               read(in,*) itimax
               read(in,*) tol_e
               read(in,*) itemax
               read(in,*) tol_n
               read(in,*) itcmax
               read(in,*) tol_Tc
               read(in,*) tipo_hr
               read(in,*) T_O
               read(in,*) T_F
               read(in,*) h_O
               read(in,*) h_F
               read(in,*) modelo

    case ( 5 )
               dos = system('notepad Gibbs_1p3_Taxa_Finita.in') ! lista dados
               open(in,file='Gibbs_1p3_Taxa_Finita.in')
               read(in,*) caso
               read(in,*) T
               read(in,*) p_si
               read(in,*) sai
               read(in,*) modelo

    end select

    read(in,*) titulo
  
  close(in)

  p = p_si / po

!-------------------------------------------------------------------------------

! abertura do arquivo de saída

  open(out,file=caso)

  k = len(trim(adjustl(titulo)))

  write(out,1) trim(adjustl(titulo)), caso

  1 format(/, &
        7x, 'título = ', a<k>,//, &
       23x, 'caso = ', a50        )

!-------------------------------------------------------------------------------

! cálculos

  select case ( calculo )

    case ( 1 )
               call GIBBS_TERMO_executa

    case ( 2 )
               call GIBBS_CONGELADO_executa

    case ( 3 )
               call GIBBS_EQUILIBRIO_executa

    case ( 4 )
               call GIBBS_CAMARA_executa

    case ( 5 )
               call GIBBS_TAXA_FINITA_executa

    end select
  
!-------------------------------------------------------------------------------

  write(out,2) tcpu
  2 format(/, f14.3, ' = tempo de processamento (segundos)')

  close (out)

  note_caso = 'notepad '//caso

  dos = system(note_caso) ! lista arquivo de resultados


!-------------------------------------------------------------------------------


contains


!-------------------------------------------------------------------------------
 
  subroutine GIBBS_TERMO_executa

    call GIBBS_TERMO_dados_recebe ( out, N, base, T )
  
    call GIBBS_TERMO_dados_escreve ( out )

    call GIBBS_le_coeficientes_cp_h_g

    call GIBBS_escreve_coeficientes_cp_h_g ( out )

    ! cpi
    call GIBBS_cpi_aloca_memoria
    if ( base == 1 ) then
       call GIBBS_cpi_massa_calculo
    else
       call GIBBS_cpi_molar_calculo
    end if
    call GIBBS_cpi_escreve ( out )
    call GIBBS_cpi_fecha_memoria

    ! hi
    call GIBBS_hi_aloca_memoria
    if ( base == 1 ) then
       call GIBBS_hi_massa_calculo
    else
       call GIBBS_hi_molar_calculo
    end if
    call GIBBS_hi_escreve ( out )
    call GIBBS_hi_fecha_memoria

    ! gi
    call GIBBS_gi_aloca_memoria
    if ( base == 1 ) then
       call GIBBS_gi_massa_calculo
    else
       call GIBBS_gi_molar_calculo
    end if
    call GIBBS_gi_escreve ( out )
    call GIBBS_gi_fecha_memoria

    call GIBBS_le_coeficientes_ki_mii

    call GIBBS_escreve_coeficientes_ki_mii ( out )

    ! mii
    call GIBBS_mii_aloca_memoria
    call GIBBS_mii_calculo
    call GIBBS_mii_escreve ( out )
    call GIBBS_mii_fecha_memoria

    ! ki
    call GIBBS_ki_aloca_memoria
    call GIBBS_ki_calculo
    call GIBBS_ki_escreve ( out )
    call GIBBS_ki_fecha_memoria

    call GIBBS_fecha_memoria_ca_especie_Moli
    call GIBBS_fecha_memoria_cb_cc

  end subroutine GIBBS_TERMO_executa

!-------------------------------------------------------------------------------

  subroutine GIBBS_CONGELADO_executa

    call GIBBS_CONGELADO_dados_recebe ( out, N, T, ni )

    call GIBBS_CONGELADO_dados_escreve ( out )

    call GIBBS_CONGELADO_calculos1

    call GIBBS_CONGELADO_escreve1 ( out )

    call GIBBS_Yi_escreve ( out )

    call GIBBS_mi_kc_inicio

    call GIBBS_CONGELADO_calculos2

    call GIBBS_CONGELADO_escreve2 ( out )

    call GIBBS_Xi_escreve ( out )

  end subroutine GIBBS_CONGELADO_executa

!-------------------------------------------------------------------------------
 
  subroutine GIBBS_EQUILIBRIO_executa

    call GIBBS_EQUILIBRIO_dados_recebe ( out,   sai,   modelo, itimax, itemax, &
                                         tol_e, tol_n, T,      p,      OF      )

    call GIBBS_EQUILIBRIO_inicio1

    call GIBBS_EQUILIBRIO_dados_escreve ( out )

    tcpu = timef()

    do k = 1, repete

       call GIBBS_EQUILIBRIO_calculos1

    end do

    tcpu = timef()

    call GIBBS_EQUILIBRIO_escreve1

    call GIBBS_EQUILIBRIO_calculos2

    call GIBBS_EQUILIBRIO_escreve2

    call GIBBS_Ci_aloca_memoria

    call GIBBS_EQUILIBRIO_calculos3

    call GIBBS_EQUILIBRIO_escreve3

  end subroutine GIBBS_EQUILIBRIO_executa

!-------------------------------------------------------------------------------
 
  subroutine GIBBS_CAMARA_executa

    call GIBBS_CAMARA_dados_recebe ( out,    sai,    modelo, tipo_hr, itimax, &
                                     itemax, itcmax, tol_e,  tol_n,   tol_Tc, &
                                     T_O,    T_F,    h_O,    h_F,     p,      OF )

    call GIBBS_CAMARA_inicio

    call GIBBS_CAMARA_dados_escreve ( out )

    tcpu = timef()

    do k = 1, repete

       call GIBBS_CAMARA_calculos1

    end do

    tcpu = timef()

    call GIBBS_EQUILIBRIO_escreve1

    call GIBBS_EQUILIBRIO_calculos2

    call GIBBS_EQUILIBRIO_escreve2

    call GIBBS_Ci_aloca_memoria

    call GIBBS_EQUILIBRIO_calculos3

    call GIBBS_EQUILIBRIO_escreve3
   
    call GIBBS_CAMARA_escreve ( out )

  end subroutine GIBBS_CAMARA_executa

!-------------------------------------------------------------------------------
 
  subroutine GIBBS_TAXA_FINITA_executa

    call GIBBS_TAXA_FINITA_dados_recebe ( out, sai, modelo, T, p )

    call GIBBS_TAXA_FINITA_inicio

    call GIBBS_TAXA_FINITA_dados_escreve ( out )

    call GIBBS_escreve_coeficientes_cp_h_g ( out )

    call GIBBS_TAXA_FINITA_calculos1

    call GIBBS_TAXA_FINITA_escreve1

    tcpu = timef()

    do k = 1, repete

       call GIBBS_TAXA_FINITA_calculos2

    end do

    tcpu = timef()

    call GIBBS_TAXA_FINITA_escreve2

  end subroutine GIBBS_TAXA_FINITA_executa

!-------------------------------------------------------------------------------

end program
