module ref_radiativa

use dados_rhg
use paredes
use portlib

implicit none

! Variáveis gerais
real*8  ::  eps_par_amb ! emissividade térmica parede-ambiente (adim.)
real*8  ::  sigma       ! constante de Stefan-Boltzmann (W/m2.K4)
real*8  ::  epsig_rad   ! produto de eps_par_amb por sigma (W/m2.K4)
real*8  ::  Tamb        ! temperatura do ambiente (K)
real*8  ::  kw_a_rad    ! coeficiente "a" de k_wall = a + b*T + c*(T**2), onde T=temperatura (K)
real*8  ::  kw_b_rad    ! coeficiente "b" de k_wall = a + b*T + c*(T**2), onde T=temperatura (K)
real*8  ::  kw_c_rad    ! coeficiente "c" de k_wall = a + b*T + c*(T**2), onde T=temperatura (K)
real*8  ::  Twall_presc ! temperatura da parede (prescrita) (K)

integer :: tipo_Twall_rad  ! tipo de parede: 1 - temperatura prescrita; 2 - temperatura variável

! propriedades geométricas no centro dos volumes de controle
real*8,dimension(:),allocatable :: xp_rad  ! coordenada espacial (m)
real*8,dimension(:),allocatable :: ep_rad  ! espessura da parede da tubeira (m)
real*8,dimension(:),allocatable :: sbp_rad ! área externa da tubeira (m2)
real*8,dimension(:),allocatable :: sp_rad  ! área no centro do volume de controle

! propriedades geométricas na face leste dos volumes de controle
real*8,dimension(:),allocatable :: xe_rad  ! coordenada espacial (m)
real*8,dimension(:),allocatable :: re_rad  ! raio da tubeira (m)
real*8,dimension(:),allocatable :: ee_rad  ! espessura da parede da tubeira (m)
real*8,dimension(:),allocatable :: se_rad  ! área na face leste do volume de controle

!  outros
character*100 :: titulo_rad         ! título do caso a simular
character*50  :: nome_caso_rad      ! nome do arquivo do caso a simular
character*50  :: nome_parede_rad    ! nome do arquivo de dados do programa PAREDE

integer :: w_T_rad       ! freqüência de escrita de T; 0=não escreve

integer :: w_lista_rad          ! 1 = cria listagem do caso; 0=não
integer :: lista_caso_rad       ! lista conteúdo do arquivo do caso
integer :: lista_rex_rad        ! lista conteúdo do arquivo de saída para o programa REX

integer,private :: itg    ! número da iteração global corrente entre MACH2D e RAD

character*50, private :: nome_caso  ! nome do caso

character*220,private :: comando    ! linha de comando

contains

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

subroutine RAD_NUMERICO_ARQUIVO_cria_arquivo ( unit, nome_arquivo )

    integer,intent(in)       :: unit          ! identificador do arquivo
    character*(*),intent(in) :: nome_arquivo  ! nome do arquivo

    comando = trim(nome_arquivo)
   
    open ( unit, file = comando )

end subroutine RAD_NUMERICO_ARQUIVO_cria_arquivo

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

subroutine RAD_NUMERICO_ARQUIVO_cronometro ( tcpu )

    real*8,intent(inout) :: tcpu ! tempo de processamento

    call cpu_time (tcpu) 

end subroutine RAD_NUMERICO_ARQUIVO_cronometro

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

subroutine RAD_NUMERICO_ARQUIVO_escreve_campo_1D ( unit,   w,      head_0, head_1, head_2, &
                                        head_3, n_ini,  n,      x,      T       )

    integer,intent(in)             :: unit   ! identificador do arquivo
    integer,intent(in)             :: w      ! freqüência de escrita; 0=não escreve
    integer,intent(in)             :: n_ini  ! ponto ou volume inicial de escrita
    character*(*),intent(in)       :: head_0 ! título do campo
    character*(*),intent(in)       :: head_1 ! título da primeira coluna
    character*(*),intent(in)       :: head_2 ! título da segunda coluna
    character*(*),intent(in)       :: head_3 ! título da terceira coluna
    integer,intent(in)             :: n      ! número de dados a escrever
    real*8,intent(in),dimension(:) :: x      ! coordenada espacial
    real*8,intent(in),dimension(:) :: T      ! variável a escrever

    real*8  :: T_min, T_max ! extremos do campo
    real*8  :: x_min, x_max ! coordenadas dos extremos do campo
    integer :: i ! número do dado

