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 ( 1 Fev 00)

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

! versão atual    1.0 ( 1 Fev 00)

! última alteração =   11 Abr 00

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

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

  ! NUM_VAR_cds_coef_temperatura
  ! NUM_VAR_cds_coef_velocidade
  ! NUM_VAR_compres
  ! NUM_VAR_correcao_com_plinha
  ! NUM_VAR_cp
  ! NUM_VAR_cp_le_coeficientes
  ! NUM_VAR_escreve_campos_auxiliares
  ! NUM_VAR_eta
  ! NUM_VAR_executa
  ! NUM_VAR_f
  ! NUM_VAR_grava_mach1d
  ! NUM_VAR_grava_reinicio
  ! NUM_VAR_grava_rhg
  ! NUM_VAR_h
  ! NUM_VAR_k
  ! NUM_VAR_k_le_coeficientes
  ! NUM_VAR_le_dados
  ! NUM_VAR_le_inicio
  ! NUM_VAR_le_mach1d
  ! NUM_VAR_le_reinicio
  ! NUM_VAR_le_rhg
  ! NUM_VAR_massa_especifica_media
  ! NUM_VAR_mi
  ! NUM_VAR_mi_le_coeficientes
  ! NUM_VAR_qc
  ! NUM_VAR_saida_rex
  ! NUM_VAR_simplec
  ! 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_mach1d
  ! 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_mach1d
  ! NUM_VAR_grava_reinicio
  ! NUM_VAR_grava_rhg
  ! NUM_VAR_saida_rex

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

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

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

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

implicit none

real*8,private :: T_in        ! temperatura na entrada do canal (K)
real*8,private :: fm_t        ! fluxo de massa em todos os canais (kg/s)
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 base do canal (K)
real*8,private :: eta_cte     ! eficiência das aletas (adimensional)
real*8,private :: compres_cte ! coeficiente de expansão térmica (1/K)
real*8,private :: rugoso      ! rugosidade absoluta da parede do canal (m)

! massa específica (ro) do refrigerante (kg/m3)
real*8,private :: ro_a ! coeficiente "a" de ro = a + b*T + c*(T**2), onde T=temperatura (K)
real*8,private :: ro_b ! coeficiente "b" de ro = a + b*T + c*(T**2), onde T=temperatura (K)
real*8,private :: ro_c ! coeficiente "c" de ro = a + b*T + c*(T**2), onde T=temperatura (K)

! condutividade térmica (k_wall) da parede do canal (W/m.K)
real*8,private :: kw_a ! coeficiente "a" de k_wall = a + b*T + c*(T**2), onde T=temperatura (K)
real*8,private :: kw_b ! coeficiente "b" de k_wall = a + b*T + c*(T**2), onde T=temperatura (K)
real*8,private :: kw_c ! coeficiente "c" de k_wall = a + b*T + c*(T**2), onde T=temperatura (K)

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 na base do canal
integer,private :: tipo_eta     ! tipo de eficiência das aletas (adimensional)
integer,private :: tipo_compres ! tipo de coeficiente de expansão térmica (1/K)

integer,private :: RHG    ! 0=executa CANAL; 1=executa CANAL acoplado ao MACH1D
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 :: ro_in ! massa específica na entrada do canal (kg/m3)
real*8,private :: u_in  ! velocidade na entrada do canal (m/s)

real*8,private :: T_ex  ! temperatura na saída do canal (K)

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

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

character*50,private :: nome_arq_do_Mach1D   ! nome do arquivo de dados do programa MACH1D
                                             ! para cálculo de T_wall

character*50,private :: nome_arq_para_Mach1D ! nome do arquivo de saída para o programa MACH1D

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

