module NUMERICO_VARIAVEL

! sigla = NUM_VAR_

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

! versão original 1.0 (20 Jan 00)

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

! versão atual    1.0 (20 Jan 00)

! última alteração =    3 Mai 00

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

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

  ! NUM_VAR_cds_coef_pressao
  ! NUM_VAR_cds_coef_temperatura
  ! NUM_VAR_cds_coef_velocidade
  ! NUM_VAR_correcao_com_plinha
  ! NUM_VAR_cp_gama
  ! NUM_VAR_cp_le_coeficientes
  ! NUM_VAR_escreve_campos_auxiliares
  ! NUM_VAR_executa
  ! NUM_VAR_f
  ! NUM_VAR_fr_C_rec
  ! NUM_VAR_grava_canal
  ! NUM_VAR_grava_reinicio
  ! NUM_VAR_h
  ! NUM_VAR_k
  ! NUM_VAR_k_le_coeficientes
  ! NUM_VAR_le_dados
  ! NUM_VAR_le_inicio
  ! NUM_VAR_le_reinicio
  ! NUM_VAR_le_rhg
  ! NUM_VAR_Mach
  ! NUM_VAR_massa_especifica_media
  ! NUM_VAR_mi
  ! NUM_VAR_mi_le_coeficientes
  ! NUM_VAR_pressao_entrada
  ! NUM_VAR_p_tot
  ! NUM_VAR_simplec
  ! NUM_VAR_T_tot
  ! NUM_VAR_T_wall
  ! NUM_VAR_uds_coef_pressao
  ! NUM_VAR_uds_coef_temperatura
  ! NUM_VAR_uds_coef_velocidade
  ! NUM_VAR_velocidade_media

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

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

  ! NUM_VAR_cp_le_coeficientes
  ! NUM_VAR_k_le_coeficientes
  ! NUM_VAR_le_dados
  ! NUM_VAR_le_inicio
  ! NUM_VAR_le_reinicio
  ! NUM_VAR_le_rhg
  ! NUM_VAR_mi_le_coeficientes
  ! NUM_VAR_T_wall

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

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

  ! NUM_VAR_executa
  ! NUM_VAR_grava_canal
  ! NUM_VAR_grava_reinicio

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

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_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_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_var    ! nome do arquivo de dados do módulo NUMERICO_CONSTANTE

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

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