!   head_0 = inicia na coluna 1  ; formato A50
!   head_1 = inicia na coluna 1  ; formato A6
!   head_2 = inicia na coluna 11 ; formato A21
!   head_3 = inicia na coluna 36 ; formato A21
!   inteiro: formato I6
!   real   : formato 1pe25.15

    if ( w == 0 ) return

    write(unit,1) head_0, head_1, head_2, head_3 
    1 format(//,t1, a50, &
             //,t1, a6,   t11, a21,   t36, a21, / )
    do i = n_ini, n
       if ( i==n_ini .or. i==n .or. mod(i,w)==0 ) &
          write(unit,2) i, x(i), T(i)
    end do

    2 format(i6,2(1pe25.15))

    T_min = T(n_ini)
    T_max = T(n_ini)

    do i = n_ini, n
       if ( T(i) <= T_min ) then
          T_min = T(i)
          x_min = x(i)
       end if
       if ( T(i) >= T_max ) then
          T_max = T(i)
          x_max = x(i)
       end if
    end do

    write(unit,3) T_min, x_min, T_max, x_max
    3 format(/,1pe25.15,' = mínimo',1pe25.15,' = X (m)',/, &
               1pe25.15,' = máximo',1pe25.15,' = X (m)') 

end subroutine RAD_NUMERICO_ARQUIVO_escreve_campo_1D

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

subroutine RAD_NUMERICO_ARQUIVO_escreve_tempo_cpu ( tcpu_i, tcpu_f, unit )

    integer,intent(in) :: unit ! arquivo
    real*8,intent(in)  :: tcpu_i ! zero do cronômetro 
    real*8,intent(in)  :: tcpu_f ! tempo de processamento

    write(unit,1) tcpu_f - tcpu_i
    1 format(/,1pe15.5,' = tempo de processamento (s)',/)

end subroutine RAD_NUMERICO_ARQUIVO_escreve_tempo_cpu

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

subroutine RAD_NUMERICO_ARQUIVO_notepad_abre_dado ( nome_arquivo )

    character*(*),intent(in) :: nome_arquivo ! nome do arquivo

    comando = 'notepad '// trim(nome_arquivo)

    dos = systemqq ( comando )

end subroutine RAD_NUMERICO_ARQUIVO_notepad_abre_dado

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

subroutine RAD_NUMERICO_ARQUIVO_passa_caso ( caso )

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

    caso = nome_caso_rad

end subroutine RAD_NUMERICO_ARQUIVO_passa_caso

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

subroutine RAD_NUMERICO_escreve_campos_auxiliares ( unit,  n,   w,    xp,  qlt,  T_wall,  k_wall, &
                                                    pl,    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) :: T_wall  ! temperatura da parede do lado externo (K)
    real*8,dimension(:),intent(in) :: k_wall  ! condutividade térmica da parede do canal (W/m.K)
    real*8,dimension(:),intent(in) :: pl      ! fluxo de calor total para a parede (W/m2)
    real*8,dimension(:),intent(in) :: T_tub   ! temperatura da parede da tubeira (K)

    call RAD_NUMERICO_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 RAD_NUMERICO_ARQUIVO_escreve_campo_1D ( unit,   & ! unidade de arquivo para escrever
                            w,                       & ! freqüência de escrita
                            '*** temperatura - parede interna ***                                ', & ! 50 colunas
                            ' ponto',                & !  6 colunas
                            'x (m)                ', & ! 21 colunas
                            'T (K)                ', & ! 21 colunas
                            1,                       & ! número do ponto inicial
                            nx,                      & ! número de dados a escrever
                            xp,                      & ! coordenada espacial
                            T_tub                    ) ! variável a escrever

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

end subroutine RAD_NUMERICO_escreve_campos_auxiliares

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

subroutine RAD_NUMERICO_executa ( n, nome_ci, xe, xp, se, sp, sbp, re, ep, ee  )

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

    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) :: sbp ! área da superfície externa da tubeira (m2)
    real*8,dimension(:),intent(in) :: re  ! raio da tubeira (m)
    real*8,dimension(:),intent(in) :: ep  ! espessura da parede da tubeira (m)
    real*8,dimension(:),intent(in) :: ee  ! espessura da parede da tubeira (m)

    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

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

    character*220  :: comando ! linha de comando

