module usuario

! *** ESTE MÓDULO DEVE SER REVISADO A CADA NOVA SIMULAÇÃO
!     PARA EVENTUAIS ADAPTAÇÕES QUE SEJAM NECESSÁRIAS ***

! Objetivos: definir a estimativa inicial das incógnitas,
!            aplicar as condições de contorno,
!            definir funções para propriedades,
!            fazer pós-processamento.

! Versão atual: 15 Mar 09

! Relação das sub-rotinas em ordem de apresentação:
!    1) temperatura_da_parede
!    2) calor_especifico_e_gama
!    3) viscosidade_laminar_nos
!    4) condutividade_termica_nos
!    5) viscosidade_laminar_faces
!    6) condutividade_termica_faces
!    7) estimativa_inicial
!    8) ccu
!    9) ccv
!   10) ccT
!   11) ccp
!   12) Tin_pin_entrada
!   13) plin_entrada_p_ficticio
!   14) u_v_nos_ficticios_com_pl
!   15) Uce_Vcn_de_dn_faces_contornos
!   16) Uce_Vcn_faces_contornos_com_pl
!   17) contornos_pos_processamento
!   18) EXATO_MACH_fluxo_massa
!   19) EXATO_MACH_newton

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

use coeficientes
use dados_rhg

implicit none

real*8 :: pr = 1.01325d5   ! pressão atmosférica ao nível do mar

real*8 :: go = 9.80665d+0  ! aceleração gravitacional ao nível do mar

real*8,parameter,private :: Ru = 8.314510d+0 ! constante universal dos gases (J/mol.K)

integer,private :: ll ! número do coeficiente

integer,private :: mm ! faixa de temperatura; 1 para T >= 1000 K; 2 para T < 1000 K

real*8,dimension(:,:),allocatable :: cam  ! coeficientes para cálculo de cp (monoespécie, prop. variáveis)
real*8,dimension(:,:),allocatable :: cbm  ! coeficientes para viscosidade (monoespécie, prop. variáveis)
real*8,dimension(:,:),allocatable :: ccm  ! coeficientes para cond. térmica (monoespécie, prop. variáveis)

contains

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

  subroutine temperatura_da_parede

    ! definição de Twin na parede em cada nó nx cuja face norte
    ! coincide com o contorno norte

    ! só é considerado se ccTw = 1
    
    if (itg_rhg == 1) then

        Twin = Tw_cte

    else

        open(55,file='temp_parede_interna.bin',form='unformatted')
        read(55) Twin
        close(55)
    
    end if

  end subroutine temperatura_da_parede

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

  subroutine calor_especifico_e_gama

    ! definição do cp no nó de cada volume de controle
    
    select case (tipo_cp)
    case (1)
       cp = cp_cte
    case (2)
       if (it == 1) call leitura_coef_calor_especifico
       
       do np = 1, nxy

          if ( T(np) >= 1.0d+3 ) then
             mm = 1
          else
             mm = 2
          end if
      	  
          cp(np) = ( cam(1,mm) + cam(2,mm)*T(np) + cam(3,mm)*(T(np)**2)   &
                 + cam(4,mm)*(T(np)**3) + cam(5,mm)*(T(np)**4) ) * Ru

       end do     

       cp = cp*1.0d3 / 18.01528d+0

    end select

!    cp = 1.030d+3 ! Para ar a T = 500 K (Incropera)

    ! cálculo algébrico do gcp no nó de cada volume de controle

    gcp = cp / ( cp - Rg )

  end subroutine calor_especifico_e_gama

!-------------------------------------------------
  
  subroutine leitura_coef_calor_especifico
  
       if (itg_rhg == 1) allocate ( cam(7,2) )

       open(80,file='Gibbs_coeficientes_cp_h_g.dat')

       read(80,*) 
       read(80,*)   
       read(80,*) 
       read(80,*)
       do mm = 1, 2
          do ll = 1, 7
             read(80,*) cam(ll,mm)
          end do
          read(80,*)
       end do

       close(80)
  
  end subroutine leitura_coef_calor_especifico

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

  subroutine viscosidade_laminar_nos

    ! definição do vl no nó de cada volume de controle

    ! SEU VALOR NÃO PODE SER ZERADO POIS CAUSARÁ INDETERMINAÇÕES
    
    select case (tipo_vl)
    case (1)
        vlp = visc_cte
    case (2)
        if (it == 1) call leitura_coef_visc
      
       ! vlp (Pa.s)
       
       do np = 1, nxy

            if ( T(np) >= 1.0d+3 ) then
               mm = 1
            else
               mm = 2
            end if

	        vlp(np) = dexp( cbm(1,mm)*dlog(T(np)) + cbm(2,mm)/T(np)    &
                    + cbm(3,mm)/(T(np)**2) + cbm(4,mm) ) * 1.0d-7

        end do
   
    end select

