module NUMERICO_EQUILIBRIO

! sigla = NUM_EQ_

! objetivo: calcular a solução numérica do modelo matemático
!           para propriedades e dados variáveis considerando
!           escoamento reativo em equilíbrio químico local

! versão original 1.0 (24 Jun 05)

! programador: Carlos Henrique Marchi
! local:       LENA/DEMEC - Curitiba - PR

! versão atual    1.0 (24 Jun 05)

! última alteração =    9 Fev 06

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

! RELAÇÃO DAS ROTINAS DESTE MÓDULO:

  ! NUM_EQ_cds_coef_pressao
  ! NUM_EQ_cds_coef_temperatura
  ! NUM_EQ_cds_coef_velocidade
  ! NUM_EQ_correcao_com_plinha
  ! NUM_EQ_cp_gama_congelados
  ! NUM_EQ_escreve_campos_auxiliares
  ! NUM_EQ_executa
  ! NUM_EQ_f
  ! NUM_EQ_fr_C_rec
  ! NUM_EQ_gama_efetivo
  ! NUM_EQ_grava_canal
  ! NUM_EQ_grava_reinicio
  ! NUM_EQ_h
  ! NUM_EQ_le_dados
  ! NUM_EQ_le_inicio
  ! NUM_EQ_le_reinicio
  ! NUM_EQ_le_rhg
  ! NUM_EQ_Mach
  ! NUM_EQ_massa_especifica_media
  ! NUM_EQ_mi_k_congelados
  ! NUM_EQ_pressao_entrada
  ! NUM_EQ_p_tot
  ! NUM_EQ_simplec
  ! NUM_EQ_T_tot
  ! NUM_EQ_T_wall
  ! NUM_EQ_uds_coef_pressao
  ! NUM_EQ_uds_coef_temperatura
  ! NUM_EQ_uds_coef_velocidade
  ! NUM_EQ_velocidade_media

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

! rotinas onde são lidos ARQUIVOS de dados neste MÓDULO:

  ! NUM_EQ_le_dados
  ! NUM_EQ_le_inicio
  ! NUM_EQ_le_reinicio
  ! NUM_EQ_le_rhg
  ! NUM_EQ_T_wall

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

! rotinas onde são gerados ARQUIVOS neste MÓDULO:

  ! NUM_EQ_executa
  ! NUM_EQ_grava_canal
  ! NUM_EQ_grava_reinicio

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

use ARQUIVO   ! versão 1.2 (19 Jan 00)

use SOLVERS   ! versão 1.0 (19 Jan 00)

use GIBBS_1p4 ! versão 1.4 ( 8 Nov 05)

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

implicit none

real*8, private :: p_cam       ! pressão na câmara de combustão (Pa)
real*8, private :: T_cam       ! temperatura na câmara de combustão (K)
real*8, private :: Rg_cte      ! constante do gás (J/kg.K)
real*8, private :: gama_cte    ! razão de calores específicos (adimensional)
real*8, private :: cp_cte      ! calor específico à pressão constante (J/kg.K)
real*8, private :: f_cte       ! fator de atrito de Darcy (adim)
real*8, private :: h_cte       ! coeficiente de convecção (W/m2.K)
real*8, private :: mi_cte      ! viscosidade absoluta (Pa.s)
real*8, private :: k_cte       ! condutividade térmica (W/m.K)
real*8, private :: T_wall_cte  ! temperatura da parede da tubeira (K)
real*8, private :: fr_cte      ! fator de recuperação (adim)
real*8, private :: epsilon_cte ! emissividade média parede-gases (adim)
real*8, private :: sigma       ! constante de Stefan-Boltzmann (W/m2.K4)
real*8, private :: rugoso      ! rugosidade absoluta da parede da tubeira (m)

integer,private :: tipo_cp     ! tipo de calor específico à pressão constante
integer,private :: tipo_f      ! tipo de fator de atrito de Darcy
integer,private :: tipo_h      ! tipo de coeficiente de convecção
integer,private :: tipo_mi     ! tipo de viscosidade absoluta
integer,private :: tipo_k      ! tipo de condutividade térmica
integer,private :: tipo_T_wall ! tipo de temperatura da parede da tubeira
integer,private :: tipo_fr     ! tipo de fator de recuperação

integer,private :: RHG    ! 0=executa MACH1D; 1=executa MACH1D acoplado ao CANAL
integer,private :: itg    ! número da iteração global corrente entre MACH1D e CANAL
integer,private :: itgmax ! número da iteração global máxima entre MACH1D e CANAL a realizar

integer,private :: i      ! número do volume de controle
integer,private :: itt    ! número da iteração no tempo

real*8, private :: u_in   ! velocidade na entrada da tubeira (m/s)
real*8, private :: p_in   ! pressão na entrada da tubeira que satisfaz QM (Pa)
real*8, private :: T_in   ! temperatura na entrada da tubeira (K)
real*8, private :: ro_in  ! massa específica na entrada da tubeira (m/s)

real*8, private :: p_ia   ! pressão antiga na entrada da tubeira (Pa)
real*8, private :: pl_in  ! variação da pressão na entrada da tubeira (Pa)

real*8, private :: p_ex   ! pressão na saída da tubeira (Pa)
real*8, private :: T_ex   ! temperatura na saída da tubeira (K)
real*8, private :: ro_ex  ! massa específica na saída da tubeira (m/s)

integer,private :: lista  ! lista conteúdo do arquivo de dados do módulo NUMERICO_EQUILIBRIO
integer,private :: le     ! lê conteúdo do arquivo de dados do módulo NUMERICO_EQUILIBRIO
integer,private :: copia  ! gera back up do arquivo de dados do módulo NUMERICO_EQUILIBRIO

character*50,private :: nome_num_var    ! nome do arquivo de dados do módulo NUMERICO_EQUILIBRIO

character*50,private :: T_wall_nome     ! nome do arquivo de dados de T_wall
character*50,private :: para_canal_nome ! nome do arquivo de saída para o programa CANAL

integer,private :: calor ! 0=sem troca de calor com a parede; 1=com troca de calor

! parte reativa

integer,private :: modelo_q ! número do modelo de reações em equilíbrio químico
integer,private :: itimax   ! número máximo de iterações para resolver e(j)
integer,private :: itemax   ! número máximo de iterações para resolver n
real*8, private :: OF       ! razão em massa oxidante/combustível (adim.)
real*8, private :: tol_e    ! tolerância para resolver e(j)
real*8, private :: tol_n    ! tolerância para resolver n
integer,private :: it_q     ! número da iteração global para iniciar os cálculos reativos
integer,private :: w_q      ! freqüência de iterações globais para fazer cálculos reativos
real*8, private :: dpY      ! variação da pressão para calcular gama

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