! *** campos auxiliares ***

    real*8,dimension(:),allocatable :: qlt ! fluxo de calor total para a parede (W/m2)
    
    real*8,dimension(:),allocatable :: T_wall  ! temperatura da parede do lado externo (K)
    real*8,dimension(:),allocatable :: k_wall  ! condutividade térmica da parede (W/m.K)

    real*8,dimension(:),allocatable :: Taw   ! temperatura de parede da tubeira (lado interno) (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 :: qr    ! 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  !

    integer :: i, iter
    integer :: unit = 40

    ! alocação de memória  
    allocate ( qlt(n), T_wall(n), k_wall(n) )
    allocate ( T_tub(n), Taw(n), qh(n), s_tub(n), hg(n), Tg(n), qr(n), Rg(n) )

    call RAD_NUMERICO_ARQUIVO_cronometro ( tcpu_rad1 )
    
    if (itg_rhg > 1) call RAD_NUMERICO_le_rhg ( unit_au )

    T_wall = Twall_presc

    ! lê arquivo de dados do programa MACH1D
    call RAD_NUMERICO_le_mach1d (n,    unit,  epsig,  Cd,  F_star, &
                                 Taw,  Tg,    T_tub,  hg,  s_tub,  qh )

    ! cálculo da taxa de transferência de calor
    call RAD_NUMERICO_qr ( n,  epsig, Taw, Tg,     T_tub, T_wall, T,  hg, h, s_tub, &
                           qh, ep,    sbp, k_wall, qr,    Rg,     re  )
 
    ! calcula temperaturas nas paredes
    call RAD_NUMERICO_T_wall ( n, Taw, qr, Rg, s_tub, ep, k_wall, T_tub, T_wall )

    deallocate ( Taw, s_tub, hg, Tg, Rg)

    ! calcula os fluxos de calor para a parede
    qlt = epsig_rad * ( T_wall**4 - Tamb**4 )


    call RAD_NUMERICO_grava_mach1d ( unit_au, n, T_tub )
 
    call RAD_NUMERICO_ARQUIVO_cronometro ( tcpu_rad2 )
    
    ! cria o arquivo dos campos auxiliares
    if ( itg_rhg == itg_rhg_max ) then
        call RAD_NUMERICO_ARQUIVO_passa_caso ( caso )
        comando = trim(caso) // '.rad.outros.txt'
        call RAD_NUMERICO_ARQUIVO_cria_arquivo ( unit_ou, comando )
        write ( unit_ou, 3 ) 
        3 format (/,t5,'*** CAMPOS AUXILIARES ***')
    end if

    if ( itg_rhg == itg_rhg_max ) then
        call RAD_NUMERICO_escreve_campos_auxiliares ( unit_ou,  n,   w_lista_rad,  xp,  qlt,  T_wall,  k_wall, &
                                                      pl,    xe,  T_tub )
    
        call RAD_NUMERICO_ARQUIVO_escreve_tempo_cpu ( tcpu_rad1, tcpu_rad2, unit_ou )

        close ( unit_ou )

        
        if ( lista_caso_rad == 1 ) then
           comando = 'notepad ' // trim(caso)//'.rad.outros.txt' 
           dos = system( comando )
        end if

    end if

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

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

    ! Tempos de processamento:
    if ( itg_rhg == itg_rhg_max) then
        tcpu_mach = tcpu_mach2 - tcpu_mach1
        tcpu_mach = tcpu_mach  * itg_rhg_max
        tcpu_rad  = tcpu_rad2  - tcpu_rad1
        tcpu_rad  = tcpu_rad   * itg_rhg_max
        tcpu_tot  = tcpu_mach  + tcpu_rad
    end if

    call RAD_NUMERICO_grava_rhg ( unit_au )

    call RAD_NUMERICO_saida_rex ( unit_au, n, Cd, F_star, T_tub_max, xe )

    deallocate ( qlt, T_wall, k_wall )
    deallocate ( T_tub, qh, qr )
    
end subroutine RAD_NUMERICO_executa

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

subroutine RAD_NUMERICO_geral

    ! dados do programa PAREDE
    if ( itg_rhg == 1 ) allocate ( xp_rad(nx),   ep_rad(nx),   sp_rad(nx),   sbp_rad(nx) )
    if ( itg_rhg == 1 ) allocate ( xe_rad(nx-1), re_rad(nx-1), ee_rad(nx-1), se_rad(nx-1) )

    xe_rad = xe_par
    xp_rad = xp_par
    se_rad = se_par
    sp_rad = sp_par
    sbp_rad = sbp_par
    re_rad = re_par
    ee_rad = ee_par
    ep_rad = ep_par

    if ( itg_rhg == 1 ) call RAD_NUMERICO_le_dados

    call RAD_NUMERICO_executa ( nx,     caso,   xe_rad, xp_rad,    se_rad,    sp_rad,    sbp_rad, &
                                re_rad, ep_rad, ee_rad  )


end subroutine RAD_NUMERICO_geral

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

  subroutine RAD_NUMERICO_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
    
    integer :: i

    real*8,dimension(:),intent(in) :: T_tub ! temperatura da parede da tubeira (K)

    open ( unit, file = 'temp_parede_interna.bin', form = 'unformatted' )

    write ( unit ) T_tub

    close ( unit )
    
  end subroutine RAD_NUMERICO_grava_mach1d

! ---------------------------------------------------------
  subroutine RAD_NUMERICO_grava_rhg ( unit )

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

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

    write ( unit, * ) itg_rhg+1
    write ( unit, * ) itg_rhg_max

    close ( unit )

  end subroutine RAD_NUMERICO_grava_rhg 

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

subroutine RAD_NUMERICO_le_dados

    integer :: dos
    
    dos = system ('notepad RHG2D_radiativo_geral.ent')

    open(5,file='RHG2D_radiativo_geral.ent')
    
    read(5,*) nome_caso_rad
    read(5,*) eps_par_amb
    read(5,*) sigma
    read(5,*) Tamb
    read(5,*) tipo_Twall_rad
    read(5,*) Twall_presc
    read(5,*) kw_a_rad
    read(5,*) kw_b_rad
    read(5,*) kw_c_rad
    read(5,*) w_lista_rad          
    read(5,*) lista_caso_rad       
    read(5,*) lista_rex_rad        
    
    close(5)

    epsig_rad = eps_par_amb * sigma

end subroutine RAD_NUMERICO_le_dados

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

subroutine RAD_NUMERICO_le_rhg ( unit )

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

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

    read ( unit, * ) !itg_rhg
    read ( unit, * ) itg_rhg_par

    close ( unit )

end subroutine RAD_NUMERICO_le_rhg 

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

subroutine RAD_NUMERICO_le_mach1d ( n,    unit,  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) :: unit ! unidade de leitura do arquivo de dados
 
    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(n),intent(out) :: Taw  ! temperatura de parede da tubeira adiabática (K)
    real*8,dimension(n),intent(out) :: Tg    ! temperatura do gás na tubeira (K)
    real*8,dimension(n),intent(out) :: T_tub ! temperatura da parede da tubeira (K)
    real*8,dimension(n),intent(out) :: hg    ! coeficiente de convecção do gás na tubeira (W/m2.K)
    real*8,dimension(n),intent(out) :: s_tub ! área da parede da tubeira (m2)
    real*8,dimension(n),intent(out) :: qh    ! taxa de transf. calor total para a parede da tubeira (W)

    integer :: n_le

    open ( unit,file = 'transf.mach.canal.bin',form='unformatted' )
    
    read ( unit ) epsig, Cd, F_star, Taw, Tg, T_tub, hg, s_tub, qh
    
    close ( unit )

  end subroutine RAD_NUMERICO_le_mach1d 

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

subroutine RAD_NUMERICO_qr ( n,  epsig, Taw, Tg,     T_tub, T_wall, T,  hg, h, s_tub, &
                             qh, ep,    sbp, k_wall, qr,    Rg,     re  )

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

    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 - lado interno (K)
    real*8,dimension(:),intent(inout) :: T_wall ! temperatura da parede do lado externo (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 (interna) 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 da parede (m)
    real*8,dimension(:),intent(in) :: sbp    ! área da parede (externa) da tubeira (m2)
    real*8,dimension(:),intent(in) :: re     ! raio da tubeira (m)

    real*8,dimension(:),intent(out) :: k_wall ! condutividade térmica média da parede (W/m.K)
    real*8,dimension(:),intent(out) :: qr     ! taxa de transf. calor total para o espaço(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)
    
    real    :: pi
    integer :: i
    
    pi = dacos(-1.0d0)

    allocate ( Rtot(n) )

    select case ( tipo_Twall_rad )

       case ( 1 )
          Rtot = 1.0d0 / ( epsig_rad * s_tub * (T_wall**4 - Tamb**4) )
          qr = ( T_wall - Tamb ) / Rtot

       case ( 2 )
       

          k_wall = kw_a_rad + kw_b_rad * ( T_wall + T_tub ) / 2.0d0          &
                 + ( kw_c_rad / 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 ) 
          
          qr = ( Taw - T_wall ) / Rtot
          T_wall = ( Tamb**4 + qr / ( epsig_rad * s_tub ) ) ** 0.25d0
          
          Rtot = Rg + ep / ( k_wall * s_tub ) + 1.0d0 / &
                 ( epsig_rad * s_tub * (T_wall**2 + Tamb**2) * (T_wall + Tamb) )  
          qr = ( Taw - Tamb ) / Rtot      
          
    end select
   
    deallocate ( Rtot )

end subroutine RAD_NUMERICO_qr 

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

subroutine RAD_NUMERICO_saida_rex ( unit, n, Cd, F_star, T_tub_max, 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,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 RAD_NUMERICO_ARQUIVO_passa_caso ( caso )

    open ( unit, file = trim(caso) //'.rad.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, 3
!    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

    ! tempo de processamento total
    write(unit,16) tcpu_mach, tcpu_can, tcpu_tot

    16 format(3x, "'tempo de processamento total (s)'", /, &
               1pe25.15,' = escoamento reativo', /, &
               1pe25.15,' = radiação na parede', /, &
               1pe25.15,' = tempo de processamento total',/)

    close(unit)

  end subroutine RAD_NUMERICO_saida_rex

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

subroutine RAD_NUMERICO_T_wall ( n, Taw, qr, Rg, s_tub, ep, k_wall, T_tub, T_wall )

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

    real*8,dimension(:),intent(in)  :: Taw    ! temperatura de parede da tubeira adiabática (K)
    real*8,dimension(:),intent(in)  :: qr     ! 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 externo (K)

    T_tub = Taw - qr * Rg

    T_tub(1) = T_tub(2)
    T_tub(n) = T_tub(n-1)

    T_wall = T_tub - qr * ep  / ( k_wall * s_tub )

    T_wall(1) = T_wall(2)
    T_wall(n) = T_wall(n-1)


end subroutine RAD_NUMERICO_T_wall

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

end module ref_radiativa