contains

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

  subroutine NUM_VAR_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_VAR_cds_coef_pressao

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

  subroutine NUM_VAR_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_VAR_cds_coef_temperatura

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

  subroutine NUM_VAR_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_VAR_cds_coef_velocidade

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

  subroutine NUM_VAR_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_cte * 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_VAR_correcao_com_plinha

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

  subroutine NUM_VAR_cp_gama ( n, T, cp, gama, a_cp )

    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) :: 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(in)  :: a_cp ! coeficientes do calor específico à pressão constante (J/kg.K) 

    integer :: j ! nível de temperatura: 1>=1000K, 2<1000K

    select case ( tipo_cp )

       case ( 1 )
          cp   = cp_cte
          gama = gama_cte

       case ( 2 )
          do i = 1, n
             if ( T(i) >= 1.0d+3) then
                j = 1
             else
                j = 2
             end if
             cp(i) = (a_cp(1,j) + a_cp(2,j)*T(i) + a_cp(3,j)*(T(i)**2)  &
                   +  a_cp(4,j)*(T(i)**3) + a_cp(5,j)*(T(i)**4)) * Rg_cte
          end do
          gama = cp / ( cp - Rg_cte )

    end select
    
  end subroutine NUM_VAR_cp_gama

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

  subroutine NUM_VAR_cp_le_coeficientes ( unit, a_cp )

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

    real*8,dimension(:,:),intent(out) :: a_cp ! coeficientes do calor específico à pressão constante (J/kg.K) 

    integer :: j ! nível de temperatura: 1>=1000K, 2<1000K
    integer :: l ! número do coeficiente

    open ( unit, file = 'propri1.dat' )

    read(unit,*)
    read(unit,*)   
    read(unit,*)   
    read(unit,*)
	    
    do j = 1, 2
       do l = 1, 5
          read(unit,*) a_cp(l,j)
       end do
       read(unit,*)
       read(unit,*)
       read(unit,*)
    end do

    close(unit)

  end subroutine NUM_VAR_cp_le_coeficientes

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

  subroutine NUM_VAR_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 )

    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) :: 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) :: gama   ! razão de calores específicos (adimensional)
    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)

    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                ', & ! 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
                            '*** calor específico à pressão constante ***       ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'cp (J/kg.K)          ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! 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
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! 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
                            '*** 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                   ', & ! 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_VAR_escreve_campos_auxiliares

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

  subroutine NUM_VAR_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 :: cp     ! calor específico à pressão constante (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)

    real*8,dimension(5,2) :: a_cp ! coeficientes do calor específico à pressão constante
    real*8,dimension(4,2) :: a_mi ! coeficientes da viscosidade absoluta 
    real*8,dimension(4,2) :: a_k  ! coeficientes da condutividade térmica 

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

    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_VAR_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_VAR_le_rhg ( unit_au )

    ! 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 (t1,'itt, cd, F_star, T_ex, u_in, u_ex')
    end if

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

    ! lê arquivo de dados para mi variável
    if ( tipo_mi == 2 ) call NUM_VAR_mi_le_coeficientes ( unit_au, a_mi )

    ! lê arquivo de dados para k variável
    if ( tipo_k == 2 ) call NUM_VAR_k_le_coeficientes ( unit_au, a_k )

    ! lê arquivo de dados para cp variável
    if ( tipo_cp == 2 ) call NUM_VAR_cp_le_coeficientes ( unit_au, a_cp )

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

    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_VAR_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
       !ro(1) = p(1) / ( Rg * T(1) )

       ! cálculo das propriedades
       call NUM_VAR_mi ( n, T, a_mi, mi )
       call NUM_VAR_k  ( n, T, a_k,  k  )
       call NUM_VAR_cp_gama ( n, T, cp, gama, a_cp )
       Red = ro * u * 2 * rp / mi
       Pr  = cp * mi / k
       call NUM_VAR_Mach ( Rg_cte, gama, u, T, M )
       call NUM_VAR_fr_C_rec ( gama, Pr, M, fr, C_rec )
       call NUM_VAR_f ( rp, Red, f )
       call NUM_VAR_h ( n, fm, rc, Sg, sp, M, gama, cp, mi, T_wall, Pr, h )

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

       ! cálculo da velocidade média
       call NUM_VAR_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_VAR_uds_coef_temperatura ( n,  dt,   T_in,  Rg_cte, epsilon_cte, sigma,  &
                                           ro, roa,  rom,   u,      um,          T,      &
                                           Ta, p,    pa,    se,     sp,          rp,     &
                                           xe, xp,   sn,    f,      h,           T_wall, &
                                           cp, gama, C_rec, aw,     ap,          ae,     &
                                           bp  )
       if ( beta > 0 ) &
          call NUM_VAR_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_cte * T )
       call NUM_VAR_massa_especifica_media ( n, beta, ro, rom )

       ! solução da correção da pressão
       call NUM_VAR_uds_coef_pressao ( n,  dt, Rg_cte, pl_in, ro, rom, roa, um, T, de, &
                                       se, sp, xe,     xp,    pl, aw,  ap,  ae, bp     )
       if ( beta > 0 ) &
          call NUM_VAR_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_VAR_correcao_com_plinha ( n, xp, dp, de, pl, se, u, p, T, ro, um )
       call NUM_VAR_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)
             2 format ( i10, 5 ( 1pe23.15 ) )
       end if

    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_VAR_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
    ro_in = p_in / ( Rg_cte * 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_cte * 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 .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 ***')
    end if

    ! 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 ) &
       call NUM_VAR_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 )

    q_sai = fm * cp(n-1) * 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

    call NUM_VAR_Mach ( Rg_cte, gama, u, T, M )
    call NUM_VAR_p_tot ( n, gama, M, p, p_tot )
    call NUM_VAR_T_tot ( n, gama, M, T, T_tot )

    ! gera arquivo de dados para o programa CANAL
    if ( RHG == 1 ) &
       call NUM_VAR_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 )
    deallocate ( gama,  cp,  f,   h,  mi, k, T_wall, fr, Red, Pr )

    itgr    = itg
    itgmaxr = itgmax

  end subroutine NUM_VAR_executa

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

  subroutine NUM_VAR_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_VAR_f

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

  subroutine NUM_VAR_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_VAR_fr_C_rec

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

  subroutine NUM_VAR_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_VAR_grava_canal 

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

  subroutine NUM_VAR_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_VAR_grava_reinicio 

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

  subroutine NUM_VAR_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_VAR_h

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

  subroutine NUM_VAR_k ( n, T, a_k, 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(in) :: a_k ! coeficientes da condutividade térmica 

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

    integer :: j ! nível de temperatura: 1>=1000K, 2<1000K

    select case ( tipo_k )

       case ( 1 )
          k = k_cte

       case ( 2 )
          do i = 1, n
             if ( T(i) >= 1.0d+3) then
                j = 1
             else
                j = 2
             end if
             k(i) = dexp((a_k(1,j)*dlog(T(i)) + a_k(2,j)/T(i)   &
                    + a_k(3,j)/(T(i)**2) + a_k(4,j))) * 1.0d-4
          end do

    end select
   
  end subroutine NUM_VAR_k

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

  subroutine NUM_VAR_k_le_coeficientes ( unit, a_k )

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

    real*8,dimension(:,:),intent(out) :: a_k ! coeficientes da condutividade térmica 

    integer :: j ! nível de temperatura: 1>=1000K, 2<1000K
    integer :: l ! número do coeficiente

    open(unit, file = 'propri2.dat' )

    read(unit,*)
    read(unit,*)   
    read(unit,*)
    read(unit,*)   
	    
    do j=1,2
       read(unit,*) ( a_k(l,j), l=1,4 )
    end do

    close(unit)

  end subroutine NUM_VAR_k_le_coeficientes

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

  subroutine NUM_VAR_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, * ) RHG
    read ( unit, * ) para_canal_nome

    close ( unit )

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

  end subroutine NUM_VAR_le_dados 

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

  subroutine NUM_VAR_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_VAR_le_dados ( unit )

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

  end subroutine NUM_VAR_le_inicio

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

  subroutine NUM_VAR_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_VAR_le_reinicio 

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

  subroutine NUM_VAR_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_VAR_le_rhg 

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

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

    real*8, 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_VAR_Mach

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

  subroutine NUM_VAR_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_VAR_massa_especifica_media

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

  subroutine NUM_VAR_mi ( n, T, a_mi, mi )

    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) :: a_mi ! coeficientes da viscosidade absoluta 

    real*8,dimension(:),intent(inout) :: mi ! viscosidade absoluta (Pa.s)

    integer :: j ! nível de temperatura: 1>=1000K, 2<1000K

    select case ( tipo_mi )

       case ( 1 )
          mi = mi_cte

       case ( 2 )
          do i = 1, n
             if ( T(i) >= 1.0d+3) then
                j = 1
             else
                j = 2
             end if
             mi(i) = dexp((a_mi(1,j)*dlog(T(i)) + a_mi(2,j)/T(i)   &
                   + a_mi(3,j)/(T(i)**2) + a_mi(4,j))) * 1.0d-7
          end do

    end select
   
  end subroutine NUM_VAR_mi

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

  subroutine NUM_VAR_mi_le_coeficientes ( unit, a_mi )

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

    real*8,dimension(:,:),intent(out) :: a_mi ! coeficientes da viscosidade absoluta 

    integer :: j ! nível de temperatura: 1>=1000K, 2<1000K
    integer :: l ! número do coeficiente

    open(unit, file = 'propri2.dat' )

    read(unit,*)
    read(unit,*)   
	    
    do j=1,2
       read(unit,*) ( a_mi(l,j), l=1,4 )
    end do

    close(unit)

  end subroutine NUM_VAR_mi_le_coeficientes

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

  subroutine NUM_VAR_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_VAR_pressao_entrada

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

  subroutine NUM_VAR_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_VAR_p_tot

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

  subroutine NUM_VAR_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_VAR_simplec

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

  subroutine NUM_VAR_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_VAR_T_tot

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

  subroutine NUM_VAR_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_VAR_T_wall

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

  subroutine NUM_VAR_uds_coef_pressao ( n,  dt, Rg_cte, 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_cte ! 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_cte * T(i-1) )

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

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

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

  subroutine NUM_VAR_uds_coef_temperatura ( n,  dt,   T_in,  Rg_cte, epsilon_cte, sigma,  &
                                            ro, roa,  rom,   u,      um,          T,      &
                                            Ta, p,    pa,    se,     sp,          rp,     &
                                            xe, xp,   sn,    f,      h,           T_wall, &
                                            cp, gama, C_rec, 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) :: T_in        ! temperatura na entrada da tubeira (K)
    real*8, intent(in) :: Rg_cte      ! constante do gás (J/kg.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)  :: gama   ! razão de calores específicos (adimensional)
    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 :: 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)

    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) + ae(i) ) &
             + h(i) * C_rec(i) * sn(i)                              &
             + epsilon_cte * sigma * sn(i) * ( T(i) ** 3 )    

       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 )

    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_VAR_uds_coef_temperatura

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

  subroutine NUM_VAR_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_VAR_uds_coef_velocidade

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

  subroutine NUM_VAR_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_VAR_velocidade_media

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

end module NUMERICO_VARIAVEL