contains

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

  subroutine NUM_EQ_cds_coef_pressao ( n, beta, ro, pl, um, se, bp, bc )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: beta ! 0=UDS; 1=CDS

    real*8,dimension(:),intent(in)    :: ro  ! massa específica (kg/m3)
    real*8,dimension(:),intent(in)    :: pl  ! variação de pressão (Pa)
    real*8,dimension(:),intent(in)    :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in)    :: se ! área na face leste do volume de controle

    real*8,dimension(:),intent(inout) :: bp  ! termo independente

    real*8,dimension(:),intent(inout) :: bc  ! termo da correção adiada

    real*8 :: oeste, leste ! auxiliares

    ! volume 1 (fictício)

    bc(1) = 0.0d0

    ! volumes internos

    do i = 2, n-1

       oeste = se(i-1) * ( ro(i) - ro(i-1) ) * um(i-1)

       leste = se(i)   * ( ro(i+1) - ro(i) ) * um(i)

       bc(i) = 0.5d0 * beta * ( oeste - leste )

    end do

    ! volume n (fictício)

    bc(n) = 0.0d0

    ! atualiza o termo independente

    bp = bp + bc

  end subroutine NUM_EQ_cds_coef_pressao

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

  subroutine NUM_EQ_cds_coef_temperatura ( n, beta, cp, rom, T, um, se, bp, bc )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: beta ! 0=UDS; 1=CDS

    real*8,dimension(:),intent(in)    :: cp  ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),intent(in)    :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in)    :: T   ! temperatura (K) 
    real*8,dimension(:),intent(in)    :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in)    :: se  ! área na face leste do volume de controle

    real*8,dimension(:),intent(inout) :: bp  ! termo independente

    real*8,dimension(:),intent(inout) :: bc  ! termo da correção adiada

    real*8 :: oeste, leste ! auxiliares

    ! volume 1 (fictício)

    bc(1) = 0.0d0

    ! volumes internos

    do i = 2, n-1

       oeste = rom(i-1) * um(i-1) * se(i-1) * ( T(i) - T(i-1) )

       leste = rom(i) * um(i) * se(i) * ( T(i+1) - T(i) )

       bc(i) = 0.5d0 * beta * cp(i) * ( oeste - leste )

    end do

    ! volume n (fictício)

    bc(n) = 0.0d0

    ! atualiza o termo independente

    bp = bp + bc

  end subroutine NUM_EQ_cds_coef_temperatura

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

  subroutine NUM_EQ_cds_coef_velocidade ( n, beta, rom, u, um, se, bp, bc )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: beta ! 0=UDS; 1=CDS

    real*8,dimension(:),intent(in)    :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in)    :: u   ! velocidade (m/s)
    real*8,dimension(:),intent(in)    :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in)    :: se  ! área na face leste do volume de controle

    real*8,dimension(:),intent(inout) :: bp  ! termo independente

    real*8,dimension(:),intent(inout) :: bc  ! termo da correção adiada

    real*8 :: oeste, leste ! auxiliares

    ! volume 1 (fictício)

    bc(1) = 0.0d0

    ! volumes internos

    do i = 2, n-1

       oeste = rom(i-1) * um(i-1) * se(i-1) * ( u(i) - u(i-1) )

       leste = rom(i) * um(i) * se(i) * ( u(i+1) - u(i) )

       bc(i) = 0.5d0 * beta * ( oeste - leste )

    end do

    ! volume n (fictício)

    bc(n) = 0.0d0

    ! atualiza o termo independente

    bp = bp + bc

  end subroutine NUM_EQ_cds_coef_velocidade

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

  subroutine NUM_EQ_correcao_com_plinha ( n, xp, dp, de, pl, Rg, se, u, p, T, ro, um )

    integer,intent(in)                :: n  ! número de volumes de controle com dois fictícios

    real*8,dimension(:),intent(in)    :: xp ! coordenada espacial no centro do volume de controle
    real*8,dimension(:),intent(in)    :: se ! área na face leste do volume de controle
    real*8,dimension(:),intent(in)    :: dp ! coeficiente SIMPLEC no centro do volume
    real*8,dimension(:),intent(in)    :: de ! coeficiente SIMPLEC na face leste do volume
    real*8,dimension(:),intent(in)    :: pl ! variação de pressão (Pa)
    real*8,dimension(:),intent(in)    :: Rg ! constante da mistura de gases (J/kg.K)

    real*8,dimension(:),intent(inout) :: u  ! velocidade (m/s)
    real*8,dimension(:),intent(inout) :: p  ! pressão (Pa)
    real*8,dimension(:),intent(inout) :: T  ! temperatura (K)
    real*8,dimension(:),intent(inout) :: ro ! massa específica (kg/m3)
    real*8,dimension(:),intent(inout) :: um ! velocidade média (m/s)

    real*8 :: fator ! auxiliar

    ! pressão
    p = p + pl

    ! massa específica nodal
    ro = ro + pl / ( Rg * T )

    ! velocidade média
    do i = 1, n-1
       um(i) = um(i) - de(i) * ( pl(i+1) - pl(i) )
    end do

    ! velocidade na entrada da tubeira
    u_in = um(1)

    ! velocidade nodal
    do i = 2,n-1
       u(i) = u(i) - 0.5d0 * dp(i) * ( pl(i+1) - pl(i-1) )
    end do

    fator = 2 * xp(2) / ( xp(3) - xp(2) )
	u(1) = u(2) - fator * ( u(3) - u(2) )

    fator = 2 * ( xp(n) - xp(n-1) ) / ( xp(n-1) - xp(n-2) )
	u(n) = u(n-1) + fator * ( u(n-1) - u(n-2) )

  end subroutine NUM_EQ_correcao_com_plinha

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

  subroutine NUM_EQ_cp_gama_congelados ( n, T, p, cp, gama, Rg, Yi, hi, Nr, ej, itt, nitm )

    integer,intent(in) :: n  ! número de volumes de controle com dois fictícios

    real*8,dimension(:),  intent(in)    :: T     ! temperatura (K)
    real*8,dimension(:),  intent(in)    :: p     ! pressão (Pa)

    real*8,dimension(:),  intent(inout) :: cp    ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),  intent(inout) :: gama  ! razão de calores específicos (adimensional)

    real*8,dimension(:),  intent(inout) :: Rg    ! constante da mistura de gases (J/kg.K)
    real*8,dimension(:,:),intent(inout) :: Yi    ! fração mássica (adim.)
    real*8,dimension(:,:),intent(inout) :: hi    ! entalpia mássica (J/kg)

    integer,intent(in) :: Nr   ! número de reações químicas do modelo adotado
    integer,intent(in) :: itt  ! iteração do ciclo principal
    integer,intent(in) :: nitm ! número máximo de iterações do ciclo principal
    real*8,dimension(:,:),intent(inout) :: ej    ! taxas de reação
    real*8,dimension(Nr) ::  e_j  ! variável auxiliar para reações químicas
    
    integer :: j ! contador (número de reações)

    integer :: aux_1, aux_2 

    !call GIBBS_MACH1D_abre_ej
    
    ! passagem da iteração executada para a subrotina Gibbs
    call GIBBS_MACH1D_recebe4 (itt)
    call GIBBS_MACH1D_recebe5 (n)
    call GIBBS_MACH1D_recebe7 (nitm)
    call GIBBS_MACH1D_alocacao_ej

    do i = 1, n

       !Passagem de dados para o código Gibbs
       call GIBBS_MACH1D_recebe3 (i)

       call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modelo_q, itimax, itemax,  &
                                            tol_e, tol_n, T(i), p(i)/101325.0d0, OF )
                                            
       if (itt > 1) call GIBBS_MACH1D_recebe6 (ej)
                                            
       call GIBBS_EQUILIBRIO_calculos1
       
       call GIBBS_MACH1D_envia7 (e_j)
       
       do j = 1, Nr
          ej (i,j) = e_j(j)
       end do
       
       call GIBBS_hi_massa_calculo
       call GIBBS_MACH1D_envia4 ( cp(i), gama(i), Rg(i), Yi(i,:), hi(i,:) )

    end do
    
    !call GIBBS_MACH1D_fecha_ej
    call GIBBS_MACH1D_desalocacao_ej
    
  end subroutine NUM_EQ_cp_gama_congelados

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

  subroutine NUM_EQ_escreve_campos_auxiliares ( unit, n,   w,      epsilon, xp, qlr,  qlt,  &
                                                f,    h,   T_wall, Taw,     cp, gama, Red,  &
                                                Pr,   mi,  k,      C_rec,   fr, xe,   pl,   &
                                                um,   rom, Ne,     especie, Yi, Rg,   ent, gamaf )

    integer,intent(in) :: unit ! unidade de leitura do arquivo de dados
    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: w    ! freqüência de escrita de campos auxiliares
    integer,intent(in) :: Ne   ! número de espécies químicas

    real*8,intent(in) :: epsilon ! emissividade média parede-gases (adim)

    real*8,dimension(:),intent(in) :: xp     ! coordenada espacial no centro do volume de controle
    real*8,dimension(:),intent(in) :: xe     ! coordenada espacial na face leste do volume de controle
    real*8,dimension(:),intent(in) :: qlr    ! fluxo de calor radiativo para a parede (W/m2)
    real*8,dimension(:),intent(in) :: qlt    ! fluxo de calor total para a parede (W/m2)
    real*8,dimension(:),intent(in) :: f      ! fator de atrito de Darcy (adim)
    real*8,dimension(:),intent(in) :: h      ! coeficiente de convecção (W/m2.K)
    real*8,dimension(:),intent(in) :: T_wall ! temperatura da parede da tubeira (K)
    real*8,dimension(:),intent(in) :: Taw    ! temperatura de parede adiabática (K)
    real*8,dimension(:),intent(in) :: cp     ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),intent(in) :: ent    ! entropia (J/kg.K)
	real*8,dimension(:),intent(in) :: gama   ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in) :: gamaf  ! gama congelado (adimensional)
    real*8,dimension(:),intent(in) :: Rg     ! constante da mistura de gases (J/kg.K)
    real*8,dimension(:),intent(in) :: Red    ! número de Reynolds (adim)
    real*8,dimension(:),intent(in) :: Pr     ! número de Prandtl (adim)
    real*8,dimension(:),intent(in) :: mi     ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),intent(in) :: k      ! condutividade térmica (W/m.K)
    real*8,dimension(:),intent(in) :: C_rec  ! coeficiente de recuperação de Taw
    real*8,dimension(:),intent(in) :: fr     ! fator de recuperação (adim)
    real*8,dimension(:),intent(in) :: pl     ! fluxo de calor total para a parede (W/m2)
    real*8,dimension(:),intent(in) :: um     ! velocidade média (m/s)
    real*8,dimension(:),intent(in) :: rom    ! massa específica média (Kg/m3)

    character*4,dimension(:),intent(in) :: especie ! nome da espécie química

    real*8,dimension(:,:),intent(in)    :: Yi      ! fração mássica (adim.)

    integer      :: ii    ! número da espécie
    character*21 :: aux

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** fluxo de calor na parede ***                   ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'q"_wh (W/m2)         ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            qlt                      ) ! variável a escrever

    if ( h(1) > 0.0d0 .or. epsilon > 0.0d0 ) then
       call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                               w,                       & ! freqüência de escrita
                               '*** fluxo radiativo / (radiativo+advectivo) ***    ', & ! 50 colunas
                               ' ponto',                & !  6 colunas
                               'x (m)                ', & ! 21 colunas
                               'razão (adim.)        ', & ! 21 colunas
                               2,                       & ! número do ponto inicial
                               n-1,                     & ! número de dados a escrever
                               xp,                      & ! coordenada espacial
                               qlr/qlt                  ) ! variável a escrever
    end if

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** fator de atrito de Darcy ***                   ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'f (adimensional)     ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            f                        ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** coeficiente de convecção ***                   ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'h (W/m2.K)           ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            h                        ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** temperatura da parede da tubeira ***           ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'T_wall (K)           ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            T_wall                   ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** temperatura de parede adiabática ***           ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'Taw (K)              ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            Taw                      ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** cp = calor específico à pressão constante ***  ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'cp (J/kg.K)          ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n,                       & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            cp                       ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** razão de calores específicos ***               ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'gama (adimensional)  ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n,                       & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            gama                     ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** gama congelado ***                             ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'gamaf (adimensional) ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n,                       & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            gamaf                    ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** constante da mistura de gases ***              ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'Rg (J/kg.K)          ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n,                       & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            Rg                       ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** número de Reynolds baseado no diâmetro ***     ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'Re (adimensional)    ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            Red                      ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** número de Prandtl ***                          ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'Pr (adimensional)    ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            Pr                       ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** viscosidade absoluta ***                       ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'mi (Pa.s)            ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            mi                       ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** condutividade térmica ***                      ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'k (W/m.K)            ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            k                        ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** coeficiente de recuperação de Taw ***          ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'C_rec (adimensional) ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            C_rec                    ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** fator de recuperação ***                       ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'fr (adimensional)    ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            fr                       ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** plinha ***                                     ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'pl (Pa)              ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n,                       & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            pl                       ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** velocidade nas faces ***                       ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'um                   ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xe,                      & ! coordenada espacial
                            um                       ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** massa específica nas faces ***                 ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'rom (kg/m3)          ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xe,                      & ! coordenada espacial
                            rom                      ) ! variável a escrever

    do ii = 1, Ne
       aux = ' *** espécie ' // especie(ii) // ' ***'
       call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                               w,                       & ! freqüência de escrita
                               aux                   ,  & ! 50 colunas
                               ' ponto',                & !  6 colunas
                               'x (m)                ', & ! 21 colunas
                               'Yi (adim.)           ', & ! 21 colunas
                               1,                       & ! número do ponto inicial
                               n,                       & ! número de dados a escrever
                               xp,                      & ! coordenada espacial
                               Yi(:,ii)                 ) ! variável a escrever
    end do

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** entropia ***                                   ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'entropia (J/kg.K)    ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n,                       & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            ent                      ) ! variável a escrever

  end subroutine NUM_EQ_escreve_campos_auxiliares

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

  subroutine NUM_EQ_executa ( n,         w_it,  beta,     nitm,    tipo_ci, nome_ci, w_aux, &
                              lista_aux, dt,    fm_exato, F_exato, rc,      Sg,      xe,    &
                              xp,        rp,    se,       sp,      sn,      fm,      F_ex,  &
                              r_qwh,     q_sai, u,        p,       T,       ro,      M,     &
                              p_tot,     T_tot, itgr,     itgmaxr  )

    integer,intent(in) :: n         ! número de volumes de controle com dois fictícios
    integer,intent(in) :: w_it      ! freqüência de escrita de parâmetros ao longo das iterações
    integer,intent(in) :: beta      ! 0=UDS; 1=CDS
    integer,intent(in) :: nitm      ! número máximo de iterações no tempo
    integer,intent(in) :: tipo_ci   ! tipo de condição inicial: 1=analítica; 2=reinício
    integer,intent(in) :: w_aux     ! freqüência de escrita de campos auxiliares
    integer,intent(in) :: lista_aux ! lista arquivo com os campos auxiliares (1); 0=não
 
    character*(*),intent(in) :: nome_ci ! nome do arquivo de dados de reinício

    real*8,intent(in) :: dt       ! intervalo de tempo (s)
    real*8,intent(in) :: fm_exato ! fluxo de massa exato (kg/s)
    real*8,intent(in) :: F_exato  ! empuxo na saída da tubeira exato (N)   
    real*8,intent(in) :: rc       ! raio de curvatura na garganta da tubeira (m)
    real*8,intent(in) :: Sg       ! área da garganta (m2)

    real*8,dimension(:),intent(in) :: xe ! coordenada espacial na face leste do volume de controle
    real*8,dimension(:),intent(in) :: xp ! coordenada espacial no centro do volume de controle
    real*8,dimension(:),intent(in) :: rp ! raio no centro do volume de controle
    real*8,dimension(:),intent(in) :: se ! área na face leste do volume de controle
    real*8,dimension(:),intent(in) :: sp ! área no centro do volume de controle
    real*8,dimension(:),intent(in) :: sn ! área da parede da tubeira

    real*8,intent(out) :: fm    ! fluxo de massa numérico (kg/s)
    real*8,intent(out) :: F_ex  ! empuxo na saída da tubeira numérico (N)
    real*8,intent(out) :: r_qwh ! razão q_wall / q_sai
    real*8,intent(out) :: q_sai ! taxa de transferência de calor advectivo na saída da tubeira (W)

    real*8,dimension(:),intent(inout) :: u     ! velocidade (m/s)
    real*8,dimension(:),intent(inout) :: p     ! pressão (Pa)
    real*8,dimension(:),intent(inout) :: T     ! temperatura (K)
    real*8,dimension(:),intent(inout) :: ro    ! massa específica (kg/m3)
    real*8,dimension(:),intent(inout) :: M     ! número de Mach (adimensional)
    real*8,dimension(:),intent(inout) :: p_tot ! pressão de estagnação (Pa)
    real*8,dimension(:),intent(inout) :: T_tot ! temperatura de estagnação (K)

    integer,intent(inout) :: itgr    ! número da iteração global corrente entre MACH1D e CANAL
    integer,intent(inout) :: itgmaxr ! número da iteração global máxima entre MACH1D e CANAL a realizar

    integer :: unit_it = 81 ! identificador do arquivo das iterações
    integer :: unit_au = 82 ! identificador do arquivo dos campos auxiliares

    real*8 :: Cd     ! coeficiente de descarga (adimensional)
    real*8 :: F_star ! empuxo dinâmico (adimensional)

    real*8 :: tcpu1 ! tempo de processamento
    real*8 :: tcpu2 ! tempo de processamento

    character*(50) :: caso    ! nome do caso

    character*220  :: comando ! linha de comando

