module NUMERICO_CONSTANTE

! sigla = NUM_CTE_

! objetivo: calcular a solução numérica do modelo matemático
!           para propriedades e dados constantes

! versão original 1.0 (03 Jun 99)

! programador: Carlos Henrique Marchi
! local:       SINMEC/UFSC - Florianópolis - SC

! versão atual    1.1 (16 Dez 99)

! última alteração =   11 Abr 00

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

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

  ! NUM_CTE_cds_coef_pressao
  ! NUM_CTE_cds_coef_temperatura
  ! NUM_CTE_cds_coef_velocidade
  ! NUM_CTE_correcao_com_plinha
  ! NUM_CTE_escreve_campos_auxiliares
  ! NUM_CTE_executa
  ! NUM_CTE_le_dados
  ! NUM_CTE_le_inicio
  ! NUM_CTE_massa_especifica_media
  ! NUM_CTE_passa_dados
  ! NUM_CTE_pressao_entrada
  ! NUM_CTE_simplec
  ! NUM_CTE_uds_coef_pressao
  ! NUM_CTE_uds_coef_temperatura
  ! NUM_CTE_uds_coef_velocidade
  ! NUM_CTE_velocidade_media

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

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

  ! NUM_CTE_le_dados
  ! NUM_CTE_le_inicio

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

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

  ! NUM_CTE_executa

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

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

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

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

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      ! constante do gás (J/kg.K)
real*8,private :: cp      ! calor específico à pressão constante (J/kg.K)
real*8,private :: f       ! fator de atrito de Darcy (adim)
real*8,private :: h       ! coeficiente de convecção (W/m2.K)
real*8,private :: T_wall  ! temperatura da parede da tubeira (K)
real*8,private :: fr      ! fator de recuperação (adim)
real*8,private :: epsilon ! emissividade média parede-gases (adim)
real*8,private :: sigma   ! constante de Stefan-Boltzmann (W/m2.K4)
real*8,private :: gama    ! razão de calores específicos (adimensional)

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_CONSTANTE
integer,private :: le    ! lê conteúdo do arquivo de dados do módulo NUMERICO_CONSTANTE
integer,private :: copia ! gera back up do arquivo de dados do módulo NUMERICO_CONSTANTE

character*50,private :: nome_num_cte ! nome do arquivo de dados do módulo NUMERICO_CONSTANTE

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