!    vlp = 2.701d-5 ! Para ar a T = 500 K (Incropera)

  end subroutine viscosidade_laminar_nos

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

  subroutine condutividade_termica_nos

    ! definição do k no nó de cada volume de controle

    ! SEU VALOR NÃO PODE SER ZERADO POIS CAUSARÁ INDETERMINAÇÕES

    select case (tipo_k)
    case (1)
        kp = k_cte
    case(2)
        if (it == 1) call leitura_coef_k
        ! kp (W/m.K)
        
        do np = 1, nxy

            if ( T(np) >= 1.0d+3 ) then
               mm = 1
            else
               mm = 2
            end if
        	  
	        kp(np)  = dexp( ccm(1,mm)*dlog(T(np)) + ccm(2,mm)/T(np)    &
                    + ccm(3,mm)/(T(np)**2) + ccm(4,mm) ) * 1.0d-4

            end do
    
    end select

!    kp = 4.07d-2 ! Para ar a T = 500 K (Incropera)

  end subroutine condutividade_termica_nos

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

  subroutine leitura_coef_k
  
    if ( itg_rhg == 1 ) allocate ( ccm(4,2) )

    open(80,file='Gibbs_coeficientes_ki_mi.dat')

    read(80,*)
    read(80,*)
    do mm = 1, 2
       read(80,*) 
    end do
    do mm = 1, 2
       read(80,*) ( ccm(ll,mm), ll = 1, 4 )
    end do
    read(80,*)

    close(80)

  end subroutine leitura_coef_k

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

  subroutine leitura_coef_visc
  
    if ( itg_rhg == 1 ) allocate ( cbm(4,2) )

    open(80,file='Gibbs_coeficientes_ki_mi.dat')

    read(80,*)
    read(80,*)
    do mm = 1, 2
       read(80,*) ( cbm(ll,mm), ll = 1, 4 )
    end do

    close(80)

  end subroutine leitura_coef_visc

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

  subroutine viscosidade_laminar_faces

    integer npe, npn

    ! Média harmônica de Patankar (1980) para malha de faces centradas

    ! Nas faces que coincidem com os contornos do domínio de cálculo,
    ! se vlp(Tp), então vl-contornos deve ser função de T nos contornos

    ! vl na face leste

    do i = 1, nx-1

       do j = 2, ny-1

          np  = (j-1)*nx + i
          npe = np + 1

          vle(np) = ( 2 * vlp(np) * vlp(npe) ) / ( vlp(np) + vlp(npe) )

       end do

    end do

    ! vl na face norte

    do i = 2, nx-1

       do j = 1, ny-1

          np  = (j-1)*nx + i
          npn = np + nx

          vln(np) = ( 2 * vlp(np) * vlp(npn) ) / ( vlp(np) + vlp(npn) )

       end do

    end do

  end subroutine viscosidade_laminar_faces

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

  subroutine condutividade_termica_faces

    integer npe, npn

    ! Média harmônica de Patankar (1980) para malha de faces centradas

    ! Nas faces que coincidem com os contornos do domínio de cálculo,
    ! se kp(Tp), então k-contornos deve ser função de T nos contornos

    ! k na face leste

    do i = 1, nx-1

       do j = 2, ny-1

          np  = (j-1)*nx + i
          npe = np + 1

          ke(np) = ( 2 * kp(np) * kp(npe) ) / ( kp(np) + kp(npe) )

       end do

    end do

    ! k na face norte

    do i = 2, nx-1

       do j = 1, ny-1

          np  = (j-1)*nx + i
          npn = np + nx

          kn(np) = ( 2 * kp(np) * kp(npn) ) / ( kp(np) + kp(npn) )

       end do

    end do

  end subroutine condutividade_termica_faces

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

  subroutine estimativa_inicial

    ! solução analítica Q1D para tubeiras circulares

    integer :: ind, npe, npw, npww
    real*8  :: dyf, dyr, ys,  xw,   dxf,  dxr,  aux 
    real*8  :: Si,  exp, cte, ue1D, Me1D, Te1D, pe1D

    ! determina o fluxo de massa da solução analítica 1D
    call EXATO_MACH_fluxo_massa ( gama, Rg, po, T0, Sg, fm1D )

    ! zeragem dos campos
    u    = 0.0d0
    v    = 0.0d0
    ue   = 0.0d0
    ve   = 0.0d0
    un   = 0.0d0
    vn   = 0.0d0
    p    = 0.0d0
    T    = 0.0d0
    ro   = 0.0d0
    pl   = 0.0d0
    Uce  = 0.0d0
    Vcn  = 0.0d0
    de   = 0.0d0
    dn   = 0.0d0
    due  = 0.0d0
    dun  = 0.0d0
    dve  = 0.0d0
    dvn  = 0.0d0
    roe  = 0.0d0
    ron  = 0.0d0
    cup  = 0.0d0
    cvp  = 0.0d0
    Sup  = 0.0d0
    Svp  = 0.0d0
    plin = 0.0d0
    cp   = 0.0d0
    vlp  = 0.0d0
    kp   = 0.0d0
    M    = 0.0d0

    ! estimativa inicial de u, p e T nos nós
    ! e cálculo da solução analítica 1D nos nós reais

    exp = gama / (gama - 1)

    do i = 2, nx-1

       j = ny-1
       np = (j-1)*nx + i

       Si = pi * ( rn(np) ** 2 ) 

       ind = 0
       if ( i > ig ) ind = 1
       call EXATO_MACH_newton ( gama, Si/Sg, ind, M1D(i) )

       cte = 1 + (gama - 1) * (M1D(i)**2) / 2

       p1D(i) = po / (cte ** exp)

       T1D(i) = T0 / cte

       u1D(i) = M1D(i) * dsqrt ( gama * Rg * T1D(i) )

       do j = 1, ny

          np = (j-1)*nx + i

          u(np) = u1D(i)
          p(np) = p1D(i)
          T(np) = T1D(i)
          if ( j > 1 .and. j < ny-1 ) Vcn(np) = v(np)*xk(np) - u(np)*yk(np)
          if ( j > 1 .and. j < ny-1 ) un(np) = u1D(i)
  
       end do

    end do

    ! calcula u e v nos nós fictícios de acordo com as condições de contorno
    call u_v_nos_ficticios_com_pl

    ! inicializa T e p na entrada da tubeira com a estimativa inicial
    call Tin_pin_entrada

    ! calcula p e T nos nós fictícios de acordo com as condições de contorno na entrada
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       p(np)  = 2 * pin(j) - p(npe)
       T(np)  = 2 * Tin(j) - T(npe)
    end do

    ! calcula p e T nos nós fictícios de acordo com as condições de contorno na saída
    i = nx
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npw = np - 1
       npww  = np - 2
       xw    = x(npw)
       dxf   = 2 * (xw - xp(npw))
       dxr   = xp(npw) - xp(npww)
       p(np) = p(npw) + (dxf/dxr) * (p(npw) - p(npww))
       T(np) = T(npw) + (dxf/dxr) * (T(npw) - T(npww))
    end do

    ! estimativa inicial de u nas faces leste
    ! e cálculo da solução analítica 1D nos contornos

    do i = 1, nx-1

       j = ny-1
       np = (j-1)*nx + i

       Si = pi * ( raio(np) ** 2 ) 

       ind = 0
       if ( i > ig ) ind = 1
       call EXATO_MACH_newton ( gama, Si/Sg, ind, Me1D )

       cte = 1 + (gama - 1) * (Me1D**2) / 2

       pe1D = po / (cte ** exp)

       Te1D = T0 / cte

       ue1D = Me1D * dsqrt ( gama * Rg * Te1D )

       do j = 2, ny-1

          np = (j-1)*nx + i

          ue(np)  = ue1D
          Uce(np) = ue1D*ye(np) - v(np)*xe(np)
  
       end do

       ! contorno de entrada
       if ( i == 1 ) then
          u1D(i) = ue1D
          p1D(i) = pe1D
          T1D(i) = Te1D
          M1D(i) = Me1D
       end if

       ! contorno de saída
       if ( i == nx-1 ) then
          u1D(nx) = ue1D
          p1D(nx) = pe1D
          T1D(nx) = Te1D
          M1D(nx) = Me1D
          Fd1D  = fm1D * ue1D
          Fpv1D = pe1D * Si
       end if

    end do

    ! inicializa as contravariantes nos contornos
    call Uce_Vcn_de_dn_faces_contornos

    ! inicializa a massa específica
    call massa_especifica_nos
    call massa_especifica_faces

    ! inicializa os coeficientes e fontes
    a   = 0.0d0
    ap  = 0.0d0
    apu = 0.0d0
    b   = 0.0d0
    do np = 1, nxy
       a(np,5)  = 1.0d0
       ap(np,3) = 1.0d0
    end do

    if ( ci == 2 ) then
       open(99,file=trim(adjustl(caso))//'.bin',form='unformatted')
       if (modf == 1) read(99) ito,norma1,tcpuo,pin,u,v,p,pl,ue,vn,T,ro,ve,un,Uce,Vcn,roe,ron
       if (modf == 2) read(99) ito,norma1,tcpuo,pin,u,v,p,pl,ue,vn,T,ro,ve,un,Uce,Vcn,roe,ron,Yi
       itmax = ito + itmax
       ito = ito + 1
       close(99)
    end if

  end subroutine estimativa_inicial

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

  subroutine ccu

    integer :: npw, npww, npe, npee
    real*8  :: xw,  dxf,  dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       npe  = np + 1
       npee = np + 2
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,5) = 1
       a(np,6) =-1
       dxf = 2 * xp(npe)
       dxr = xp(npee) - xp(npe)
       b(np)   = (dxf/dxr) * (u(npe) - u(npee))
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,4) =-1
       a(np,5) = 1
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (u(npw) - u(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,5) = 1
       a(np,8) =-1
       b(np)   = 0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       if ( modvis == 1 ) then ! viscoso
          a(np,2) = 1
       else                    ! invíscido
          a(np,2) =-1
       end if
       a(np,5) = 1
       b(np)   = 0
    end do
  
  end subroutine ccu

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

  subroutine ccv

    integer :: npw, npww
    real*8  :: xw,  dxf,  dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,5) = 1
       a(np,6) = 1
       b(np)   = 0
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,4) =-1
       a(np,5) = 1
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (v(npw) - v(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,5) = 1
       a(np,8) = 1
       b(np)   = 0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       if ( modvis == 1 ) then ! viscoso
          a(np,2) = 1
       else                    ! invíscido
          a(np,2) =-1
       end if
       a(np,5) = 1
       b(np)   = 0
    end do
  
  end subroutine ccv
!-------------------------------------------------

  subroutine ccT

    integer :: npw, npww
    real*8  :: xw,  dxf,  dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,5) = 1
       a(np,6) = 1
       b(np) = 2 * Tin(j)
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,4) =-1
       a(np,5) = 1
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (T(npw) - T(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 1, nx
       np = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       a(np,5) = 1
       a(np,8) =-1
       b(np)   = 0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 1, nx
       np = (j-1)*nx + i
       do k = 1, 9
          a(np,k) = 0
       end do
       if ( (ccTw == 1).or.(ccTw == 2) ) then ! T prescrito
          a(np,2) = 1
          b(np)   = 2 * Twin(i)
       else                  ! adiabático
          a(np,2) =-1
          b(np)   = 0
       end if
       a(np,5) = 1
    end do
  
  end subroutine ccT

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

  subroutine ccp

    integer :: npn, npnn, nps, npss, npw,  npww
    real*8  :: xw, dxf, dxr

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       ap(np,1) = 0
       ap(np,2) = 0
       ap(np,3) = 1
       ap(np,4) = 1
       ap(np,5) = 0
       b(np)    = 2 * plin(j)
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       ap(np,1) = 0
       ap(np,2) =-1
       ap(np,3) = 1
       ap(np,4) = 0
       ap(np,5) = 0
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np) = (dxf/dxr) * (pl(npw) - pl(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np   = (j-1)*nx + i
       ap(np,1) = 0
       ap(np,2) = 0
       ap(np,3) = 1
       ap(np,4) = 0
       ap(np,5) =-1
       b(np)    = 0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       ap(np,1) =-1
       ap(np,2) = 0
       ap(np,3) = 1
       ap(np,4) = 0
       ap(np,5) = 0
       b(np)    = 0
    end do
  
  end subroutine ccp

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

  subroutine Tin_pin_entrada

    integer :: npe
    real*8  :: aux

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       uin(j) = ( u(np) + u(npe) ) / 2
       vin(j) = 0.0d0
       Tin(j) = T0 - ( gama - 1 ) * ( uin(j)**2 + vin(j)**2 ) / ( 2 * gama * Rg )
       Mw(j)  = dsqrt ( (uin(j)**2 + vin(j)**2) / (gama * Rg * Tin(j)) )
       aux    = 1 + ( gama - 1 ) * ( Mw(j) ** 2 ) / 2
       pin(j) = po / ( aux ** ( gama / ( gama - 1 ) ) )
    end do

  end subroutine Tin_pin_entrada

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

  subroutine plin_entrada_p_ficticio

    integer :: npe

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       p(np)   = 2 * pin(j) - p(npe)
       plin(j) = pin(j) - pina(j)
    end do

    end subroutine plin_entrada_p_ficticio

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

  subroutine u_v_nos_ficticios_com_pl

    integer :: npn, nps, npw, npww, npe, npee
    real*8  :: xw,  dxf, dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       npe  = np + 1
       npee = np + 2
       dxf = 2 * xp(npe)
       dxr = xp(npee) - xp(npe)
       u(np) = u(npe) + (dxf/dxr) * (u(npe) - u(npee))
       v(np) = - v(npe)
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       u(np) = u(npw) + (dxf/dxr) * (u(npw) - u(npww))
       v(np) = v(npw) + (dxf/dxr) * (v(npw) - v(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np  = (j-1)*nx + i
       npn = np + nx
       u(np) = u(npn)
       v(np) = -v(npn)
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np  = (j-1)*nx + i
       nps = np - nx
       if ( modvis == 1 ) then ! viscoso
          u(np) = -u(nps)
          v(np) = -v(nps)
       else                    ! invíscido
          u(np) = u(nps)
          v(np) = v(nps)
       end if
    end do
  
  end subroutine u_v_nos_ficticios_com_pl

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

  subroutine Uce_Vcn_de_dn_faces_contornos

    integer :: npw, npe
    real*8  :: ue,  uin

    ! velocidades Uce e de

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       uin = ( u(np) + u(npe) ) / 2
       Uce(np) = uin * ye(np)
       de(np)  = de(npe)
    end do

    ! contorno direito (saída)
    i = nx-1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npw = np - 1
       npe = np + 1
       ue  = ( u(np) + u(npe) ) / 2
       Uce(np) = ue * ye(np)
       de(np)  = de(npw)
    end do

    ! velocidades Vcn e dn

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np = (j-1)*nx + i
       Vcn(np) = 0.0d0
       dn(np)  = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny-1
    do i = 2, nx-1
       np  = (j-1)*nx + i
       Vcn(np) = 0.0d0
       dn(np)  = 0.0d0
    end do

  end subroutine Uce_Vcn_de_dn_faces_contornos

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

  subroutine Uce_Vcn_faces_contornos_com_pl

    ! velocidades Uce

    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       Uce(np) = Uce(np) + de(np) * (pl(np) - pl(np+1))
    end do

    i = nx-1
    do j = 2, ny-1
       np = (j-1)*nx + i
       Uce(np) = Uce(np) + de(np) * (pl(np) - pl(np+1))
    end do

  end subroutine Uce_Vcn_faces_contornos_com_pl

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

  subroutine contornos_pos_processamento

    integer :: npn, npw, nps, npnn, npe
    real*8  :: raz 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       nps = np - nx
       u(np) = uin(j)
       v(np) = vin(j)
       p(np) = pin(j)
       T(np) = Tin(j)
       xp(np) = x(np)
       yp(np) = ( y(nps) + y(np) ) / 2
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       u(np) = ( u(npw) + u(np) ) / 2
       v(np) = ( v(npw) + v(np) ) / 2
       p(np) = ( p(npw) + p(np) ) / 2
       T(np) = ( T(npw) + T(np) ) / 2
       xp(np) = x(npw)
       yp(np) = ( y(npw-nx) + y(npw) ) / 2
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np   = (j-1)*nx + i
       npw  = np - 1
       npn  = np + nx
       u(np) = u(npn)
       v(np) = 0.0d0
       p(np) = p(npn)
       T(np) = T(npn)
       xp(np) = ( x(npw) + x(np) ) / 2
       yp(np) = ( y(npw) + y(np) ) / 2
    end do

    ! canto esquerdo inferior
    j = 1
    i = 1
       np   = (j-1)*nx + i
       npn  = np + nx
       u(np) = u(npn)
       v(np) = 0.0d0
       p(np) = p(npn)
       T(np) = T(npn)
       xp(np) = x(np)
       yp(np) = y(np)

    ! canto direito inferior
    j = 1
    i = nx
       np   = (j-1)*nx + i
       npw  = np - 1
       npn  = np + nx
       u(np) = u(npn)
       v(np) = 0.0d0
       p(np) = p(npn)
       T(np) = T(npn)
       xp(np) = x(npw)
       yp(np) = y(npw)

    ! contorno superior (parede)
    j = ny
    do i = 1, nx
       np  = (j-1)*nx + i
       nps = np - nx
       u(np) = (u(np) + u(nps)) / 2
       v(np) = (v(np) + v(nps)) / 2
       p(np) = p(nps)
       if ( (ccTw == 1).or.(ccTw == 2) ) then ! T prescrito
          T(np) = Twin(i)
       else                  ! adiabático
          T(np) = T(nps)
       end if
       xp(np) = ( x(nps) + x(nps-1) ) / 2
       yp(np) = ( y(nps) + y(nps-1) ) / 2
    end do

    ! cantos superiores
    j = ny
    do i = 1, nx, nx-1
       np  = (j-1)*nx + i
       nps = np - nx
       npw = np - 1
       npe = np + 1
       p(np) = p(nps)
       if ( i == 1 ) then
          xp(np) = x(nps)
          yp(np) = y(nps)
          u(np) = u(npe)
          v(np) = v(npe)
       end if
       if ( i == nx ) then
          xp(np) = x(nps-1)
          yp(np) = y(nps-1)
          u(np) = u(npw)
          v(np) = v(npw)
       end if
    end do
  
  end subroutine contornos_pos_processamento

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

  subroutine EXATO_MACH_fluxo_massa ( gama, Rg, p_cam, T_cam, Sg, fm_exato )

    real*8,intent(in)  :: gama     ! razão de calores específicos (adimensional)
    real*8,intent(in)  :: Rg       ! constante do gás (J/kg.K)
    real*8,intent(in)  :: p_cam    ! pressão na câmara de combustão (Pa)
    real*8,intent(in)  :: T_cam    ! temperatura na câmara de combustão (K)
    real*8,intent(in)  :: Sg       ! área da garganta (m2)
    real*8,intent(out) :: fm_exato ! fluxo de massa exato (kg/s)

    real*8 :: a1, a2, a3 ! auxiliares

    a1 = gama / ( Rg * T_cam )

    a2 = 2.0d0 / ( gama + 1.0d0 )

    a3 = (gama + 1.0d0 ) / ( gama - 1.0d0 )

    fm_exato = p_cam * Sg * dsqrt ( a1 * (a2**a3) )

  end subroutine EXATO_MACH_fluxo_massa

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

  subroutine EXATO_MACH_newton ( gama, a, ind, Mp )

    real*8,intent(in)   :: a      ! razão de áreas (adimensional)
    real*8,intent(in)   :: gama   ! razão de calores específicos (adimensional)
    integer,intent(in)  :: ind    ! 0=escoamento subsônico; 1=supersônico
    real*8,intent(out)  :: Mp     ! solução do número de Mach analítico (adimensional)
    real*8 :: c, d, e, den, fat1, flinha ! variáveis auxiliares
    real*8 :: deltam ! estimativa do erro na solução de Mp

    integer :: nitm = 500        ! número máximo de iterações no método de Newton-Raphson
    real*8  :: tol  =   1.0d-15  ! tolerância no método de Newton-Raphson
    real*8  :: M_inicio_sub      = 1.0d-01 ! Mach inicial subsônico
    real*8  :: M_inicio_super    = 1.0d+01 ! Mach inicial supersônico
    real*8  :: M_inicio_negativo = 1.0d-02 ! Mach inicial se negativo
 
    if ( ind == 1 ) Mp = M_inicio_super
    if ( ind == 0 ) Mp = M_inicio_sub

    e = (gama - 1.0d0) / 2.0d0
    c = (gama + 1.0d0) / ( 2.0d0 - 2.0D0 * gama)
    d = ((gama + 1.0d0) / 2.0d0) ** c

    do j = 1, nitm
      fat1 = Mp * ( c * (1.0d0+e*Mp*Mp) ** (c - 1.0d0)    &
             * 2.0d0 * e * Mp)
      den = Mp * ( 1.0d0 + e * Mp * Mp) ** c
      flinha = - d * ((1.0d0+e*Mp*Mp) ** c + fat1) / (den * den)
      deltam = - (d/den-a) / flinha
      Mp = Mp + deltam
      if ( Mp < 0.0d0 ) Mp = M_inicio_negativo
      if ( dabs(deltam) < tol ) exit
    end do

  end subroutine EXATO_MACH_newton

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

  subroutine transf_parede
    
    real*8  :: pi, dx, dy, dl, perim, fmac, Pr, rec
    real*8  :: dr
    integer :: ij
    
    pi = dacos(-1.0d0)

    Taw = 0.0d0
    h = 0.0d0

    do i = 2, nx-1
        np = nxy - 2*nx + i
        perim = pi*( y(np)+y(np-1) )
        dx = x(np) - x(np-1)
        dy = y(np-1) - y(np)
        dl = dsqrt( dx**2 + dy**2 ) 
        s_tub(i) = perim*dl


        do j = 2, ny-1
            ij = (j-1)*nx + i
            dr = 0.5d0*(y(ij-1) + y(ij) - y(ij-nx-1) - y(ij-nx))
            Pr = vlp(i)*cp(ij)/kp(ij)
            rec = Pr ** (1.0d0/3.0d0)
            fmac = 0.5d0*( gcp(ij)-1.0d0 ) * M(ij)**2
            Taw(i) = Taw(i) + dr*yp(ij)*( T0 * ( ( 1.0d0 + rec*fmac ) / ( 1.0d0 + fmac ) ) )
            h(i) = h(i) + dr*yp(ij)*( 0.026d0*( (ro(ij)*u(ij))**0.8d0 ) * ( Pr**0.4d0 ) * kp(ij) / &
                   ( (vlp(ij)**0.8d0) * ((y(np)+y(np-1))**0.2d0) ) )
        end do
        
        Taw(i) = 8.0d0*Taw(i) / ( y(np) + y(np-1) )**2
        h(i) = 8.0d0*h(i) / ( y(np) + y(np-1) )**2
        

!        Pr = vlp(i)*cp(i)/kp(i)
!        rec = Pr ** (1.0d0/3.0d0)
!        fmac = 0.5d0*( gcp(i)-1.0d0 ) * M(i)**2
!        Taw(i) = T0 * ( ( 1.0d0 + rec*fmac ) / ( 1.0d0 + fmac ) )
!        h(i) = 0.026d0*( (ro(i)*u(i))**0.8d0 ) * ( Pr**0.4d0 ) * kp(i) / &
!               ( (vlp(i)**0.8d0) * ((y(np)+y(np-1))**0.2d0) )
        
        qlc(i) = h(i)*s_tub(i)*( Twin(i) - Taw(i) )
        qlr(i) = stef*eps*s_tub(i)*( Twin(i)**4 - T(i)**4 )
        qlt(i) = qlc(i) + qlr(i)
      
    end do
  
  end subroutine transf_parede

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

  subroutine grava_canal 
  
    integer :: unit = 99

    open ( unit, file = 'transf.mach.canal.bin', form='unformatted' )

!    write ( unit ) Cd_mach, F_star_mach, T(nxy-2*nx+1:nxy-nx), Twin, Taw, qlt, h, s_tub, 0.0d0 !, eps*stef
!    write ( unit ) Cd_mach, F_star_mach, T(nxy-2*nx+1:nxy-nx), Twin, qlt, h, s_tub, 0.0d0 ! eps*stef
    
    write ( unit ) eps*stef, Cd_mach, F_star_mach, Taw, T(1:nx), Twin, h, s_tub, qlt
        
    close ( unit )

  end subroutine grava_canal 

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

  subroutine grava_rad 
  
    integer :: unit = 99

    open ( unit, file = 'transf.mach.canal.bin', form='unformatted' )

!    write ( unit ) Cd_mach, F_star_mach, T(nxy-2*nx+1:nxy-nx), Twin, Taw, qlt, h, s_tub, 0.0d0 !, eps*stef
!    write ( unit ) Cd_mach, F_star_mach, T(nxy-2*nx+1:nxy-nx), Twin, qlt, h, s_tub, 0.0d0 ! eps*stef
    
    write ( unit ) eps*stef, Cd_mach, F_star_mach, Taw, T(1:nx), Twin, h, s_tub, qlt
        
    close ( unit )

  end subroutine grava_rad 

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

end module usuario