! *** campos auxiliares ***

    real*8,dimension(:),allocatable :: aw ! coeficiente oeste
    real*8,dimension(:),allocatable :: ap ! coeficiente central
    real*8,dimension(:),allocatable :: ae ! coeficiente leste
    real*8,dimension(:),allocatable :: bp ! termo independente
    real*8,dimension(:),allocatable :: bc ! termo da correção adiada
    real*8,dimension(:),allocatable :: bf ! termo de fricção para um

    real*8,dimension(:),allocatable :: ua  ! velocidade anterior (m/s)
    real*8,dimension(:),allocatable :: pa  ! pressão anterior (Pa)
    real*8,dimension(:),allocatable :: ta  ! temperatura anterior (K)
    real*8,dimension(:),allocatable :: roa ! massa específica anterior (kg/m3)
    real*8,dimension(:),allocatable :: uma ! velocidade média anterior (m/s)
    real*8,dimension(:),allocatable :: um  ! velocidade média (m/s)
    real*8,dimension(:),allocatable :: rom ! massa específica média (Kg/m3)
    real*8,dimension(:),allocatable :: pl  ! variação de pressão (Pa)

    real*8,dimension(:),allocatable :: dp ! coeficiente SIMPLEC no centro do volume
    real*8,dimension(:),allocatable :: de ! coeficiente SIMPLEC na face leste do volume

    real*8,dimension(:),allocatable :: C_rec ! coeficiente de recuperação de Taw
    real*8,dimension(:),allocatable :: Taw   ! temperatura de parede adiabática (K)
    real*8,dimension(:),allocatable :: qlt   ! fluxo de calor total para a parede (W/m2)
    real*8,dimension(:),allocatable :: qlr   ! fluxo de calor radiativo para a parede (W/m2)

	real*8,dimension(:),allocatable :: gama   ! razão de calores específicos (adimensional)
	real*8,dimension(:),allocatable :: gamaf  ! gama congelado (adimensional)
    real*8,dimension(:),allocatable :: cp     ! calor específico médio à pressão constante (J/kg.K)
    real*8,dimension(:),allocatable :: Rg     ! constante da mistura de gases (J/kg.K)
    real*8,dimension(:),allocatable :: f      ! fator de atrito de Darcy (adim)
    real*8,dimension(:),allocatable :: h      ! coeficiente de convecção (W/m2.K)
    real*8,dimension(:),allocatable :: mi     ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),allocatable :: k      ! condutividade térmica (W/m.K)
    real*8,dimension(:),allocatable :: T_wall ! temperatura da parede da tubeira (K)
    real*8,dimension(:),allocatable :: fr     ! fator de recuperação (adim)
    real*8,dimension(:),allocatable :: Red    ! número de Reynolds (adim)
    real*8,dimension(:),allocatable :: Pr     ! número de Prandtl (adim)

    integer :: Ne ! número de reações químicas
    integer :: ii ! número da espécie
    character*4,dimension(:),allocatable :: especie ! nome da espécie química
    real*8,dimension(:,:),allocatable :: Yi  ! fração mássica (adim.)
    real*8,dimension(:),allocatable   :: ent ! entropia (J/kg.K)
    real*8,dimension(:,:),allocatable :: hi  ! entalpia mássica (J/kg)

    integer :: Nr ! número de reações químicas
    
    real*8,dimension(:,:),allocatable :: ej   ! taxa de dissociação (GIBBS)


    real*8  :: ds  ! variação da entropia entre a câmara e a entrada do domínio
    integer :: ite
    
    ! alocação de memória  
    allocate ( ua(n), pa(n), Ta(n), pl(n), roa(n), dp(n) )
    allocate ( um(n-1), uma(n-1), de(n-1), rom(n-1) )
    allocate ( aw(n), ap(n), ae(n), bp(n) )
    allocate ( bc(n), bf(n) )
    allocate ( C_rec(n) )
    allocate ( gama(n), cp(n),  f(n),  h(n), mi(n), k(n), T_wall(n), &
               fr(n),   Red(n), Pr(n),  Rg(n), ent(0:n), gamaf(n) )

    fm = fm_exato

    ! inicialização de campos
    p     = p_cam/((1.0d0+(gama_cte-1.0d0)*(M**2)/2.0d0)**(gama_cte/(gama_cte-1.0d0)))
    ro    = p / ( Rg_cte * T )
    um    = u(1:n-1)
    pl    = 0.0d0
    rom   = ro(1:n-1)
    pl_in = 0.0d0
    bc    = 0.0d0

    ! inicialização na entrada da tubeira
    u_in  = um(1)
    T_in  = T_cam - 0.5d0*(gama_cte-1.0d0)*(u_in**2)/(gama_cte*Rg_cte)
    call NUM_EQ_pressao_entrada ( gama_cte, Rg_cte, p_cam, u_in, T_in, p_in )
    p(1)  = 2.0d0*p_in - p(2)
    u(1)  = 2*u(2) - u(3)
    T(1)  = -T(2) + 2*T_in
    ro(1) = p(1) / ( Rg_cte * T(1) )

    ! inicialização na saída da tubeira
    p(n)  = 2.0d0*p(n-1) - p(n-2)
    u(n)  = 2.0d0*u(n-1) - u(n-2)
    T(n)  = 2.0d0*T(n-1) - T(n-2)
    ro(n) = p(n) / ( Rg_cte * T(n) )

    itg    = itgr
    itgmax = itgmaxr
    if ( RHG == 1 ) call NUM_EQ_le_rhg ( unit_au )

    ! reinicialização dos campos
    if ( tipo_ci==2 .or. itg > 1 ) &
       call NUM_EQ_le_reinicio ( n, unit_au, nome_ci, u,  p,  T,  ro,  um, rom, pl )

    ! cria o arquivo das iterações
    if ( w_it /= 0 .and. itg == itgmax ) then
       call ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.mach1d.itn'
       call ARQUIVO_cria_arquivo ( unit_it, comando )
       write ( unit_it, 1 ) 
       1 format (t8,'itt', t13,'cd',   t36, 'F_star', t59, 'T_ex', &
                           t82,'u_in', t105,'u_ex',   t128,'Y(H2O,n-1)' )
    end if

    ! lê arquivo de dados do programa CANAL com T_wall
    if ( itg > 1 ) tipo_T_wall = 2
    call NUM_EQ_T_wall ( unit_au, n, T_wall )

    ! inicialização do módulo químico e cálculo de propriedades (Yi etc) na câmara
    if ( tipo_cp == 2 ) then
       call GIBBS_EQUILIBRIO_dados_recebe ( unit_au, 0, modelo_q, itimax, itemax,  &
                                            tol_e, tol_n, T_cam, p_cam/101325.0d0, OF )
       
       call GIBBS_EQUILIBRIO_inicio1
       call GIBBS_EQUILIBRIO_calculos1
       call GIBBS_MACH1D_envia2 ( Ne )
       call GIBBS_MACH1D_envia8 ( Nr )
       allocate ( Yi(n,Ne), especie(Ne), hi(n,Ne), ej(n,Nr) )
       call GIBBS_MACH1D_inicio1
       call GIBBS_hi_massa_calculo
       call GIBBS_MACH1D_envia3 ( especie )
       call GIBBS_MACH1D_envia4 ( cp_cte, gama_cte, Rg_cte, Yi(1,:), hi(1,:) )
       call NUM_EQ_cp_gama_congelados ( n, T, p, cp, gama, Rg, Yi, hi, Nr, ej, itt, nitm )
       call GIBBS_mi_kc_inicio
       gamaf = gama
       call NUM_EQ_gama_efetivo ( n, T, p, gamaf, Rg, gama )
    else
       allocate ( Yi(n,1), especie(1), hi(n,1) )
       Yi = 0.0d0
       Ne = 1
       especie = 'nada'
       cp   = cp_cte
       gama = gama_cte
       gamaf = gama
       Rg = Rg_cte
    end if

    ! inicialização das propriedades
    call NUM_EQ_mi_k_congelados ( n, T, mi, k )
    Red = ro * u * 2 * rp / mi
    call NUM_EQ_f ( rp, Red, f )
    Pr  = cp * mi / k
    call NUM_EQ_Mach ( Rg, gama, u, T, M )
    call NUM_EQ_fr_C_rec ( gama, Pr, M, fr, C_rec )
    call NUM_EQ_h ( n, fm, rc, Sg, sp, M, gama, cp, mi, T_wall, Pr, h )

    call ARQUIVO_cronometro ( tcpu1 )

    ! início do ciclo iterativo

    do itt = 1, nitm 

       ! atualização das incógnitas do instante anterior
       ua  = u
       pa  = p
       Ta  = T
       roa = ro
       uma = um

       ! atualização da pressão na entrada da tubeira
       p_ia = p_in
       call NUM_EQ_pressao_entrada ( gama_cte, Rg_cte, p_cam, u_in, T_in, p_in )
       p(1) = 2.0d0*p_in - p(2)
       pl_in = p_in - p_ia
     
       ! cálculo das propriedades

       if ( itt >= it_q .and. ( itt == nitm .or. mod(itt,w_q)==0 ) ) then

          if ( tipo_cp == 2 ) call NUM_EQ_cp_gama_congelados ( n, T, p, cp, gama, Rg, Yi, hi, Nr, ej, itt, nitm )

          call NUM_EQ_mi_k_congelados ( n, T, mi, k )

          Red = ro * u * 2 * rp / mi

          call NUM_EQ_f ( rp, Red, f )

          if ( calor == 1 ) then
             Pr  = cp * mi / k
             gamaf = gama
             call NUM_EQ_gama_efetivo ( n, T, p, gamaf, Rg, gama )
             call NUM_EQ_Mach ( Rg, gama, u, T, M )
             call NUM_EQ_fr_C_rec ( gama, Pr, M, fr, C_rec )
             call NUM_EQ_h ( n, fm, rc, Sg, sp, M, gama, cp, mi, T_wall, Pr, h )
          end if

       end if

       ! solução da velocidade nodal
       call NUM_EQ_uds_coef_velocidade ( n,  dt, f,  ro, roa, u,  ua, um, p,  se, &
                                         sp, rp, xe, xp, rom, aw, ap, ae, bp, bf  )
       if ( beta > 0 ) &
          call NUM_EQ_cds_coef_velocidade ( n, beta, rom, u, um, se, bp, bc )
       call SOLVERS_tdma ( n, aw, ap, ae, bp, u )

       ! cálculo dos coeficientes do SIMPLEC para pl
       call NUM_EQ_simplec ( n, sp, aw, ap, ae, dp, de )

       ! cálculo da velocidade média
       call NUM_EQ_velocidade_media ( n,  dt, roa, u,  uma, p,  se, &
                                      sp, xe, aw,  ap, ae,  bc, bf, um )

       ! solução da temperatura
       T_in = T_cam - 0.5d0*(gama_cte-1.0d0)*(u_in**2)/(gama_cte*Rg_cte)

       call NUM_EQ_uds_coef_temperatura ( n,  Ne, dt, T_in, epsilon_cte, sigma, &
                                          ro, roa,   rom,  u,      um,  T,      &
                                          Ta, p,     pa,   se,     sp,  rp,     &
                                          xe, xp,    sn,   f,      h,   T_wall, &
                                          cp, C_rec, aw, ap, ae, bp, Yi, hi )

       if ( beta > 0 ) &
          call NUM_EQ_cds_coef_temperatura ( n, beta, cp, rom, T, um, se, bp, bc )
       call SOLVERS_tdma ( n, aw, ap, ae, bp, T )

       ! cálculo da massa específica
       ro = p / ( Rg * T )
       call NUM_EQ_massa_especifica_media ( n, beta, ro, rom )

       ! solução da correção da pressão
       call NUM_EQ_uds_coef_pressao ( n,  dt, Rg, pl_in, ro, rom, roa, um, T, de, &
                                      se, sp, xe, xp,    pl, aw,  ap,  ae, bp     )
       if ( beta > 0 ) &
          call NUM_EQ_cds_coef_pressao ( n, beta, ro, pl, um, se, bp, bc )
       call SOLVERS_tdma ( n, aw, ap, ae, bp, pl )

       ! correções das variáves com plinha
       call NUM_EQ_correcao_com_plinha ( n, xp, dp, de, pl, Rg, se, u, p, T, ro, um )
       call NUM_EQ_massa_especifica_media ( n, beta, ro, rom )

       ! cálculo das variáveis monitoradas ao longo das iterações
       fm     = rom(1) * u_in * se(1)
       F_ex   = fm * um(n-1)
       cd     = fm / fm_exato
       F_star = F_ex / F_exato
       T_ex   = 0.5d0 * ( T(n-1) + T(n) )

       ! escreve variáveis monitoradas
       if ( w_it /= 0 .and. itg == itgmax ) then
          if ( itt==1 .or. itt==nitm .or. mod(itt,w_it)==0 ) &
             write ( unit_it, 2 ) itt, cd, F_star, T_ex, u_in, um(n-1), Yi(n-1,1)
             2 format ( i10, 6 ( 1pe23.15 ) )
       end if
       
       ! renomeação dos arquivos auxiliares
       !call GIBBS_MACH1D_renomeia_ej

    end do

    call ARQUIVO_cronometro ( tcpu2 )

    if ( w_it /= 0 .and. itg == itgmax ) then
       call ARQUIVO_escreve_tempo_cpu ( tcpu1, tcpu2, unit_it )
       close ( unit_it )
    end if

    call NUM_EQ_grava_reinicio ( n, unit_au, nome_ci, u,  p,  T,  ro,  um, rom, pl )

    ! entrada da tubeira
    u(1)  = u_in
    p(1)  = p_in
    T(1)  = T_in

    ! saída da tubeira
    T_ex  = 0.5d0 * ( T(n-1) + T(n) )
    p_ex  = 0.5d0 * ( p(n-1) + p(n) )
    u(n)  = um(n-1)
    p(n)  = p_ex
    T(n)  = T_ex

    deallocate ( ua, pa, ta, roa, uma, dp, de )
    deallocate ( aw, ap, ae, bp,  bf )

    ! cria o arquivo dos campos auxiliares
    if ( w_aux /= 0 .and. itg == itgmax ) then
       call ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.mach1d.outros.txt'
       call ARQUIVO_cria_arquivo ( unit_au, comando )
       write ( unit_au, 3 ) 
       3 format (/,t5,'*** CAMPOS AUXILIARES ***')
       ! escreve dados do módulo químico
       if ( tipo_cp > 1 ) then
          write ( unit_au, 4 ) 
          4 format (/,t5,'*** Condições químicas na câmara ***')
          call GIBBS_EQUILIBRIO_dados_recebe ( unit_au, 0, modelo_q, itimax, itemax,  &
                                               tol_e, tol_n, T_cam, p_cam/101325.0d0, OF )
          call GIBBS_MACH1D_alocacao_ej
          call GIBBS_EQUILIBRIO_calculos1
          call GIBBS_MACH1D_desalocacao_ej
          call GIBBS_si_aloca_memoria
          call GIBBS_MACH1D_camara_entropia ( ent(0) )
          call GIBBS_EQUILIBRIO_dados_escreve ( unit_au )
          call GIBBS_EQUILIBRIO_escreve1
          call GIBBS_CONGELADO_calculos2
          call GIBBS_CONGELADO_escreve2 ( unit_au )
          call NUM_EQ_cp_gama_congelados ( n, T, p, cp, gama, Rg, Yi, hi, Nr, ej, itt, nitm )
          gamaf = gama
          call NUM_EQ_gama_efetivo ( n, T, p, gamaf, Rg, gama )
          do i = 1, n
             call GIBBS_MACH1D_entropia ( T(i), p(i)/101325.0d0, Rg(i), Yi(i,:), ent(i) )
          end do
       end if
    end if

    ro_in = p_in / ( Rg(1) * T_in )
    ro(1) = ro_in
    ro_ex = p_ex / ( Rg(n) * T_ex )
    ro(n) = ro_ex

    ! calcula os fluxos de calor para a parede
    allocate ( Taw(n), qlt(n), qlr(n) )
    Taw = C_rec * T
    qlr = epsilon_cte * sigma * ( (T_wall**4) - (T**4) )
    qlt = qlr + h*(T_wall-Taw)

    if ( w_aux /= 0 .and. itg == itgmax ) then

       ds = cp_cte * dlog ( T_in / T_cam ) - Rg_cte * dlog ( p_in / p_cam )

       write ( unit_au, 5 ) ds, ent(0)
       5 format ( //, '*** variação da entropia entre a estagnação e a entrada = ', 1pe23.15, &
                  //, '*** entropia na estagnação = ', 1pe23.15, / )

       call NUM_EQ_escreve_campos_auxiliares ( unit_au, n, w_aux,  epsilon_cte, xp, qlr, &
                                               qlt, f,  h,   T_wall, Taw, cp, gama, Red, &
                                               Pr, mi,  k,  C_rec, fr, xe, pl, um, rom,  &
                                               Ne, especie, Yi, Rg, ent, gamaf )

    end if

    close ( unit_au )

    ! Yi e hi são baseados em T_ex
    q_sai = 0.0d0
    do ii = 1, Ne
       q_sai = q_sai + hi(n,ii) * Yi(n,ii)
    end do
    q_sai = - fm * q_sai

    ! cálculo da taxa total de transferência de calor para a parede
    r_qwh = 0.0d0
    do i = 2, n-1
       r_qwh = r_qwh + qlt(i)*sn(i)
    end do

    r_qwh = r_qwh / q_sai

    call NUM_EQ_Mach  ( Rg, gama, u, T, M )
    call NUM_EQ_p_tot ( n,  gama, M, p, p_tot )
    call NUM_EQ_T_tot ( n,  gama, M, T, T_tot )

    ! gera arquivo de dados para o programa CANAL
    if ( RHG == 1 ) &
       call NUM_EQ_grava_canal ( n, unit_au, para_canal_nome, Cd, F_star, qlt, &
                                    T_wall,  Taw,             sn, h,      T    )
  
    deallocate ( C_rec, Taw, qlt, qlr, pl, um, rom, Yi, especie, Rg, ent, hi, gamaf )
    deallocate ( gama,  cp,  f,   h,  mi, k, T_wall, fr, Red, Pr )

    itgr    = itg
    itgmaxr = itgmax

  end subroutine NUM_EQ_executa

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

  subroutine NUM_EQ_f ( rp, Red, f )

    real*8,dimension(:),intent(in)  :: rp  ! raio no centro do volume de controle
    real*8,dimension(:),intent(in)  :: Red ! número de Reynolds (adim)

    real*8,dimension(:),intent(out) :: f   ! fator de atrito de Darcy (adim)

    select case ( tipo_f )

       case ( 1 )
          f = f_cte

       case ( 2 )
          f = 1.0d0 / ( 4.0d0 * ( ( dlog ( rugoso/(2.0d0*3.7d0*rp) + 5.74d0/(Red**0.9d0) ) ) ** 2 ) )

    end select
    
  end subroutine NUM_EQ_f

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

  subroutine NUM_EQ_fr_C_rec ( gama, Pr, M, fr, C_rec )

    real*8,dimension(:),intent(in)  :: gama  ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in)  :: Pr    ! número de Prandtl (adim)
    real*8,dimension(:),intent(in)  :: M     ! número de Mach (adimensional)

    real*8,dimension(:),intent(out) :: fr    ! fator de recuperação (adim)
    real*8,dimension(:),intent(out) :: C_rec ! coeficiente de recuperação de Taw

    select case ( tipo_fr )

       case ( 1 )
          fr = fr_cte

       case ( 2 )
          fr = Pr**(1/3.0d0)

    end select

    C_rec = 1.0d0 + fr*(gama-1.0d0)*(M**2)/2.0d0
    
  end subroutine NUM_EQ_fr_C_rec

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

  subroutine NUM_EQ_gama_efetivo ( n, T, p, gamaf, Rg, gama )

    integer,intent(in) :: n ! número de volumes de controle com dois fictícios

    real*8,dimension(:),intent(in)    :: T     ! temperatura (K)
    real*8,dimension(:),intent(in)    :: p     ! pressão (Pa)
    real*8,dimension(:),intent(in)    :: gamaf ! gama congelado (adimensional)
    real*8,dimension(:),intent(in)    :: Rg    ! constante da mistura de gases (J/kg.K)

    real*8,dimension(:),intent(inout) :: gama  ! razão de calores específicos (adimensional)

    do i = 1, n
       call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modelo_q, itimax, itemax, tol_e, &
                                            tol_n, T(i), (p(i)+dpY*p(i)/2.0d0)/101325.0d0, OF )
       call GIBBS_EQUILIBRIO_calculos1
       call GIBBS_MACH1D_reativo5
       call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modelo_q, itimax, itemax, tol_e, &
                                            tol_n, T(i), (p(i)-dpY*p(i)/2.0d0)/101325.0d0, OF )
       call GIBBS_EQUILIBRIO_calculos1
       call GIBBS_MACH1D_reativo6
       call GIBBS_MACH1D_reativo4 ( dpY*p(i), gamaf(i), 8314.51d0/Rg(i), p(i), gama(i) )
    end do
    
  end subroutine NUM_EQ_gama_efetivo

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

  subroutine NUM_EQ_grava_canal ( n, unit, para_canal_nome, Cd,  F_star, qlt, &
                                            T_wall,          Taw, sn,     h,   T )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: unit ! unidade de leitura do arquivo de dados
 
    character*(*),intent(in) :: para_canal_nome ! nome do arquivo de saída para o programa CANAL

    real*8,intent(in) :: Cd     ! coeficiente de descarga (adimensional)
    real*8,intent(in) :: F_star ! empuxo dinâmico (adimensional)

    real*8,dimension(:),intent(in) :: qlt    ! fluxo de calor total para a parede (W/m2)
    real*8,dimension(:),intent(in) :: T_wall ! temperatura da parede da tubeira (K)
    real*8,dimension(:),intent(in) :: Taw    ! temperatura de parede adiabática (K)
    real*8,dimension(:),intent(in) :: sn     ! área da parede da tubeira
    real*8,dimension(:),intent(in) :: h      ! coeficiente de convecção (W/m2.K)
    real*8,dimension(:),intent(in) :: T      ! temperatura (K)

    real*8,dimension(:),allocatable :: qt ! taxa de transferência de calor total para a parede (W)

    real*8 :: epsig

    allocate ( qt(n) )

    epsig = epsilon_cte * sigma

    qt = qlt * sn

    open ( unit, file = trim(para_canal_nome), form='unformatted' )

    write ( unit ) n, epsig, Cd, F_star, Taw, T, T_wall, h, sn, qt

    close ( unit )

    deallocate ( qt )

  end subroutine NUM_EQ_grava_canal 

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

  subroutine NUM_EQ_grava_reinicio ( n, unit, nome_ci, u,  p,  T,  ro,  um, rom, pl )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: unit ! unidade de leitura do arquivo de dados
 
    character*(*),intent(in) :: nome_ci ! nome do arquivo de dados de reinício

    real*8,dimension(:),intent(in) :: u   ! velocidade (m/s)
    real*8,dimension(:),intent(in) :: p   ! pressão (Pa)
    real*8,dimension(:),intent(in) :: T   ! temperatura (K)
    real*8,dimension(:),intent(in) :: ro  ! massa específica (kg/m3)
    real*8,dimension(:),intent(in) :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in) :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in) :: pl  ! variação de pressão (Pa)

    open ( unit, file = trim(nome_ci), form='unformatted' )

    write ( unit ) n, u_in, p_in, T_in, pl_in, p_ia, u, p, T, ro, um, rom, pl

    close ( unit )

  end subroutine NUM_EQ_grava_reinicio 

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

  subroutine NUM_EQ_h ( n, fm, rc, Sg, sp, M, gama, cp, mi, T_wall, Pr, h )

    integer,intent(in) :: n  ! número de volumes de controle com dois fictícios

    real*8, intent(in) :: fm ! fluxo de massa numérico (kg/s)
    real*8, intent(in) :: rc ! raio de curvatura na garganta da tubeira (m)
    real*8, intent(in) :: Sg ! área da garganta (m2)

    real*8,dimension(:),intent(in) :: sp     ! área no centro do volume de controle
    real*8,dimension(:),intent(in) :: M      ! número de Mach (adimensional)
	real*8,dimension(:),intent(in) :: gama   ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in) :: cp     ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),intent(in) :: mi     ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),intent(in) :: T_wall ! temperatura da parede da tubeira (K)
    real*8,dimension(:),intent(in) :: Pr     ! número de Prandtl (adim)

    real*8,dimension(:),intent(inout) :: h   ! coeficiente de convecção (W/m2.K)

    real*8 :: Dg ! diâmetro da garganta da tubeira (m)
    real*8 :: pi ! constante matemática
    real*8 :: cte1, cte2, cte3, cte4, cte5, cte6 ! auxiliares

    select case ( tipo_h )

       case ( 1 )
          h = h_cte

       case ( 2 )
          pi = dacos(-1.0d0)
          Dg = 4.0d0 * Sg / pi
          cte1 = 0.026d0*((mi(1)**0.2d0)*cp(1)/(Pr(1)**0.6d0))/(Dg**0.2d0)
          cte2 = ((fm/Sg)**0.8d0)*((Dg/rc)**0.1d0)
          h(1) = 0.0d0
          do i = 2, n-1
             cte3 = 1.0d0 + 0.5d0*(gama(i)-1.0d0)*(M(i)**2)
             cte4 = (0.5d0*(1.0d0+T_wall(i)*cte3/T_cam))**0.68d0
             cte5 = cte3**0.12d0
             cte6 = 1.0d0 / ( cte4 * cte5 )
             h(i) = cte1*cte2*((Sg/sp(i))**0.9d0)*cte6
          end do
          h(n) = 0.0d0

    end select
    
  end subroutine NUM_EQ_h

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

  subroutine NUM_EQ_le_dados ( unit )

    integer,intent(in) :: unit ! unidade de leitura do arquivo de dados

    open ( unit, file = trim(nome_num_var) )

    read ( unit, * ) p_cam
    read ( unit, * ) T_cam

    read ( unit, * ) Rg_cte
    read ( unit, * ) gama_cte

    read ( unit, * ) tipo_cp

    read ( unit, * ) tipo_f
    read ( unit, * ) f_cte
    read ( unit, * ) rugoso

    read ( unit, * ) tipo_h
    read ( unit, * ) h_cte

    read ( unit, * ) tipo_mi
    read ( unit, * ) mi_cte

    read ( unit, * ) tipo_k
    read ( unit, * ) k_cte

    read ( unit, * ) tipo_T_wall
    read ( unit, * ) T_wall_cte
    read ( unit, * ) T_wall_nome

    read ( unit, * ) tipo_fr
    read ( unit, * ) fr_cte

    read ( unit, * ) epsilon_cte
    read ( unit, * ) sigma

    read ( unit, * ) OF
    read ( unit, * ) itimax
    read ( unit, * ) tol_e
    read ( unit, * ) itemax
    read ( unit, * ) tol_n
    read ( unit, * ) modelo_q
    read ( unit, * ) dpY
    read ( unit, * ) it_q
    read ( unit, * ) w_q

    read ( unit, * ) RHG
    read ( unit, * ) para_canal_nome

    read ( unit, * ) calor

    close ( unit )

    cp_cte = gama_cte * Rg_cte / ( gama_cte - 1.0d0 )

  end subroutine NUM_EQ_le_dados 

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

  subroutine NUM_EQ_le_inicio ( unit, nome )

    integer,      intent(in) :: unit ! unidade do arquivo
    character*(*),intent(in) :: nome ! nome do arquivo de dados de configuração do módulo

    open ( unit, file = trim(nome) )

    read ( unit, * ) nome_num_var
    read ( unit, * ) lista
    read ( unit, * ) le
    read ( unit, * ) copia

    close ( unit )

    if ( lista == 1 ) call ARQUIVO_notepad_abre_dado ( nome_num_var )

    if ( le    == 1 ) call NUM_EQ_le_dados ( unit )

    if ( copia == 1 ) call ARQUIVO_copia_arquivo ( nome_num_var )

  end subroutine NUM_EQ_le_inicio

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

  subroutine NUM_EQ_le_reinicio ( n, unit, nome_ci, u,  p,  T,  ro,  um, rom, pl )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: unit ! unidade de leitura do arquivo de dados
 
    character*(*),intent(in) :: nome_ci ! nome do arquivo de dados de reinício

    real*8,dimension(:),intent(out) :: u   ! velocidade (m/s)
    real*8,dimension(:),intent(out) :: p   ! pressão (Pa)
    real*8,dimension(:),intent(out) :: T   ! temperatura (K)
    real*8,dimension(:),intent(out) :: ro  ! massa específica (kg/m3)
    real*8,dimension(:),intent(out) :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(out) :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(out) :: pl  ! variação de pressão (Pa)

    integer :: n_le

    open ( unit, file = trim(nome_ci), form='unformatted' )

    read ( unit ) n_le, u_in, p_in, T_in, pl_in, p_ia, u, p, T, ro, um, rom, pl

    if ( n_le /= n ) stop

    close ( unit )

  end subroutine NUM_EQ_le_reinicio 

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

  subroutine NUM_EQ_le_rhg ( unit )

    integer,intent(in) :: unit ! unidade de leitura do arquivo de dados

    open ( unit, file = 'rhg.itg' )

    read ( unit, * ) itg
    read ( unit, * ) itgmax

    close ( unit )

  end subroutine NUM_EQ_le_rhg 

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

  subroutine NUM_EQ_Mach ( Rg, gama, u, T, M )

    real*8,dimension(:),intent(in)  :: Rg   ! constante do gás (J/kg.K)
    real*8,dimension(:),intent(in)  :: gama ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in)  :: u    ! velocidade (m/s)
    real*8,dimension(:),intent(in)  :: T    ! temperatura (K)

    real*8,dimension(:),intent(out) :: M    ! número de Mach (adimensional)

    M = u / dsqrt ( gama * Rg * T )
  
  end subroutine NUM_EQ_Mach

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

  subroutine NUM_EQ_massa_especifica_media ( n, beta, ro, rom )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: beta ! 0=UDS; 1=CDS

    real*8,dimension(:),intent(in)  :: ro  ! massa específica (kg/m3)

    real*8,dimension(:),intent(out) :: rom ! massa específica média (kg/m3)

    do i = 1, n-1

       rom(i) = ro(i) + beta * 0.5d0 * ( ro(i+1) - ro(i) )

    end do

  end subroutine NUM_EQ_massa_especifica_media

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

  subroutine NUM_EQ_mi_k_congelados ( n, T, mi, k )

    integer,intent(in) :: n ! número de volumes de controle com dois fictícios

    real*8,dimension(:),intent(in)    :: T  ! temperatura (K)

    real*8,dimension(:),intent(inout) :: mi ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),intent(inout) :: k  ! condutividade térmica (W/m.K)

    if ( tipo_mi == 2 .or. tipo_k  == 2 ) then
          
       do i = 1, n

          call GIBBS_MACH1D_recebe2 ( T(i) )

          call GIBBS_CONGELADO_calculos2

          call GIBBS_CONGELADO_envia2 ( mi(i), k(i) )

       end do

    end if

    if ( tipo_mi == 1 ) mi = mi_cte

    if ( tipo_k  == 1 ) k  = k_cte
   
  end subroutine NUM_EQ_mi_k_congelados

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

  subroutine NUM_EQ_pressao_entrada ( gama_cte, Rg_cte, p_cam, u_in, T_in, p_in )

    real*8,intent(in)  :: gama_cte ! razão de calores específicos (adimensional)
    real*8,intent(in)  :: Rg_cte   ! constante do gás (J/kg.K)
    real*8,intent(in)  :: p_cam    ! pressão na câmara de combustão (Pa)
    real*8,intent(in)  :: u_in     ! velocidade na entrada da tubeira (m/s)
    real*8,intent(in)  :: T_in     ! temperatura na entrada da tubeira (K)

    real*8,intent(out) :: p_in  ! pressão na entrada da tubeira (Pa)

    real*8 :: k    ! auxiliar
    real*8 :: M_in ! número de Mach na entrada

    M_in = u_in / dsqrt ( gama_cte * Rg_cte * T_in )

    k = 1.0d0 + 0.5d0 * ( gama_cte - 1.0d0 ) * ( M_in ** 2 )

    p_in = p_cam / ( k ** ( gama_cte / ( gama_cte - 1.0d0 ) ) )

  end subroutine NUM_EQ_pressao_entrada

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

  subroutine NUM_EQ_p_tot ( n, gama, M, p, p_tot )

    integer,intent(in) :: n ! número de pontos a calcular

    real*8,dimension(:),intent(in)  :: gama  ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in)  :: M     ! número de Mach (adimensional)
    real*8,dimension(:),intent(in)  :: p     ! pressão (Pa)

    real*8,dimension(:),intent(out) :: p_tot ! pressão de estagnação (Pa)

    integer :: i  ! número do ponto

    real*8 :: cte ! variável auxiliar
    real*8 :: exp ! variável auxiliar

    do i = 1, n

       exp = gama(i) / (gama(i) - 1.0d0)

       cte = 1.0d0 + (gama(i) - 1.0d0) * (M(i)**2) / 2.0d0
  
       p_tot(i) = p(i) * (cte ** exp)

    end do

  end subroutine NUM_EQ_p_tot

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

  subroutine NUM_EQ_simplec ( n, sp, aw, ap, ae, dp, de )

    integer,intent(in) :: n ! número de volumes de controle com dois fictícios

    real*8,dimension(:),intent(in)  :: sp ! área no centro do volume de controle
    real*8,dimension(:),intent(in)  :: aw ! coeficiente oeste
    real*8,dimension(:),intent(in)  :: ap ! coeficiente central
    real*8,dimension(:),intent(in)  :: ae ! coeficiente leste

    real*8,dimension(:),intent(out) :: dp ! coeficiente SIMPLEC no centro do volume
    real*8,dimension(:),intent(out) :: de ! coeficiente SIMPLEC na face leste do volume

    dp(1) = 0.0d0
    dp(n) = 0.0d0

    do i = 2, n-1
       dp(i) = sp(i) / ( ap(i) + aw(i) + ae(i) )
    end do

    de(1)   = dp(2)
    de(n-1) = dp(n-1)

    do i = 2,n-2
       de(i) = 0.5d0 * ( dp(i) + dp(i+1) )
    end do

  end subroutine NUM_EQ_simplec

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

  subroutine NUM_EQ_T_tot ( n, gama, M, T, T_tot )

    integer,intent(in) :: n ! número de pontos a calcular

    real*8,dimension(:),intent(in)  :: gama  ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in)  :: M     ! número de Mach (adimensional)
    real*8,dimension(:),intent(in)  :: T     ! temperatura (K)

    real*8,dimension(:),intent(out) :: T_tot ! temperatura de estagnação (K)

    integer :: i  ! número do ponto

    real*8 :: cte ! variável auxiliar

    do i = 1, n

       cte = 1.0d0 + (gama(i) - 1.0d0) * (M(i)**2) / 2.0d0
  
       T_tot(i) = T(i) * cte
  
    end do

  end subroutine NUM_EQ_T_tot

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

  subroutine NUM_EQ_T_wall ( unit, n, T_wall )

    integer,intent(in) :: unit ! unidade do arquivo
    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios

    real*8,dimension(:),intent(out) :: T_wall ! temperatura da parede da tubeira (K)

    integer :: n_le

    select case ( tipo_T_wall )

       case ( 1 )
          T_wall = T_wall_cte

       case ( 2 )
          open ( unit, file = trim(T_wall_nome), form = 'unformatted' )
          read ( unit ) n_le, T_wall
          close ( unit )
          if ( n_le /= n ) then
             stop
          end if