contains

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

  subroutine NUM_CTE_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_CTE_cds_coef_pressao

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

  subroutine NUM_CTE_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, 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 * ( oeste - leste )

    end do

    ! volume n (fictício)

    bc(n) = 0.0d0

    ! atualiza o termo independente

    bp = bp + bc

  end subroutine NUM_CTE_cds_coef_temperatura

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

  subroutine NUM_CTE_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_CTE_cds_coef_velocidade

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

  subroutine NUM_CTE_correcao_com_plinha ( n, xp, dp, de, pl, 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(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_CTE_correcao_com_plinha

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

  subroutine NUM_CTE_escreve_campos_auxiliares ( unit, n,   w,   h,   epsilon, xp, xe, &
                                                 Taw,  qlr, qlt, pl,  um,      rom     )

    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

    real*8,intent(in) :: h       ! coeficiente de convecção (W/m2.K)
    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) :: Taw ! temperatura de parede adiabática (K)
    real*8,dimension(:),intent(in) :: qlt ! fluxo de calor total para a parede (W/m2)
    real*8,dimension(:),intent(in) :: qlr ! fluxo de calor radiativo para a parede (W/m2)
    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)

    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
                            '*** 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 > 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                ', & ! 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
                            '*** plinha ***                                     ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'pl                   ', & ! 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                  ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xe,                      & ! coordenada espacial
                            rom                      ) ! variável a escrever

    close ( unit )

  end subroutine NUM_CTE_escreve_campos_auxiliares

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

  subroutine NUM_CTE_executa ( n,     w_it,  beta, nitm, dt, fm_exato, F_exato, &
                               xe,    xp,    rp,   se,   sp, sn,       fm,      &
                               F_ex,  u,     p,    T,    ro, r_qwh,    q_sai,   &
                               w_aux, lista_aux    )

    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) :: w_aux ! freqüência de escrita de campos auxiliares
    integer,intent(in) :: lista_aux ! lista arquivo com os campos auxiliares (1); 0=não

    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,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)

    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 arquivo do caso a simular
    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)

    ! 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) )

    ! inicialização de campos
    um    = u(1:n-1)
    rom   = ro(1:n-1)
    pl    = 0.0d0
    pl_in = 0.0d0
    bc    = 0.0d0

    ! inicialização na entrada da tubeira
    u_in = um(1)
    T_in = T_cam - 0.5d0*(gama-1.0d0)*(u_in**2)/(gama*Rg)
    call NUM_CTE_pressao_entrada ( gama, Rg, 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 * 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 * T(n) )

    ! cria o arquivo das iterações
    if ( w_it /= 0 ) then
       call ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.mach1d.itn'
       call ARQUIVO_cria_arquivo ( unit_it, comando )
       write ( unit_it, 1 ) 
       1 format (t1,'itt, cd, F_star, T_ex, u_in, u_ex')
    end if

    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_CTE_pressao_entrada ( gama, Rg, p_cam, u_in, T_in, p_in )
       p(1) = 2.0d0*p_in - p(2)
       pl_in = p_in - p_ia

       ! solução da velocidade nodal
       call NUM_CTE_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_CTE_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_CTE_simplec ( n, sp, aw, ap, ae, dp, de )

       ! cálculo da velocidade média
       call NUM_CTE_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-1.0d0)*(u_in**2)/(gama*Rg)
       call NUM_CTE_uds_coef_temperatura ( n,   dt,      T_in,  f,   h,      gama, Rg,  &
                                           fr,  epsilon, sigma, cp,  T_wall, ro,   roa, &
                                           rom, u,       um,    T,   Ta,     p,    pa,  &
                                           se,  sp,      rp,    xe,  xp,     sn,   aw,  &
                                           ap,  ae,      bp,    C_rec        )
       if ( beta > 0 ) &
          call NUM_CTE_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_CTE_massa_especifica_media ( n, beta, ro, rom )

       ! solução da correção da pressão
       call NUM_CTE_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_CTE_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_CTE_correcao_com_plinha ( n, xp, dp, de, pl, se, u, p, T, ro, um )
       call NUM_CTE_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 ( 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)
       2 format ( i10, 5 ( 1pe23.15 ) )

    end do

    call ARQUIVO_cronometro ( tcpu2 )

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

    ! entrada da tubeira
    u(1)  = u_in
    p(1)  = p_in
    T(1)  = T_in
    ro_in = p_in / ( Rg * T_in )
    ro(1) = ro_in

    ! saída da tubeira
    T_ex  = 0.5d0 * ( T(n-1) + T(n) )
    p_ex  = 0.5d0 * ( p(n-1) + p(n) )
    ro_ex = p_ex / ( Rg * T_ex )
    u(n)  = um(n-1)
    p(n)  = p_ex
    T(n)  = T_ex
    ro(n) = ro_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 ) 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 ***')
    end if

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

    if ( w_aux /= 0 ) &
    call NUM_CTE_escreve_campos_auxiliares ( unit_au, n,   w_aux, h,  epsilon, xp, xe, &
                                             Taw,     qlr, qlt,   pl, um,      rom     )

    q_sai = fm * cp * T_ex

    ! 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
  
    deallocate ( C_rec, Taw, qlt, qlr, pl, um, rom )

  end subroutine NUM_CTE_executa

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

  subroutine NUM_CTE_le_dados ( unit )

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

    open ( unit, file = trim(nome_num_cte) )

    read ( unit, * ) p_cam
    read ( unit, * ) T_cam
    read ( unit, * ) Rg
    read ( unit, * ) gama
    read ( unit, * ) f
    read ( unit, * ) h
    read ( unit, * ) T_wall
    read ( unit, * ) fr
    read ( unit, * ) epsilon
    read ( unit, * ) sigma

    close ( unit )

    cp = gama * Rg / ( gama - 1.0d0 )

  end subroutine NUM_CTE_le_dados 

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

  subroutine NUM_CTE_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_cte
    read ( unit, * ) lista
    read ( unit, * ) le
    read ( unit, * ) copia

    close ( unit )

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

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

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

  end subroutine NUM_CTE_le_inicio

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

  subroutine NUM_CTE_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_CTE_massa_especifica_media

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

  subroutine NUM_CTE_passa_dados ( gama_r, Rg_r, p_cam_r, T_cam_r )

    real*8, intent(out) :: gama_r    ! razão de calores específicos (adimensional)
    real*8, intent(out) :: Rg_r      ! constante do gás (J/kg.K)
    real*8, intent(out) :: p_cam_r   ! pressão na câmara de combustão (Pa)
    real*8, intent(out) :: T_cam_r   ! temperatura na câmara de combustão (K)

    gama_r = gama
    Rg_r   = Rg
    p_cam_r = p_cam
    T_cam_r = T_cam

  end subroutine NUM_CTE_passa_dados

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

  subroutine NUM_CTE_pressao_entrada ( gama, Rg, p_cam, u_in, T_in, p_in )

    real*8,intent(in)  :: gama  ! razão de calores específicos (adimensional)
    real*8,intent(in)  :: Rg    ! 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 * Rg * T_in )

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

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

  end subroutine NUM_CTE_pressao_entrada

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

  subroutine NUM_CTE_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_CTE_simplec

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

  subroutine NUM_CTE_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) :: Rg ! constante do gás (J/kg.K)
    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(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 * T(i-1) )

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

       ap(i) = ( sp(i) * dx / dt + um(i) * se(i) ) / ( Rg * 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_CTE_uds_coef_pressao

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

  subroutine NUM_CTE_uds_coef_temperatura ( n,   dt,      T_in,  f,   h,      gama, Rg,  &
                                            fr,  epsilon, sigma, cp,  T_wall, ro,   roa, &
                                            rom, u,       um,    T,   Ta,     p,    pa,  &
                                            se,  sp,      rp,    xe,  xp,     sn,   aw,  &
                                            ap,  ae,      bp,    C_rec        )

    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) :: T_in    ! temperatura na entrada da tubeira (K)
    real*8, intent(in) :: f       ! fator de atrito de Darcy (adim)
    real*8, intent(in) :: h       ! coeficiente de convecção (W/m2.K)
    real*8, intent(in) :: gama    ! razão de calores específicos (adimensional)
    real*8, intent(in) :: Rg      ! constante do gás (J/kg.K)
    real*8, intent(in) :: fr      ! fator de recuperação (adim)
    real*8, intent(in) :: epsilon ! emissividade média parede-gases (adim)
    real*8, intent(in) :: sigma   ! constante de Stefan-Boltzmann (W/m2.K4)
    real*8, intent(in) :: cp      ! calor específico à pressão constante (J/kg.K)
    real*8, intent(in) :: T_wall  ! temperatura da parede da tubeira (K)

    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(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) :: C_rec ! coeficiente de recuperação de Taw

    real*8 :: pi    ! constante matemática
    real*8 :: dx    ! tamanho do volume de controle elementar
    real*8 :: fator ! auxiliar
    real*8 :: M_aux ! número de Mach auxiliar

    pi = dacos(-1.0d0)

    C_rec(1) = 0.0d0
    do i = 2, n-1
       M_aux    = u(i) / dsqrt ( gama * Rg * T(i) )
       C_rec(i) = 1.0d0 + fr*(gama-1.0d0)*(M_aux**2)/2.0d0
    end do
    C_rec(n) = 0.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 * rom(i-1) * um(i-1) * se(i-1)

       ae(i) = 0.0d0

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

       bp(i) = cp * 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 * ro(i) * (u(i)**3) * rp(i) * dx / 4    &
             + h * T_wall * sn(i)                             &
             + epsilon * sigma * sn(i) * ( T_wall ** 4 )

    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_CTE_uds_coef_temperatura

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

  subroutine NUM_CTE_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, intent(in) :: f  ! fator de atrito de Darcy (adim)

    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(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 * 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_CTE_uds_coef_velocidade

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

  subroutine NUM_CTE_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_CTE_velocidade_media

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

end module NUMERICO_CONSTANTE