module CANAL_NUMERICO_CONSTANTE

! sigla = NUM_CTE_

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

! versão original 1.0 (28 Jan 00)

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

! versão atual    1.0 (28 Jan 00)

! última alteração =   11 Abr 00

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

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

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

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

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

  ! NUM_CTE_le_dados
  ! NUM_CTE_le_inicio

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

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

  ! NUM_CTE_executa

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

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

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

use dados_rhg

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

implicit none

real*8,private :: T_in    ! temperatura na entrada do canal (K)
real*8,private :: ro_in   ! massa específica na entrada do canal (kg/m3)
real*8,private :: fm_t    ! fluxo de massa em todos os canais (kg/s)
real*8,private :: cp      ! calor específico à pressão constante (J/kg.K)
real*8,private :: f       ! fator de atrito de Darcy (adim)
real*8,private :: h       ! coeficiente de convecção (W/m2.K)
real*8,private :: T_wall  ! temperatura da parede na base do canal (K)
real*8,private :: eta     ! eficiência das aletas (adimensional)
real*8,private :: compres ! coeficiente de expansão térmica (1/K)

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)

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

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

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

contains

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

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

    integer,intent(in) :: n    ! número de volumes de controle com dois fictícios
    integer,intent(in) :: beta ! 0=UDS; 1=CDS
    real*8, intent(in) :: cp   ! calor específico à pressão constante (J/kg.K)

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

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

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

    real*8 :: oeste, leste ! auxiliares

    ! volume 1 (fictício)

    bc(1) = 0.0d0

    ! volumes internos

    do i = 2, n-1

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

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

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

    end do

    ! volume n (fictício)

    bc(n) = 0.0d0

    ! atualiza o termo independente

    bp = bp + bc

  end subroutine CANAL_NUM_CTE_cds_coef_temperatura

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

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

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

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

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

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

    real*8 :: oeste, leste ! auxiliares

    ! volume 1 (fictício)

    bc(1) = 0.0d0

    ! volumes internos

    do i = 2, n-1

       oeste = rom(i-1) * um(i-1) * se(i-1) * ( u(i-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

    bp = bp + bc

  end subroutine CANAL_NUM_CTE_cds_coef_velocidade

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

  subroutine CANAL_NUM_CTE_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 CANAL_NUM_CTE_correcao_com_plinha

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

  subroutine CANAL_NUM_CTE_escreve_campos_auxiliares ( unit, n, w, xp, xe, qlt, pl, um, rom )

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

    real*8,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) :: 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 CANAL_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 CANAL_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 CANAL_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 CANAL_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 CANAL_NUM_CTE_escreve_campos_auxiliares

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

  subroutine CANAL_NUM_CTE_executa ( n,   m,   w_it, beta,  nitm,  w_aux, lista_aux, &
                                     dt,  xe,  xp,   se,    sp,    dp,    le,        &
                                     sbp, sap, fm,   r_qwc, q_sai, u,     p,         &
                                     T,   ro   )

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

    real*8,intent(in) :: dt ! intervalo de tempo (s)

    real*8,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,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 :: unit_it = 81 ! identificador do arquivo das iterações
    integer :: unit_au = 82 ! identificador do arquivo dos campos auxiliares

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

    character*(50) :: caso    ! nome do arquivo do caso a simular
    character*220  :: comando ! linha de comando

! *** campos auxiliares ***

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

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

    real*8,dimension(:),allocatable :: 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)

    ! 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), bp(n) )
    allocate ( bc(n), bf(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

    ! cria o arquivo das iterações
    if ( w_it /= 0 ) then
       call CANAL_ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.canal.itn'
       call CANAL_ARQUIVO_cria_arquivo ( unit_it, comando )
       write ( unit_it, 1 ) 
       1 format (t1,'itt, p_sai-p_entra, T_saída')
    end if

    call CANAL_ARQUIVO_cronometro ( tcpu_can1 )

    ! 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

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

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

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

       ! solução da temperatura
       call CANAL_NUM_CTE_uds_coef_temperatura ( n,   beta,    dt,  f,   h,   cp, T_wall, &
                                                 eta, compres, ro,  roa, rom, u,  um,     &
                                                 T,   Ta,      p,   pa,  se,  sp, dp,     &
                                                 le,  sbp,     sap, xp,  aw,  ap, ae,     bp )
       if ( beta > 0 ) &
          call CANAL_NUM_CTE_cds_coef_temperatura ( n, beta, cp, rom, T, um, se, bp, bc )
       call TDMA ( n, aw, ap, ae, bp, T )

       ! cálculo da massa específica
       ro = ro_a + ro_b * T + ro_c * (T**2)
       call CANAL_NUM_CTE_massa_especifica_media ( n, beta, ro, rom )

       ! solução da correção da pressão
       call CANAL_NUM_CTE_uds_coef_pressao ( n,  dt, ro, rom, roa, um, de, se, &
                                             sp, le, aw,  ap,  ae, bp, fm  )
       call TDMA ( n, aw, ap, ae, bp, pl )

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

       ! escreve variáveis monitoradas
       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 do

    call CANAL_ARQUIVO_cronometro ( tcpu_can2 )

    if ( w_it /= 0 ) then
       call CANAL_ARQUIVO_escreve_tempo_cpu ( tcpu_can1, tcpu_can2, unit_it )
       close ( unit_it )
    end if

    ! 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, bp )
    deallocate ( bc, bf )

    ! cria o arquivo dos campos auxiliares
    if ( w_aux /= 0 ) then
       call CANAL_ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.canal.outros.txt'
       call CANAL_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 ( qlt(n) )
    qlt = h * ( T_wall - T )

    if ( w_aux /= 0 ) &
       call CANAL_NUM_CTE_escreve_campos_auxiliares ( unit_au, n, w_aux, xp, xe, qlt, pl, um, rom )

    q_sai = - fm * cp * 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 * sap(i) )
    end do

    r_qwc = r_qwc / q_sai
  
    deallocate ( qlt, pl, um, rom )

  end subroutine CANAL_NUM_CTE_executa

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

  subroutine CANAL_NUM_CTE_le_dados ( unit )

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

    open ( unit, file = trim(nome_num_cte) )

    read ( unit, * ) T_in 
    read ( unit, * ) ro_a
    read ( unit, * ) ro_b
    read ( unit, * ) ro_c
    read ( unit, * ) fm_t
    read ( unit, * ) cp
    read ( unit, * ) f
    read ( unit, * ) h
    read ( unit, * ) T_wall
    read ( unit, * ) eta
    read ( unit, * ) compres

    close ( unit )

  end subroutine CANAL_NUM_CTE_le_dados 

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

  subroutine CANAL_NUM_CTE_le_inicio ( unit, nome )

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

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

    open ( unit, file = trim(nome) )

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

    close ( unit )

    if ( (itg_rhg == 1).and.(lista == 1) ) call CANAL_ARQUIVO_notepad_abre_dado ( nome_num_cte )

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

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

  end subroutine CANAL_NUM_CTE_le_inicio

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

  subroutine CANAL_NUM_CTE_massa_especifica_media ( n, beta, ro, rom )

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

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

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

    do i = 1, n-2

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

    end do

    rom(n-1) = ro_in

  end subroutine CANAL_NUM_CTE_massa_especifica_media

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

  subroutine CANAL_NUM_CTE_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 CANAL_NUM_CTE_simplec

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

  subroutine CANAL_NUM_CTE_uds_coef_pressao ( n,  dt, ro, rom, roa, um, de, se, &
                                              sp, le, aw,  ap,  ae, bp, 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) :: bp  ! 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
    bp(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) )

       bp(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) )

    bp(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
    bp(n) =  0.0d0

  end subroutine CANAL_NUM_CTE_uds_coef_pressao

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

  subroutine CANAL_NUM_CTE_uds_coef_temperatura ( n,   beta,    dt,  f,   h,   cp, T_wall, &
                                                  eta, compres, ro,  roa, rom, u,  um,     &
                                                  T,   Ta,      p,   pa,  se,  sp, dp,     &
                                                  le,  sbp,     sap, xp,  aw,  ap, ae,     bp )

    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,intent(in) :: f       ! fator de atrito de Darcy (adim)
    real*8,intent(in) :: h       ! coeficiente de convecção (W/m2.K)
    real*8,intent(in) :: cp      ! calor específico à pressão constante (J/kg.K)
    real*8,intent(in) :: T_wall  ! temperatura da parede da tubeira (K)
    real*8,intent(in) :: eta     ! eficiência das aletas (adimensional)
    real*8,intent(in) :: compres ! coeficiente de expansão térmica (1/K)

    real*8,dimension(:),intent(in)  :: ro  ! massa específica (kg/m3)
    real*8,dimension(:),intent(in)  :: roa ! massa específica anterior (kg/m3)
    real*8,dimension(:),intent(in)  :: rom ! massa específica média (kg/m3)
    real*8,dimension(:),intent(in)  :: u   ! velocidade (m/s)
    real*8,dimension(:),intent(in)  :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in)  :: T   ! temperatura (K)
    real*8,dimension(:),intent(in)  :: Ta  ! temperatura anterior (K)
    real*8,dimension(:),intent(in)  :: p   ! pressão (Pa)
    real*8,dimension(:),intent(in)  :: pa  ! pressão anterior (Pa)
    real*8,dimension(:),intent(in)  :: se  ! área na face leste do volume de controle
    real*8,dimension(:),intent(in)  :: sp  ! área no centro do volume de controle
    real*8,dimension(:),intent(in)  :: 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)  :: xp  ! coordenada espacial no centro do volume de controle

    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 :: 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
    bp(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 * rom(i) * um(i) * se(i)

       ap(i) = cp * roa(i) * sp(i) * dl / dt - ( aw(i) + ae(i) ) &
             + h * ( sbp(i) + eta * sap(i) )

       bp(i) = cp * roa(i) * sp(i) * dl * Ta(i) / dt                  &
             + compres * T(i) * sp(i) * ( dl * ( p(i) - pa(i) ) / dt  &
             + u(i) * ( p(i+1) - p(i-1) ) * 0.5d0 )                   &
             - pi * f * ro(i) * u(i) * (u(i)**2) * dp(i) * dl / 8     &
             + h * T_wall * ( sbp(i) + eta * sap(i) )

    end do

    ! volume n (fictício): entrada

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

  end subroutine CANAL_NUM_CTE_uds_coef_temperatura

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

  subroutine CANAL_NUM_CTE_uds_coef_velocidade ( n,  beta, dt, f,  ro,  roa, u,  ua, um, p,  &
                                                 se, sp,   dp, le, rom, xp,  aw, ap, ae, bp, 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, intent(in) :: f  ! fator de atrito de Darcy (adim)

    real*8,dimension(:),intent(in)  :: ro  ! massa específica (kg/m3)
    real*8,dimension(:),intent(in)  :: roa ! massa específica anterior (kg/m3)
    real*8,dimension(:),intent(in)  :: u   ! velocidade (m/s)
    real*8,dimension(:),intent(in)  :: ua  ! velocidade anterior (m/s)
    real*8,dimension(:),intent(in)  :: um  ! velocidade média (m/s)
    real*8,dimension(:),intent(in)  :: p   ! pressão (Pa)
    real*8,dimension(:),intent(in)  :: se  ! área na face leste do volume de controle (m2)
    real*8,dimension(:),intent(in)  :: sp  ! área no centro do volume de controle (m2)
    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(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 :: 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
    bp(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 * ro(i) * u(i) * dabs(u(i)) * dp(i) * dl / 8   

       bp(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
    bp(n) =  u_in + beta * ( u_in - u(n-1) )
    bf(n) =  0.0d0

  end subroutine CANAL_NUM_CTE_uds_coef_velocidade

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

  subroutine CANAL_NUM_CTE_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 CANAL_NUM_CTE_velocidade_media

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

end module CANAL_NUMERICO_CONSTANTE