!          call ARQUIVO_copia_arquivo ( T_wall_nome )

    end select

  end subroutine NUM_EQ_T_wall

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

  subroutine NUM_EQ_uds_coef_pressao ( n,  dt, Rg, pl_in, ro, rom, roa, um, T, de, &
                                       se, sp, xe, xp,    pl, aw,  ap,  ae, bp )

    integer,intent(in) :: n      ! número de volumes de controle com dois fictícios
    real*8, intent(in) :: dt     ! intervalo de tempo (s)
    real*8, intent(in) :: pl_in  ! variação da pressão na entrada da tubeira (Pa)

    real*8,dimension(:),intent(in)  :: ro  ! massa específica (kg/m3)
    real*8,dimension(:),intent(in)  :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in)  :: roa ! massa específica anterior (kg/m3)
    real*8,dimension(:),intent(in)  :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in)  :: T   ! temperatura (K)
    real*8,dimension(:),intent(in)  :: de  ! coeficiente SIMPLEC na face leste do volume
    real*8,dimension(:),intent(in)  :: se  ! área na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: sp  ! área no centro do volume de controle
    real*8,dimension(:),intent(in)  :: xe  ! coordenada espacial na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: xp  ! coordenada espacial no centro do volume de controle
    real*8,dimension(:),intent(in)  :: pl  ! variação de pressão (Pa)
    real*8,dimension(:),intent(in)  :: Rg  ! constante do gás (J/kg.K)

    real*8,dimension(:),intent(out) :: aw  ! coeficiente oeste
    real*8,dimension(:),intent(out) :: ap  ! coeficiente central
    real*8,dimension(:),intent(out) :: ae  ! coeficiente leste
    real*8,dimension(:),intent(out) :: bp  ! termo independente

    real*8 :: dx    ! tamanho do volume de controle elementar
    real*8 :: fator ! auxiliar

    ! volume 1 (fictício)

    aw(1) = 0.0d0
    ap(1) = 1.0d0
    ae(1) = 1.0d0
    bp(1) = 2.0d0 * pl_in

    ! volumes internos

    do i = 2, n-1

       dx = xe(i) - xe(i-1)

       aw(i) = - rom(i-1) * de(i-1) * se(i-1)          &
               - um(i-1)  * se(i-1) / ( Rg(i-1) * T(i-1) )

       ae(i) = - rom(i) * de(i) * se(i)

       ap(i) = ( sp(i) * dx / dt + um(i) * se(i) ) / ( Rg(i) * T(i) )   &
             + rom(i-1) * de(i-1) * se(i-1)                              &
             + rom(i)   * de(i)   * se(i)

       bp(i) = - ( ( ro(i) - roa(i) ) * sp(i) * dx / dt              &
             + ro(i) * um(i) * se(i) - ro(i-1) * um(i-1) * se(i-1) )

    end do

    ! volume n (fictício)

    fator = 2 * ( xp(n) - xp(n-1) ) / ( xp(n-1) - xp(n-2) )
    aw(n) = -1.0d0
    ap(n) =  1.0d0
    ae(n) =  0.0d0
    bp(n) =  fator * ( pl(n-1) - pl(n-2) )

  end subroutine NUM_EQ_uds_coef_pressao

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

  subroutine NUM_EQ_uds_coef_temperatura ( n,  Ne, dt, T_in, epsilon_cte, sigma,   &
                                           ro, roa,   rom,    u,  um,      T,      &
                                           Ta, p,     pa,     se, sp,      rp,     &
                                           xe, xp,    sn,     f,  h,       T_wall, &
                                           cp, C_rec, aw, ap, ae,  bp, Yi, hi )

    integer,intent(in) :: n           ! número de volumes de controle com dois fictícios
    integer,intent(in) :: Ne          ! número de espécies químicas

    real*8, intent(in) :: dt          ! intervalo de tempo (s)
    real*8, intent(in) :: T_in        ! temperatura na entrada da tubeira (K)
    real*8, intent(in) :: epsilon_cte ! emissividade média parede-gases (adim)
    real*8, intent(in) :: sigma       ! constante de Stefan-Boltzmann (W/m2.K4)

    real*8,dimension(:),intent(in)  :: ro     ! massa específica (kg/m3)
    real*8,dimension(:),intent(in)  :: roa    ! massa específica anterior (kg/m3)
    real*8,dimension(:),intent(in)  :: rom    ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in)  :: u      ! velocidade (m/s)
    real*8,dimension(:),intent(in)  :: um     ! velocidade média (m/s)
    real*8,dimension(:),intent(in)  :: T      ! temperatura (K)
    real*8,dimension(:),intent(in)  :: Ta     ! temperatura anterior (K)
    real*8,dimension(:),intent(in)  :: p      ! pressão (Pa)
    real*8,dimension(:),intent(in)  :: pa     ! pressão anterior (Pa)
    real*8,dimension(:),intent(in)  :: se     ! área na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: sp     ! área no centro do volume de controle
    real*8,dimension(:),intent(in)  :: rp     ! raio no centro do volume de controle
    real*8,dimension(:),intent(in)  :: xe     ! coordenada espacial na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: xp     ! coordenada espacial no centro do volume de controle
    real*8,dimension(:),intent(in)  :: sn     ! área da parede da tubeira
    real*8,dimension(:),intent(in)  :: f      ! fator de atrito de Darcy (adim)
    real*8,dimension(:),intent(in)  :: h      ! coeficiente de convecção (W/m2.K)
    real*8,dimension(:),intent(in)  :: T_wall ! temperatura da parede da tubeira (K)
    real*8,dimension(:),intent(in)  :: cp     ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),intent(in)  :: C_rec  ! coeficiente de recuperação de Taw

    real*8,dimension(:),intent(out) :: aw     ! coeficiente oeste
    real*8,dimension(:),intent(out) :: ap     ! coeficiente central
    real*8,dimension(:),intent(out) :: ae     ! coeficiente leste
    real*8,dimension(:),intent(out) :: bp     ! termo independente

    real*8,dimension(:,:),intent(in) :: Yi  ! fração mássica (adim.)
    real*8,dimension(:,:),intent(in) :: hi  ! entalpia mássica (J/kg)

    real*8  :: pi    ! constante matemática
    real*8  :: dx    ! tamanho do volume de controle elementar
    real*8  :: fator ! auxiliar

    integer :: ii    ! número da espécie
    real*8  :: br, fmw ! auxiliares reativos

    pi = dacos(-1.0d0)

    ! volume 1 (fictício)

    aw(1) = 0.0d0
    ap(1) = 1.0d0
    ae(1) = 1.0d0
    bp(1) = 2.0d0 * T_in

    ! volumes internos

    do i = 2, n-1

       dx = xe(i) - xe(i-1)

       aw(i) = - cp(i) * rom(i-1) * um(i-1) * se(i-1)

       ae(i) = 0.0d0

       ap(i) = cp(i) * roa(i) * sp(i) * dx / dt - aw(i)  &
             + h(i) * C_rec(i) * sn(i)                   &
             + epsilon_cte * sigma * sn(i) * ( T(i) ** 3 )    

       br  = 0.0d0
       if ( itt >= it_q .and. mod(itt,w_q)==0 ) then
          fmw = rom(i-1) * um(i-1) * se(i-1)
          do ii = 1, Ne
             br = br - hi(i,ii) * ( Yi(i,ii) - Yi(i-1,ii) )
          end do
          br = br * fmw
       end if

       bp(i) = cp(i) * roa(i) * sp(i) * dx * Ta(i) / dt        &
             + sp(i) * dx * ( p(i) - pa(i) ) / dt              &
             + sp(i) * u(i) * ( p(i+1) - p(i-1) ) * 0.5d0      &
             + pi * f(i) * ro(i) * (u(i)**3) * rp(i) * dx / 4  &
             + h(i) * T_wall(i) * sn(i)                        &
             + epsilon_cte * sigma * sn(i) * ( T_wall(i) ** 4 ) + br

    end do

    ! volume n (fictício)

    fator = 2 * ( xp(n) - xp(n-1) ) / ( xp(n-1) - xp(n-2) )
    aw(n) = -1.0d0
    ap(n) =  1.0d0
    ae(n) =  0.0d0
    bp(n) =  fator * ( T(n-1) - T(n-2) )

  end subroutine NUM_EQ_uds_coef_temperatura

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

  subroutine NUM_EQ_uds_coef_velocidade ( n,  dt, f,  ro, roa, u,  ua, um, p,  se, &
                                          sp, rp, xe, xp, rom, aw, ap, ae, bp, bf  )

    integer,intent(in) :: n  ! número de volumes de controle com dois fictícios
    real*8, intent(in) :: dt ! intervalo de tempo (s)

    real*8,dimension(:),intent(in)  :: ro  ! massa específica (kg/m3)
    real*8,dimension(:),intent(in)  :: roa ! massa específica anterior (kg/m3)
    real*8,dimension(:),intent(in)  :: u   ! velocidade (m/s)
    real*8,dimension(:),intent(in)  :: ua  ! velocidade anterior (m/s)
    real*8,dimension(:),intent(in)  :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in)  :: p   ! pressão (Pa)
    real*8,dimension(:),intent(in)  :: se  ! área na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: sp  ! área no centro do volume de controle
    real*8,dimension(:),intent(in)  :: rp  ! raio no centro do volume de controle
    real*8,dimension(:),intent(in)  :: xe  ! coordenada espacial na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: xp  ! coordenada espacial no centro do volume de controle
    real*8,dimension(:),intent(in)  :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in)  :: f   ! fator de atrito de Darcy (adim)

    real*8,dimension(:),intent(out) :: aw  ! coeficiente oeste
    real*8,dimension(:),intent(out) :: ap  ! coeficiente central
    real*8,dimension(:),intent(out) :: ae  ! coeficiente leste
    real*8,dimension(:),intent(out) :: bp  ! termo independente
    real*8,dimension(:),intent(out) :: bf  ! termo de fricção para um

    real*8 :: pi    ! constante matemática
    real*8 :: dx    ! tamanho do volume de controle elementar
    real*8 :: fator ! auxiliar

    pi = dacos(-1.0d0)

    ! volume 1 (fictício)

    fator = 2 * xp(2) / ( xp(3) - xp(2) )
    aw(1) =  0.0d0
    ap(1) =  1.0d0
    ae(1) = -1.0d0
    bp(1) = -fator * ( u(3) - u(2) )
    bf(1) = 0.0d0

    ! volumes internos

    do i = 2, n-1

       dx = xe(i) - xe(i-1)

       aw(i) = - rom(i-1) * um(i-1) * se(i-1)

       ae(i) = 0.0d0

       ap(i) = roa(i) * sp(i) * dx / dt - ( aw(i) + ae(i) )

       bf(i) = - pi * f(i) * ro(i) * (u(i)**2) * rp(i) * dx / 4   

       bp(i) = roa(i) * sp(i) * dx * ua(i) / dt + bf(i)   &
             + 0.5d0  * sp(i) * ( p(i-1) - p(i+1) )

    end do

    ! volume n (fictício)

    fator = 2 * ( xp(n) - xp(n-1) ) / ( xp(n-1) - xp(n-2) )
    aw(n) = -1.0d0
    ap(n) =  1.0d0
    ae(n) =  0.0d0
    bp(n) =  fator * ( u(n-1) - u(n-2) )
    bf(n) =  0.0d0

  end subroutine NUM_EQ_uds_coef_velocidade

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

  subroutine NUM_EQ_velocidade_media ( n,  dt, roa, u,  uma, p,  se, &
                                       sp, xe, aw,  ap, ae,  bc, bf, um )

    integer,intent(in) :: n  ! número de volumes de controle com dois fictícios
    real*8, intent(in) :: dt ! intervalo de tempo (s)

    real*8,dimension(:),intent(in)    :: roa ! massa específica anterior (kg/m3)
    real*8,dimension(:),intent(in)    :: u   ! velocidade (m/s)
    real*8,dimension(:),intent(in)    :: uma ! velocidade média anterior (m/s)
    real*8,dimension(:),intent(in)    :: p   ! pressão (Pa)
    real*8,dimension(:),intent(in)    :: se  ! área na face leste do volume de controle
    real*8,dimension(:),intent(in)    :: sp  ! área no centro do volume de controle
    real*8,dimension(:),intent(in)    :: xe  ! coordenada espacial na face leste do volume de controle
    real*8,dimension(:),intent(in)    :: aw  ! coeficiente oeste
    real*8,dimension(:),intent(in)    :: ap  ! coeficiente central
    real*8,dimension(:),intent(in)    :: ae  ! coeficiente leste
    real*8,dimension(:),intent(in)    :: bc  ! termo da correção adiada
    real*8,dimension(:),intent(in)    :: bf  ! termo de fricção para um

    real*8,dimension(:),intent(inout) :: um  ! velocidade média (m/s)

    real*8 :: massa_e ! massa do volume E
    real*8 :: massa_p ! massa do volume P

    real*8 :: somap, somae

    ! média antiga

    do i = 2, n-2

       massa_p = roa(i) * sp(i) * (xe(i) - xe(i-1))

       massa_e = roa(i+1) * sp(i+1) * (xe(i+1) - xe(i))

       somap = aw(i)*u(i-1) + ae(i)*u(i+1)
 
       somae = aw(i+1)*u(i) + ae(i+1)*u(i+2)

       um(i) = (-somap - somae + bc(i) + bc(i+1) + bf(i) + bf(i+1)       &
             + (massa_p+massa_e)*uma(i)/dt - 2.0d0*se(i)*(p(i+1)-p(i)))  &
             / (ap(i)+ap(i+1))

    end do

    um(1) = 0.5d0 * ( u(1) + u(2) )

    um(n-1) = 0.5d0 * ( u(n) + u(n-1) )

  end subroutine NUM_EQ_velocidade_media

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

end module NUMERICO_EQUILIBRIO