contains

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

  subroutine NUM_VAR_cds_coef_temperatura ( n, beta, cp, rom, T, um, se, bpp, 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) :: bpp ! 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-1) - T(i) )

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

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

    end do

    ! volume n (fictício)

    bc(n) = 0.0d0

    ! atualiza o termo independente

    bpp = bpp + bc

  end subroutine NUM_VAR_cds_coef_temperatura

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

  subroutine NUM_VAR_cds_coef_velocidade ( n, beta, rom, u, um, se, bpp, 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) :: bpp ! 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-1) - u(i) )

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

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

    end do

    ! volume n (fictício)

    bc(n) = 0.0d0

    ! atualiza o termo independente

    bpp = bpp + bc

  end subroutine NUM_VAR_cds_coef_velocidade

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

  subroutine NUM_VAR_compres ( ro, T, compres )

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

    real*8,dimension(:),intent(out) :: compres ! coeficiente de expansão térmica (1/K)

    select case ( tipo_compres )

       case ( 1 )
          compres = compres_cte

       case ( 2 )
          compres = - ( ro_b + 2.0d0 * ro_c * T ) / ro

    end select
    
  end subroutine NUM_VAR_compres

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

  subroutine NUM_VAR_correcao_com_plinha ( n, beta, dpp, de, xp, pl, u, p, um )

    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)    :: dpp ! coeficiente SIMPLEC no centro do volume
    real*8,dimension(:),intent(in)    :: de  ! coeficiente SIMPLEC na face leste do volume
    real*8,dimension(:),intent(in)    :: xp  ! coordenada espacial no centro do volume de controle

    real*8,dimension(:),intent(inout) :: 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) :: um  ! velocidade média (m/s)

    real*8 :: fator ! auxiliar

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

    ! pressão
    p = p + pl

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

    ! velocidade na entrada do canal
    um(n-1) = u_in

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

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

    u(n) =  u_in + beta * ( u_in - u(n-1) )

  end subroutine NUM_VAR_correcao_com_plinha

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

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

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

    real*8 :: Rg_cte = 461.5254d+0

    select case ( tipo_cp )

       case ( 1 )
          cp = cp_cte

       case ( 2 )
          cp = (a_cp(1) + a_cp(2)*T + a_cp(3)*(T**2)  &
             +  a_cp(4)*(T**3) + a_cp(5)*(T**4)) * Rg_cte

    end select
    
  end subroutine NUM_VAR_cp

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

  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
    integer :: sub ! substância: 7 = H2O líquido

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

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

    close(unit)

  end subroutine NUM_VAR_cp_le_coeficientes

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

  subroutine NUM_VAR_escreve_campos_auxiliares ( unit,   n,  w,   xp,      qlt,    f,   h,  &
                                                 T_wall, cp, eta, compres, k_wall, Red, Pr, &
                                                 Nu,     mi, k,   pl,      um,     rom, xe, &
                                                 T_tub   )

    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,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) :: 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 do lado do refrigerante (K)
    real*8,dimension(:),intent(in) :: cp      ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),intent(in) :: eta     ! eficiência das aletas (adimensional)
    real*8,dimension(:),intent(in) :: compres ! coeficiente de expansão térmica (1/K)
    real*8,dimension(:),intent(in) :: k_wall  ! condutividade térmica da parede do canal (W/m.K)
    real*8,dimension(:),intent(in) :: Red     ! número de Reynolds (adim)
    real*8,dimension(:),intent(in) :: Pr      ! número de Prandtl (adim)
    real*8,dimension(:),intent(in) :: Nu      ! número de Nusselt (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) :: 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)
    real*8,dimension(:),intent(in) :: T_tub   ! temperatura da parede da tubeira (K)

    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"_wc (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

    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 do lado do gás ***       ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'T_wh (K)             ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            T_tub                    ) ! variável a escrever

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** temperatura da parede do lado do refrigerante *', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'T_wc (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
                            '*** 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
                            '*** eficiência das aletas ***                      ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'eta (adimensional)   ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            eta                      ) ! variável a escrever

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

    call ARQUIVO_escreve_campo_1D ( unit,            & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** condutividade térmica da parede do canal ***   ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'k_wall (W/m.K)       ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            k_wall                   ) ! 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
                            '*** número de Nusselt ***                          ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'Nu (adimensional)    ', & ! 21 colunas
                            2,                       & ! número do ponto inicial
                            n-1,                     & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            Nu                       ) ! 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
                            '*** 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_eta ( n, xe, le, tp, re, ee, bp, h, k_wall, eta )

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

    real*8,dimension(:),intent(in)  :: xe     ! coordenada espacial na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: le     ! comprimento médio do escoamento no canal (m)
    real*8,dimension(:),intent(in)  :: tp     ! espessura das aletas do canal (m)
    real*8,dimension(:),intent(in)  :: re     ! raio da tubeira (m)
    real*8,dimension(:),intent(in)  :: ee     ! espessura entre tubeira e canal (m)
    real*8,dimension(:),intent(in)  :: bp     ! altura do canal (m)
    real*8,dimension(:),intent(in)  :: h      ! coeficiente de convecção (W/m2.K)
    real*8,dimension(:),intent(in)  :: k_wall ! condutividade térmica da parede do canal (W/m.K)

    real*8,dimension(:),intent(out) :: eta    ! eficiência das aletas (adimensional)

    real*8 :: u, base

    select case ( tipo_eta )

       case ( 1 )
          eta = eta_cte

       case ( 2 )
          do i = 2, n-1
             base = dsqrt ( ((xe(i)-xe(i-1))**2)              &
                  + ((re(i-1)+ee(i-1)-re(i)-ee(i))**2) )
             u = dsqrt ( h(i) * 2.0d0 * ( le(i) - le(i-1) )   & 
               / ( k_wall(i) * tp(i) * base ) )
             eta(i) = dtanh ( u * bp(i) ) / ( u * bp(i) )
          end do
          eta(1) = 0.0d0
          eta(n) = 0.0d0

    end select
    
  end subroutine NUM_VAR_eta

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

  subroutine NUM_VAR_executa ( n,         m,       w_it, beta,  nitm,  tipo_ci, w_aux, &
                               lista_aux, nome_ci, dt,   xe,    xp,    se,      sp,    &
                               dp,        le,      sbp,  sap,   tp,    re,      ep,    &
                               ee,        bp,      fm,   r_qwc, q_sai, u,       p,     &
                               T,         ro,      itgr, itgmaxr  )

    integer,intent(in) :: n         ! número de volumes de controle com dois fictícios
    integer,intent(in) :: m         ! número de canais em torno da tubeira
    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,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) :: 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) :: dp  ! diâmetro hidráulico da seção transversal do canal (m)
    real*8,dimension(:),intent(in) :: le  ! comprimento médio do escoamento no canal (m)
    real*8,dimension(:),intent(in) :: sbp ! área da base do canal (m2)
    real*8,dimension(:),intent(in) :: sap ! área molhada das aletas do canal (m2)
    real*8,dimension(:),intent(in) :: tp  ! espessura das aletas do canal (m)
    real*8,dimension(:),intent(in) :: re  ! raio da tubeira (m)
    real*8,dimension(:),intent(in) :: ep  ! espessura entre tubeira e canal (m)
    real*8,dimension(:),intent(in) :: ee  ! espessura entre tubeira e canal (m)
    real*8,dimension(:),intent(in) :: bp  ! altura do canal (m)

    real*8,intent(out) :: fm     ! fluxo de massa num único canal (kg/s)
    real*8,intent(out) :: r_qwc  ! razão q_wall / q_sai
    real*8,intent(out) :: q_sai  ! taxa de transferência de calor advectivo na saída do canal (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,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 de arquivos auxiliares
    integer :: unit_ou = 83 ! identificador do arquivo dos campos auxiliares

    integer :: i_itg

    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 :: bpp ! 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 :: dpp ! coeficiente SIMPLEC no centro do volume
    real*8,dimension(:),allocatable :: de  ! coeficiente SIMPLEC na face leste do volume

    real*8,dimension(:),allocatable :: qlt ! fluxo de calor total para a parede (W/m2)
    
    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 do lado do refrigerante (K)
    real*8,dimension(:),allocatable :: eta     ! eficiência das aletas (adimensional)
    real*8,dimension(:),allocatable :: compres ! coeficiente de expansão térmica (1/K)
    real*8,dimension(:),allocatable :: k_wall  ! condutividade térmica da parede do lado do refrigerante (W/m.K)
    real*8,dimension(:),allocatable :: Red     ! número de Reynolds (adim)
    real*8,dimension(:),allocatable :: Pr      ! número de Prandtl (adim)
    real*8,dimension(:),allocatable :: Nu      ! número de Nusselt (adim)

    real*8,dimension(:),allocatable :: Taw   ! temperatura de parede da tubeira adiabática (K)
    real*8,dimension(:),allocatable :: Tg    ! temperatura do gás na tubeira (K)
    real*8,dimension(:),allocatable :: T_tub ! temperatura da parede da tubeira (K)
    real*8,dimension(:),allocatable :: qh    ! taxa de transf. calor total para a parede da tubeira (W)
    real*8,dimension(:),allocatable :: s_tub ! área da parede da tubeira (m2)
    real*8,dimension(:),allocatable :: hg    ! coeficiente de convecção do gás na tubeira (W/m2.K)
    real*8,dimension(:),allocatable :: qc    ! taxa de transf. calor total para o refrigerante (W)
    real*8,dimension(:),allocatable :: Rg    ! resistência térmica do gás na tubeira (K/W)

    real*8 :: epsig     ! constante do MACH1D
    real*8 :: Cd        ! coeficiente de descarga (adimensional)
    real*8 :: F_star    ! empuxo dinâmico (adimensional)
    real*8 :: T_tub_max ! temperatura máxima da parede da tubeira (K)
    real*8 :: dq_max    ! diferença total entre qh e qc (%)
    real*8 :: auxiliar  !

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

    ! alocação de memória  
    allocate ( ua(n), pa(n), Ta(n), pl(n), roa(n), dpp(n) )
    allocate ( um(n-1), uma(n-1), de(n-1), rom(n-1) )
    allocate ( aw(n), ap(n), ae(n), bpp(n) )
    allocate ( bc(n), bf(n) )
    allocate ( cp(n), f(n), h(n), mi(n), k(n), T_wall(n), eta(n), compres(n), &
                                               k_wall(n), Red(n), Pr(n),      Nu(n) )
    allocate ( T_tub(n), Taw(n), qh(n), s_tub(n), hg(n), Tg(n), qc(n), Rg(n) )

    fm = fm_t / m

    ro_in = ro_a + ro_b * T_in + ro_c * (T_in**2)

    u_in = fm / ( ro_in * se(n-1) )

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

    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) // '.canal.itn'
       call ARQUIVO_cria_arquivo ( unit_it, comando )
       write ( unit_it, 1 ) 
       1 format (t1,'itt, p_sai-p_entra, T_saída')
    end if

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

    T_wall = T_wall_cte

    ! lê arquivo de dados do programa MACH1D
    if ( tipo_T_wall == 2 ) &
       call NUM_VAR_le_mach1d ( n,   m,  unit_au, nome_arq_do_Mach1D, epsig, Cd, F_star, &
                                Taw, Tg, T_tub,   hg,                 s_tub, qh          )

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

    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

       ! cálculo das propriedades
       call NUM_VAR_mi ( T, a_mi, mi )
       call NUM_VAR_k  ( T, a_k,  k  )
       call NUM_VAR_cp ( n, T, a_cp, cp )
       Red = ro * dabs(u) * dp / mi
       Pr  = cp * mi / k
       call NUM_VAR_f ( dp, Red, f )
       call NUM_VAR_h ( dp, Red, f, Pr, k, Nu, h )
       call NUM_VAR_compres ( ro, T, compres )
       k_wall = kw_a + kw_b * T_wall + kw_c * (T_wall**2)
       call NUM_VAR_eta ( n, xe, le, tp, re, ee, bp, h, k_wall, eta )

       ! solução da velocidade nodal
       call NUM_VAR_uds_coef_velocidade ( n,  beta, dt, ro,  roa, u, ua, um, p,  se,  &
                                          sp, dp,   le, rom, xp,  f, aw, ap, ae, bpp, bf )
       if ( beta > 0 ) &
          call NUM_VAR_cds_coef_velocidade ( n, beta, rom, u, um, se, bpp, bc )
       call SOLVERS_tdma ( n, aw, ap, ae, bpp, u )

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

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

       ! solução da temperatura
       call NUM_VAR_qc ( n,  m,  epsig, Taw, Tg,  T_tub,  T_wall, T, hg, h, s_tub, &
                         qh, ep, sbp,   sap, eta, k_wall, qc,     Rg )
       call NUM_VAR_uds_coef_temperatura ( n,  beta,    dt, ro, roa, rom, u,  um, T,  &
                                           Ta, p,       pa, se, sp,  dp,  le, xp, f,  &
                                           cp, compres, qc, aw, ap,  ae,  bpp )
       if ( beta > 0 ) &
          call NUM_VAR_cds_coef_temperatura ( n, beta, cp, rom, T, um, se, bpp, bc )
       call SOLVERS_tdma ( n, aw, ap, ae, bpp, T )

       ! calcula temperaturas nas paredes
       if ( tipo_T_wall == 2 ) &
          call NUM_VAR_T_wall ( n, m, Taw, qc, Rg, s_tub, ep, k_wall, T_tub, T_wall )

       ! cálculo da massa específica
       ro = ro_a + ro_b * T + ro_c * (T**2)
       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, ro, rom, roa, um,  de, se, &
                                       sp, le, aw,  ap,  ae, bpp, fm  )
       call SOLVERS_tdma ( n, aw, ap, ae, bpp, pl )

       ! correções das variáves com plinha
       call NUM_VAR_correcao_com_plinha ( n, beta, dpp, de, xp, pl, u, p, um )

       ! escreve variáveis monitoradas
       if ( w_it /= 0 .and. itg == itgmax ) then
          if ( itt==1 .or. itt==nitm .or. mod(itt,w_it)==0 ) then
             T_ex  = 0.5d0 * ( T(1) + T(2) )
             write ( unit_it, 2 ) itt, 0.5d0*(p(1)+p(2)-p(n-1)-p(n)), T_ex
             2 format ( i10, 2 ( 1pe23.15 ) )
          end if
       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,      T_wall   )

    ! entrada do canal
    u(n)  = u_in
    T(n)  = T_in
    ro(n) = ro_in
    p(n)  = 0.5d0*(p(n-1)+p(n))

    ! saída do canal
    u(1)  = um(1)
    T(1)  = 0.5d0 * ( T(1) + T(2) )
    ro(1) = rom(1)
    p(1)  = 0.5d0*(p(1)+p(2))

    deallocate ( ua, pa, ta, roa, uma, dpp, de )
    deallocate ( aw, ap, ae, bpp )
    deallocate ( bc, bf )
    deallocate ( Taw, s_tub, hg, Tg, Rg)

    ! cria o arquivo dos campos auxiliares
    if ( w_aux /= 0 .and. itg == itgmax ) then
       call ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.canal.outros.txt'
       call ARQUIVO_cria_arquivo ( unit_ou, comando )
       write ( unit_ou, 3 ) 
       3 format (/,t5,'*** CAMPOS AUXILIARES ***')
    end if

    ! calcula os fluxos de calor para a parede
    allocate ( qlt(n) )
    qlt = h * ( T_wall - T )

    call NUM_VAR_cp ( n, T, a_cp, cp )
    q_sai = - fm * cp(1) * T_ex

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

    r_qwc = r_qwc / q_sai

    if ( RHG == 1 ) call NUM_VAR_grava_mach1d ( unit_au, n, T_tub )

    if ( w_aux /= 0 .and. itg == itgmax ) &
       call NUM_VAR_escreve_campos_auxiliares ( unit_ou, n,  w_aux, xp,      qlt,    f,   h,  &
                                                T_wall,  cp, eta,   compres, k_wall, Red, Pr, &
                                                Nu,      mi, k,     pl,      um,     rom, xe, &
                                                T_tub    )

    qc(n)  = 0.0d0
    dq_max = 0.0d0
    T_tub_max = 0.0d0
    do i = 2, n-1
       auxiliar = dabs ( dabs ( qh(i) ) / m - dabs ( qc(i) ) )
       dq_max = dq_max + auxiliar
       qc(n) = qc(n) + dabs ( qc(i) )
       if ( T_tub(i) > T_tub_max ) T_tub_max = T_tub(i)
    end do

    if ( RHG == 1 ) then

       call ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.RHG.ITN'
       open ( unit_au, file = trim(comando) )
       do i_itg = 1, itg-1
          read ( unit_au, * )
       end do
       write ( unit_au, 4 ) itg, 1.0d+2 * dq_max / qc(n)
       4 format (i6,' = iteração', 10x, 1pe25.15,' = diferença total entre qh e qc (%)')
       close ( unit_au )

       call NUM_VAR_grava_rhg ( unit_au )

    end if

    call NUM_VAR_saida_rex ( unit_au, n, Cd, F_star, T_tub_max, p(1)-p(n), T_ex, xe )

    deallocate ( qlt, pl, um, rom )
    deallocate ( cp,  f,  h,  mi, k, T_wall, Red, Pr, Nu, eta, compres, k_wall )
    deallocate ( T_tub, qh, qc )

    itgr    = itg
    itgmaxr = itgmax

  end subroutine NUM_VAR_executa

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

  subroutine NUM_VAR_f ( dp, Red, f )

    real*8,dimension(:),intent(in)  :: dp  ! diâmetro hidráulico da seção transversal do canal (m)
    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/(3.7d0*dp) + 5.74d0/(Red**0.9d0) ) ) ** 2 ) )

    end select
    
  end subroutine NUM_VAR_f

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

  subroutine NUM_VAR_grava_mach1d ( unit, n, T_tub )

    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(in) :: T_tub ! temperatura da parede da tubeira (K)

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

    write ( unit ) n, T_tub

    close ( unit )

  end subroutine NUM_VAR_grava_mach1d

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

  subroutine NUM_VAR_grava_reinicio ( n,   unit, nome_ci, u,  p,  T,  ro,  um, &
                                      rom, pl,   T_wall   )

    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)
    real*8,dimension(:),intent(in) :: T_wall ! temperatura da parede do lado do refrigerante (K)

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

    write ( unit ) n, u, p, T, ro, um, rom, pl, T_wall

    close ( unit )

  end subroutine NUM_VAR_grava_reinicio 

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

  subroutine NUM_VAR_grava_rhg ( unit )

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

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

    write ( unit, * ) itg+1
    write ( unit, * ) itgmax

    close ( unit )

  end subroutine NUM_VAR_grava_rhg 

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

  subroutine NUM_VAR_h ( dp, Red, f, Pr, k, Nu, h )

    real*8,dimension(:),intent(in)  :: dp  ! diâmetro hidráulico da seção transversal do canal (m)
    real*8,dimension(:),intent(in)  :: Red ! número de Reynolds (adim)
    real*8,dimension(:),intent(in)  :: f   ! fator de atrito de Darcy (adim)
    real*8,dimension(:),intent(in)  :: Pr  ! número de Prandtl (adim)
    real*8,dimension(:),intent(in)  :: k   ! condutividade térmica (W/m.K)

    real*8,dimension(:),intent(out) :: Nu  ! número de Nusselt (adim)
    real*8,dimension(:),intent(out) :: h   ! coeficiente de convecção (W/m2.K)

    select case ( tipo_h )

       case ( 1 )
          h = h_cte

       case ( 2 )
          Nu = ( f / 8.0d0 ) * ( Red - 1.0d+3 ) * Pr  &
             / ( 1.0d0 + 12.7d0 * dsqrt(f/8.0d0) * ( (Pr**(2.0/3)) - 1.0d0 ) )
          h = k * dabs(Nu) / dp

    end select
    
  end subroutine NUM_VAR_h

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

  subroutine NUM_VAR_k ( T, a_k, k )

    real*8,dimension(:),intent(in)  :: T   ! temperatura (K)
    real*8,dimension(:),intent(in)  :: a_k ! coeficientes da condutividade térmica 

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

    select case ( tipo_k )

       case ( 1 )
          k = k_cte

       case ( 2 )
          k = a_k(1) + a_k(2)*T + a_k(3)*(T**2)

    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 :: l ! número do coeficiente

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

    do l = 1,9
       read(unit,*)
    end do   
	    
    do l=1,3
       read(unit,*) a_k(l)
    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, * ) T_in

    read ( unit, * ) ro_a
    read ( unit, * ) ro_b
    read ( unit, * ) ro_c

    read ( unit, * ) fm_t

    read ( unit, * ) tipo_cp
    read ( unit, * ) cp_cte

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

    read ( unit, * ) tipo_eta
    read ( unit, * ) eta_cte

    read ( unit, * ) tipo_compres
    read ( unit, * ) compres_cte

    read ( unit, * ) kw_a
    read ( unit, * ) kw_b
    read ( unit, * ) kw_c

    read ( unit, * ) RHG
    read ( unit, * ) nome_arq_para_Mach1D

    close ( unit )

  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_mach1d ( n,   m,     unit,  nome_arq_do_Mach1D, epsig, Cd, F_star, &
                                 Taw, Tg,    T_tub, hg,                 s_tub, qh          )

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

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

    real*8,dimension(:),intent(out) :: Taw   ! temperatura de parede da tubeira adiabática (K)
    real*8,dimension(:),intent(out) :: Tg    ! temperatura do gás na tubeira (K)
    real*8,dimension(:),intent(out) :: T_tub ! temperatura da parede da tubeira (K)
    real*8,dimension(:),intent(out) :: hg    ! coeficiente de convecção do gás na tubeira (W/m2.K)
    real*8,dimension(:),intent(out) :: s_tub ! área da parede da tubeira (m2)
    real*8,dimension(:),intent(out) :: qh    ! taxa de transf. calor total para a parede da tubeira (W)

    integer :: n_le

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

    read ( unit ) n_le, epsig, Cd, F_star, Taw, Tg, T_tub, hg, s_tub, qh

    if ( n_le /= n ) stop

    close ( unit )

  end subroutine NUM_VAR_le_mach1d 

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

  subroutine NUM_VAR_le_reinicio ( n,   unit, nome_ci, u,  p,  T,  ro, um,  &
                                   rom, pl,   T_wall   )

    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)
    real*8,dimension(:),intent(out) :: T_wall ! temperatura da parede do lado do refrigerante (K)

    integer :: n_le

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

    read ( unit ) n_le, u, p, T, ro, um, rom, pl, T_wall

    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_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-2

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

    end do

    rom(n-1) = ro_in

  end subroutine NUM_VAR_massa_especifica_media

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

  subroutine NUM_VAR_mi ( T, a_mi, mi )

    real*8,dimension(:),intent(in)  :: T    ! temperatura (K)
    real*8,dimension(:),intent(in)  :: a_mi ! coeficientes da viscosidade absoluta 

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

    select case ( tipo_mi )

       case ( 1 )
          mi = mi_cte

       case ( 2 )
          mi = 1.0d-3 * dexp ( a_mi(1) + a_mi(2)/T   &
             + a_mi(3)*T + a_mi(4)*(T**2) )

    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 :: l ! número do coeficiente

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

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

    close(unit)

  end subroutine NUM_VAR_mi_le_coeficientes

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

  subroutine NUM_VAR_qc ( n,  m,  epsig, Taw, Tg,  T_tub,  T_wall, T, hg, h, s_tub, &
                          qh, ep, sbp,   sap, eta, k_wall, qc,     Rg )

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: m    ! número de canais em torno da tubeira

    real*8,intent(in) :: epsig
 
    real*8,dimension(:),intent(in) :: Taw    ! temperatura de parede da tubeira adiabática (K)
    real*8,dimension(:),intent(in) :: Tg     ! temperatura do gás na tubeira (K)
    real*8,dimension(:),intent(in) :: T_tub  ! temperatura da parede da tubeira (K)
    real*8,dimension(:),intent(in) :: T_wall ! temperatura da parede do lado do refrigerante (K)
    real*8,dimension(:),intent(in) :: T      ! temperatura (K)
    real*8,dimension(:),intent(in) :: hg     ! coeficiente de convecção do gás na tubeira (W/m2.K)
    real*8,dimension(:),intent(in) :: h      ! coeficiente de convecção (W/m2.K)
    real*8,dimension(:),intent(in) :: s_tub  ! área da parede da tubeira (m2)
    real*8,dimension(:),intent(in) :: qh     ! taxa de transf. calor total para a parede da tubeira (W)
    real*8,dimension(:),intent(in) :: ep     ! espessura entre tubeira e canal (m)
    real*8,dimension(:),intent(in) :: sbp    ! área da base do canal (m2)
    real*8,dimension(:),intent(in) :: sap    ! área molhada das aletas do canal (m2)
    real*8,dimension(:),intent(in) :: eta    ! eficiência das aletas (adimensional)

    real*8,dimension(:),intent(out) :: k_wall ! condutividade térmica média da parede (W/m.K)
    real*8,dimension(:),intent(out) :: qc     ! taxa de transf. calor total para o refrigerante (W)
    real*8,dimension(:),intent(out) :: Rg     ! resistência térmica do gás na tubeira (K/W)

    real*8,dimension(:),allocatable :: Rtot   ! resistência térmica total (K/W)

    allocate ( Rtot(n) )

    select case ( tipo_T_wall )

       case ( 1 )
          Rtot = 1.0d0 / ( h * ( sbp + eta * sap ) )
          qc = ( T_wall - T ) / Rtot

       case ( 2 )
          k_wall = kw_a + kw_b * ( T_wall + T_tub ) / 2.0d0          &
                 + ( kw_c / 3.0d0 ) * ( (T_tub**3) - (T_wall**3) )   &
                 / ( T_tub - T_wall )
          Rg = ( ( Taw - T_tub ) / s_tub ) / ( hg * ( Taw - T_tub )   &
             + epsig * ( (Tg**4) - (T_tub**4) ) )
          Rtot = Rg + ep / ( k_wall * s_tub ) + 1.0d0   &
               / ( h * m * ( sbp + eta * sap ) )
          qc = ( Taw - T ) / Rtot

    end select

    qc = qc / m

    deallocate ( Rtot )

  end subroutine NUM_VAR_qc 

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

  subroutine NUM_VAR_saida_rex ( unit, n, Cd, F_star, T_tub_max, dp, T_ex, xe )

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

    real*8,intent(in) :: Cd        ! coeficiente de descarga (adimensional)
    real*8,intent(in) :: F_star    ! empuxo dinâmico (adimensional)
    real*8,intent(in) :: T_tub_max ! temperatura máxima da parede da tubeira (K)
    real*8,intent(in) :: dp        ! queda de pressão do escoamento (Pa)
    real*8,intent(in) :: T_ex      ! temperatura do escoamento na saída dos canais (K)

    real*8,dimension(:),intent(in) :: xe  ! coordenada espacial na face leste do volume de controle

    real*8 :: dx ! tamanho médio dos volumes de controle (m)

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

    call ARQUIVO_passa_caso ( caso )

    open ( unit, file = trim(caso) //'.canal.num' )

    dx = ( xe(n-1) - xe(1) ) / ( n - 2 )

    ! número total de variáveis = 2 ou 5
    if ( RHG == 0 ) then
       write(unit,1) dx, 2
    else
       write(unit,1) dx, 5
    end if
    1 format(1pd25.15,' = dx',//, &
                   i5,' = número total de variáveis',/)

    if ( RHG == 1 ) then

       ! parâmetro 1 = coeficiente de descarga
       write(unit,11) Cd
       11 format(3x,"'coeficiente de descarga da tubeira (adimensional)'",/, &
                  1pd25.15,' = solução numérica',/)

       ! parâmetro 2 = empuxo dinâmico (adimensional)
       write(unit,12) F_star
       12 format(3x,"'empuxo dinâmico na saída da tubeira (adimensional)'",/, &
                  1pd25.15,' = solução numérica',/)

       ! parâmetro 3 = temperatura máxima da parede da tubeira (K)
       write(unit,13) T_tub_max
       13 format(3x,"'temperatura máxima da parede da tubeira (K)'",/, &
                  1pd25.15,' = solução numérica',/)

    end if

    ! parâmetro 4 = pressão do escoamento: saída - entrada (Pa)
    write(unit,14) dp
    14 format(3x,"'pressão do escoamento nos canais: saída - entrada (Pa)'",/, &
               1pd25.15,' = solução numérica',/)

    ! parâmetro 5 = temperatura do escoamento na saída dos canais (K)
    write(unit,15) T_ex
    15 format(3x,"'temperatura do escoamento na saída dos canais (K)'",/, &
               1pd25.15,' = solução numérica',/)

    close(unit)

  end subroutine NUM_VAR_saida_rex

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

  subroutine NUM_VAR_simplec ( n, sp, aw, ap, ae, dpp, 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) :: dpp ! coeficiente SIMPLEC no centro do volume
    real*8,dimension(:),intent(out) :: de  ! coeficiente SIMPLEC na face leste do volume

    dpp(1) = 0.0d0
    dpp(n) = 0.0d0

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

    de(1)   = dpp(2)
    de(n-1) = 0.0d0

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

  end subroutine NUM_VAR_simplec

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

  subroutine NUM_VAR_T_wall ( n, m, Taw, qc, Rg, s_tub, ep, k_wall, T_tub, T_wall )

    integer,intent(in) :: n ! número de volumes de controle com dois fictícios
    integer,intent(in) :: m ! número de canais em torno da tubeira

    real*8,dimension(:),intent(in)  :: Taw    ! temperatura de parede da tubeira adiabática (K)
    real*8,dimension(:),intent(in)  :: qc     ! taxa de transf. calor total para o refrigerante (W)
    real*8,dimension(:),intent(in)  :: Rg     ! resistência térmica do gás na tubeira (K/W)
    real*8,dimension(:),intent(in)  :: s_tub  ! área da parede da tubeira (m2)
    real*8,dimension(:),intent(in)  :: ep     ! espessura entre tubeira e canal (m)
    real*8,dimension(:),intent(in)  :: k_wall ! condutividade térmica média da parede (W/m.K)

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

    T_tub = Taw - qc * m * Rg

    T_wall = T_tub - qc * m * ep  / ( k_wall * s_tub )

  end subroutine NUM_VAR_T_wall

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

  subroutine NUM_VAR_uds_coef_pressao ( n,  dt, ro, rom, roa, um,  de, se, &
                                        sp, le, aw,  ap,  ae, bpp, fm  )

    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)  :: 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)  :: 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)  :: le  ! comprimento médio do escoamento no canal (m)

    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) :: bpp ! termo independente

    real*8,intent(in) :: fm ! fluxo de massa num único canal (kg/s)

    real*8 :: dl ! tamanho do volume de controle elementar

    ! volume 1 (fictício): saída

    aw(1)  = 0.0d0
    ap(1)  = 1.0d0
    ae(1)  = 1.0d0
    bpp(1) = 0.0d0

    ! volumes internos

    do i = 2, n-2

       dl = le(i) - le(i-1)

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

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

       ap(i) = - ( aw(i) + ae(i) )

       bpp(i) = - ( ( ro(i) - roa(i) ) * sp(i) * dl / dt     &
              + rom(i) * um(i) * se(i) - rom(i-1) * um(i-1) * se(i-1) )

    end do

    ! volume interno adjacente ao contorno de entrada

    i = n-1

    dl = le(i) - le(i-1)

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

    ae(i) = 0.0d0

    ap(i) = - ( aw(i) + ae(i) )

    bpp(i) = - ( ( ro(i) - roa(i) ) * sp(i) * dl / dt     &
           + fm - rom(i-1) * um(i-1) * se(i-1) )

    ! volume n (fictício): entrada

    aw(n)  =  0.0d0
    ap(n)  =  1.0d0
    ae(n)  =  0.0d0
    bpp(n) =  0.0d0

  end subroutine NUM_VAR_uds_coef_pressao

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

  subroutine NUM_VAR_uds_coef_temperatura ( n,  beta,    dt, ro, roa, rom, u,  um, T,  &
                                            Ta, p,       pa, se, sp,  dp,  le, xp, f,  &
                                            cp, compres, qc, aw, ap,  ae,  bpp )

    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) :: 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)  :: 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)  :: dp      ! diâmetro hidráulico da seção transversal do canal (m)
    real*8,dimension(:),intent(in)  :: le      ! comprimento médio do escoamento no canal (m)
    real*8,dimension(:),intent(in)  :: xp      ! coordenada espacial no centro do volume de controle
    real*8,dimension(:),intent(in)  :: f       ! fator de atrito de Darcy (adim)
    real*8,dimension(:),intent(in)  :: cp      ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),intent(in)  :: compres ! coeficiente de expansão térmica (1/K)
    real*8,dimension(:),intent(in)  :: qc      ! taxa de transf. calor total para o refrigerante (W)

    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) :: bpp     ! termo independente

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

    pi = dacos(-1.0d0)

    ! volume 1 (fictício): saída

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

    ! volumes internos

    do i = 2, n-1

       dl = le(i) - le(i-1)

       aw(i) = 0.0d0

       ae(i) = cp(i) * rom(i) * um(i) * se(i)

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

       bpp(i) = cp(i) * roa(i) * sp(i) * dl * Ta(i) / dt                  &
              + compres(i) * T(i) * sp(i) * ( dl * ( p(i) - pa(i) ) / dt  &
              + u(i) * ( p(i+1) - p(i-1) ) * 0.5d0 )                   &
              - pi * f(i) * ro(i) * u(i) * (u(i)**2) * dp(i) * dl / 8     &
              + qc(i)

    end do

    ! volume n (fictício): entrada

    aw(n)  =  0.0d0
    ap(n)  =  1.0d0
    ae(n)  =  0.0d0
    bpp(n) =  T_in + beta * ( T_in - T(n-1) )

  end subroutine NUM_VAR_uds_coef_temperatura

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

  subroutine NUM_VAR_uds_coef_velocidade ( n,  beta, dt, ro,  roa, u, ua, um, p,  se,  &
                                           sp, dp,   le, rom, xp,  f, aw, ap, ae, bpp, bf )

    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) :: 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)  :: dp  ! diâmetro hidráulico da seção transversal do canal (m)
    real*8,dimension(:),intent(in)  :: le  ! comprimento médio do escoamento no canal (m)
    real*8,dimension(:),intent(in)  :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in)  :: xp  ! coordenada espacial no centro do volume de controle
    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) :: bpp ! termo independente
    real*8,dimension(:),intent(out) :: bf  ! termo de fricção para um

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

    pi = dacos(-1.0d0)

    ! volume 1 (fictício): saída

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

    ! volumes internos

    do i = 2, n-1

       dl = le(i) - le(i-1)

       aw(i) = 0.0d0

       ae(i) = rom(i) * um(i) * se(i)

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

       bf(i) = - pi * f(i) * ro(i) * u(i) * dabs(u(i)) * dp(i) * dl / 8   

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

    end do

    ! volume n (fictício): entrada

    aw(n)  =  0.0d0
    ap(n)  =  1.0d0
    ae(n)  =  0.0d0
    bpp(n) =  u_in + beta * ( u_in - u(n-1) )
    bf(n)  =  0.0d0

  end subroutine NUM_VAR_uds_coef_velocidade

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

  subroutine NUM_VAR_velocidade_media ( n,  dt, roa, u,  uma, p,  se, sp, le, &
                                        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)    :: le  ! comprimento médio do escoamento no canal (m)
    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) * (le(i) - le(i-1))

       massa_e = roa(i+1) * sp(i+1) * (le(i+1) - le(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) = u_in

  end subroutine NUM_VAR_velocidade_media

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

end module NUMERICO_VARIAVEL