module usuario

! Objetivos: definir a estimativa inicial das incógnitas,
!            aplicar as condições de contorno,
!            definir funções para propriedades,
!            calcular a função de corrente,
!            escrever resultados,
!            apresentar gráficos e
!            listar resultados.

! Versão atual:  11 Out 06

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

use coeficientes

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

contains

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

  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_cte, Rg_cte, p0, T0, Sg, fm1D )

    ! zeragem dos campos
    u    = 0.0d0
    v    = 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
    roe  = 0.0d0
    ron  = 0.0d0
    cup  = 0.0d0
    cvp  = 0.0d0
    plin = 0.0d0
    cp   = 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_cte / (gama_cte - 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_cte, Si/Sg, ind, M1D(i) )

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

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

       T1D(i) = T0 / cte

       u1D(i) = M1D(i) * dsqrt ( gama_cte * Rg_cte * 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)
  
       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_cte, Si/Sg, ind, Me1D )

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

       pe1D = p0 / (cte ** exp)

       Te1D = T0 / cte

       ue1D = Me1D * dsqrt ( gama_cte * Rg_cte * Te1D )

       do j = 2, ny-1

          np = (j-1)*nx + i

          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
    b = 0.0d0
    do np = 1, nxy
       a(np,3) = 1.0d0
    end do

  end subroutine estimativa_inicial

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

  subroutine EXATO_MACH_fluxo_massa ( gama_cte, Rg_cte, p_cam, T_cam, Sg, fm_exato )

    real*8,intent(in)  :: gama_cte  ! razão de calores específicos (adimensional)
    real*8,intent(in)  :: Rg_cte    ! constante do gás (J/kg.K)
    real*8,intent(in)  :: p_cam     ! pressão na câmara de combustão (Pa)
    real*8,intent(in)  :: 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_cte / ( Rg_cte * T_cam )

    a2 = 2.0d0 / ( gama_cte + 1.0d0 )

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

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

  end subroutine EXATO_MACH_fluxo_massa

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

  subroutine EXATO_MACH_newton ( gama_cte, a, ind, Mp )

    real*8,intent(in)   :: a        ! razão de áreas (adimensional)
    real*8,intent(in)   :: gama_cte ! 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_cte - 1.0d0) / 2.0d0
    c = (gama_cte + 1.0d0) / ( 2.0d0 - 2.0D0 * gama_cte)
    d = ((gama_cte + 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 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
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) =-1.0d0
       a(np,5) = 0.0d0
       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
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       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
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) =-1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    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
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 1.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       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
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do
  
  end subroutine ccv
!-------------------------------------------------

  subroutine ccT

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

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 1.0d0
       a(np,5) = 0.0d0
       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
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       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
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) =-1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 1, nx
       np = (j-1)*nx + i
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do
  
  end subroutine ccT

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

  subroutine ccp

    integer :: npn, npnn, nps, npss, npw,  npww, npe
    real*8  :: dyf, dyr,  ys,  xw,   dxf,  dxr

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 1.0d0
       a(np,5) = 0.0d0
       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
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       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
       npn  = np + nx
       npnn = np + 2*nx
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) =-1.0d0
       dyf = 2 * yp(npn)
       dyr = yp(npnn) - yp(npn)
       b(np)   = (dyf/dyr) * (pl(npn) - pl(npnn))
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       nps  = np - nx
       npss = np - 2*nx
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       ys  = (y(nps)+y(nps-1))/2
       dyf = 2 * (ys - yp(nps))
       dyr = yp(nps) - yp(npss)
       b(np)   = (dyf/dyr) * (pl(nps) - pl(npss))
    end do
  
  end subroutine ccp

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

  subroutine ccYi

    integer :: npn, npnn, nps, npss, npw,  npww, npe, npee
    real*8  :: dyf, dyr,  ys,  xw,   dxf,  dxr

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       npe  = np + 1
       npee = np + 2
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 1.0d0
       a(np,5) = 0.0d0
       dxf = 2 * xp(npe)
       dxr = xp(npee) - xp(npe)
!       b(np)   = (dxf/dxr) * (Yi(npe,ii) - Yi(npee,ii))
       b(np)   = 2 * Yi_in(j,ii)
	end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (Yi(npw,ii) - Yi(npww,ii))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np   = (j-1)*nx + i
       npn  = np + nx
       npnn = np + 2*nx
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) =-1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       nps  = np - nx
       npss = np - 2*nx
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       ys  = (y(nps)+y(nps-1))/2
       dyf = 2 * (ys - yp(nps))
       dyr = yp(nps) - yp(npss)
       b(np)   = (dyf/dyr) * (Yi(nps,ii) - Yi(npss,ii))
    end do
  
  end subroutine ccYi

!-------------------------------------------------
  
  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) = ( v(np) + v(npe) ) / 2
       Tin(j) = T0 - ( gama_cte - 1 ) * ( uin(j)**2 + vin(j)**2 ) / ( 2 * gama_cte * Rg_cte )
       Mw(j)  = dsqrt ( (uin(j)**2 + vin(j)**2) / (gama_cte * Rg_cte * Tin(j)) )
       aux    = 1 + ( gama_cte - 1 ) * ( Mw(j) ** 2 ) / 2
       pin(j) = p0 / ( aux ** ( gama_cte / ( gama_cte - 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 pressao_e_massa_especifica_com_pl

    ! correção da pressão com pl

    p = p + pl

    ! correção da massa especifica nodal com pl

    ro = p / ( Rg * T )
  
  end subroutine pressao_e_massa_especifica_com_pl

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

  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
       u(np) = u(nps)
       v(np) = v(nps)
    end do
  
  end subroutine u_v_nos_ficticios_com_pl

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

  subroutine Uce_Vcn_de_dn_faces_contornos

    integer :: npw, npe
    real*8  :: uw,  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
       uw  = ( u(np) + u(npw) ) / 2
       Uce(np) = uw * 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 calor_especifico
    ! Utilizado apenas em caso de monogás

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

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

	select case (tipo_cp)
	
	! Utilização de cp constante
	case (1) 
	
	cp = cp_cte

	! Propriedades (variáveis) para ar
	case (2)

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

				! Polinômio obtido a partir dos dados de Incropera (p. 458)
				! Dados para ar, para temperaturas entre 100 e 3000 K.
				cp(np) = 1061.33426d0 - 0.43283*T(np) + 0.00102d0*T(np)**2 - &
						 6.47481d-7*T(np)**3 + 1.38642d-10*T(np)**4

			end do
		end do

	! Propriedades (variáveis) para H2O
	case (3)

		open (501, file = 'propri1.dat') 
		read(501,*)
		read(501,*)   
		read(501,*)   
		read(501,*)
				
		do l = 1, 2
		   do m = 1, 5
			  read(501,*) a_cp(m,l)
		   end do
		   read(501,*)
		   read(501,*)
		   read(501,*)
		end do

		close(501)

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

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

				cp(np) = (a_cp(1,l) + a_cp(2,l)*T(np) + a_cp(3,l)*(T(np)**2)  &
					     +  a_cp(4,l)*(T(np)**3) + a_cp(5,l)*(T(np)**4)) * Rg_cte

			  end do
		end do	
	
	end select

  end subroutine calor_especifico

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

  subroutine chain 

    real*8 :: dpsi

    psi = 0.0d0

    do i = 1, nx-1

       do j = 2, ny-1

          np = (j-1)*nx + i

          dpsi = 2 * pi * re(np) * roe(np) * Uce(np)

          psi(np) = psi(np-nx) + dpsi

       end do

    end do

  end subroutine chain 

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

  subroutine globais

    real*8 :: dfm, ue

    i   = 1
    fmi = 0.0d0
    do j = 2, ny-1
       np = (j-1)*nx + i
       fmi = fmi + re(np) * roe(np) * Uce(np)
    end do
    fmi = fmi * 2 * pi

    i   = nx-1
    fme = 0.0d0
    Fd  = 0.0d0
    do j = 2, ny-1
       np = (j-1)*nx + i
       dfm = 2 * pi * re(np) * roe(np) * Uce(np)
       fme = fme + dfm
       ue  = (u(np) + u(np+1)) / 2
       Fd  = Fd  + dfm * ue
    end do

    write(8,1) it, fmi/fme, fme/fm1D, Fd/Fd1D

    1 format(i8,2x,3(1pe24.15))

  end subroutine globais

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

  subroutine resultados_principais

    real*8  :: Fp1D, F1D, Fv1D, Fo, ce1D, c1D, cv1D, CF1D, CFv1D, Is1D, Isv1D, Sex
    real*8  :: Fp,   F  , Fv  ,     ce  , c  , cv  , CF  , CFv  , Is  , Isv,   Fpv  
    real*8  :: CdKL, KL1, KL2, KL3, KL4
    integer :: npe
	integer :: ii

    ! soluções analíticas

    i   = nx-1
    j   = ny-1
    np  = (j-1)*nx + i
    Sex = pi * ( y(np) ** 2 )

    Fp1D  = Fpv1D - pr*Sex

    F1D   = Fd1D + Fp1D

    Fv1D  = Fd1D + Fpv1D

    Fo    = p0 * Sg

    CF1D  = F1D / Fo

    CFv1D = Fv1D / Fo

    ce1D  = Fo / fm1D

    c1D   = F1D / fm1D

    cv1D  = Fv1D / fm1D

    Is1D  = c1D / go

    Isv1D = cv1D / go

    write(10,1) fm1D, Fd1D, Fp1D, Fpv1D, F1D, Fv1D, Fo, CF1D, CFv1D, ce1D, &
                c1D,  cv1D, Is1D, Isv1D
 
    1 format(/,1x,'*** Solução analítica do escoamento Q1D isentrópico ***',//, &
                  1pe25.15,' = fm1D:  fluxo de massa (kg/s)',/, &
                  1pe25.15,' = Fd1D:  empuxo dinâmico (N)',/, &
                  1pe25.15,' = Fp1D:  empuxo de pressão ao nível do mar (p = 101325 Pa) (N)',/, &
                  1pe25.15,' = Fpv1D: empuxo de pressão no vácuo (N)',/, &
                  1pe25.15,' = F1D:   empuxo total ao nível do mar (N)',/, &
                  1pe25.15,' = Fv1D:  empuxo total no vácuo (N)',/, &
                  1pe25.15,' = Fo:    empuxo padrão (N)',/, &
                  1pe25.15,' = CF1D:  coeficiente de empuxo ao nível do mar (adimensional)',/, &
                  1pe25.15,' = CFv1D: coeficiente de empuxo no vácuo (adimensional)',/, &
                  1pe25.15,' = ce1D:  velocidade característica (m/s)',/, &
                  1pe25.15,' = c1D:   velocidade de ejeção efetiva ao nível do mar (m/s)',/, &
                  1pe25.15,' = cv1D:  velocidade de ejeção efetiva no vácuo (m/s)',/, &
                  1pe25.15,' = Is1D:  impulso específico ao nível do mar (s)',/, &
                  1pe25.15,' = Isv1D: impulso específico no vácuo (s)')

    KL1 = ( gama_cte + 1 ) / ( ( 1 + rcag ) ** 2 )
    KL2 = ( 8 * gama_cte - 27 ) / ( 2304 * ( 1 + rcag ) )
    KL3 = 754 * ( gama_cte ** 2 ) - ( 757 * gama_cte ) + 3633
    KL4 = 276480 * ( ( 1 + rcag ) ** 2 )

    CdKL = 1 - KL1 * ( ( 1 / 9.6d+1 ) - KL2 + ( KL3 / KL4 ) )

    write(10,2) CdKL
 
    2 format(/,1x,'*** Solução analítica 2D ***',//,  &
                  1pe25.15,' = CdKL: coeficiente de descarga de Kliegel e Levine (adimensional)')

    ! soluções numéricas

    i   = nx-1
    Fd  = 0.0d0
    Fp  = 0.0d0
    Fpv = 0.0d0
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       Fd  = Fd  + re(np) * roe(np) * Uce(np) * u(npe)
       Fp  = Fp  + (p(npe) - pr) * re(np) * ye(np)
       Fpv = Fpv + p(npe) * re(np) * ye(np)
    end do
    Fd  = Fd  * 2 * pi
    Fp  = Fp  * 2 * pi
    Fpv = Fpv * 2 * pi

    F   = Fd + Fp

    Fv  = Fd + Fpv

    CF  = F / Fo

    CFv = Fv / Fo

    ce  = Fo / fme

    c   = F / fme

    cv  = Fv / fme

    Is  = c / go

    Isv = cv / go

    write(10,3) fme, Fd, Fp, Fpv, F, Fv, Fo, CF, CFv, ce, c, cv, Is, Isv
 
    3 format(/,1x,'*** Solução numérica 2D ***',//, &
                  1pe25.15,' = fme: fluxo de massa (kg/s)',/, &
                  1pe25.15,' = Fd:  empuxo dinâmico (N)',/, &
                  1pe25.15,' = Fp:  empuxo de pressão ao nível do mar (p = 101325 Pa) (N)',/, &
                  1pe25.15,' = Fpv: empuxo de pressão no vácuo (N)',/, &
                  1pe25.15,' = F:   empuxo total ao nível do mar (N)',/, &
                  1pe25.15,' = Fv:  empuxo total no vácuo (N)',/, &
                  1pe25.15,' = Fo:  empuxo padrão (N)',/, &
                  1pe25.15,' = CF:  coeficiente de empuxo ao nível do mar (adimensional)',/, &
                  1pe25.15,' = CFv: coeficiente de empuxo no vácuo (adimensional)',/, &
                  1pe25.15,' = ce:  velocidade característica (m/s)',/, &
                  1pe25.15,' = c:   velocidade de ejeção efetiva ao nível do mar (m/s)',/, &
                  1pe25.15,' = cv:  velocidade de ejeção efetiva no vácuo (m/s)',/, &
                  1pe25.15,' = Is:  impulso específico ao nível do mar (s)',/, &
                  1pe25.15,' = Isv: impulso específico no vácuo (s)')

    write(10,4) fme/fm1D, Fd/Fd1D,   Fp/Fp1D, Fpv/Fpv1D, F/F1D,   Fv/Fv1D, &
                CF/CF1D,  CFv/CFv1D, ce/ce1D, c/c1D,     cv/cv1D, Is/Is1D, Isv/Isv1D
 
    4 format(/,1x,'*** Eficiências: solução numérica 2D / analítica Q1D (adimensionais) ***',//, &
                  1pe25.15,' = coeficiente de descarga',/, &
                  1pe25.15,' = empuxo dinâmico',/, &
                  1pe25.15,' = empuxo de pressão ao nível do mar (p = 101325 Pa)',/, &
                  1pe25.15,' = empuxo de pressão no vácuo',/, &
                  1pe25.15,' = empuxo total ao nível do mar',/, &
                  1pe25.15,' = empuxo total no vácuo',/, &
                  1pe25.15,' = coeficiente de empuxo ao nível do mar',/, &
                  1pe25.15,' = coeficiente de empuxo no vácuo',/, &
                  1pe25.15,' = velocidade característica',/, &
                  1pe25.15,' = velocidade de ejeção efetiva ao nível do mar',/, &
                  1pe25.15,' = velocidade de ejeção efetiva no vácuo',/, &
                  1pe25.15,' = impulso específico ao nível do mar',/, &
                  1pe25.15,' = impulso específico no vácuo')

  end subroutine resultados_principais

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

  subroutine curva_residuo

    integer :: k ! auxiliar

    nome = 'Mach2D_6p0_residuo.gnu'

    open(7,file=trim(adjustl(nome)))

    do k = 1, 6
       read(7,*)
    end do

    write(7,1) trim(adjustl(caso)) // '.residuo.dat'
    1 format("plot '", a, "' using 1:2 notitle")

    write(7,2) trim(adjustl(title))
    2 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_residuo.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine curva_residuo

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

  subroutine contornos_pos_processamento

    integer :: npn, npw, nps, npnn, npww
    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
       npnn = np + 2*nx
       raz  = ( yp(npnn) - yp(npn) ) / yp(npn)
       u(np) = ( ((1+raz)**2)*u(npn) - u(npnn) ) / ( raz*(2+raz) )
       v(np) = 0.0d0
       p(np) = ( p(npn) + p(np) ) / 2
       T(np) = ( ((1+raz)**2)*T(npn) - T(npnn) ) / ( raz*(2+raz) )
       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
       npnn = np + 2*nx
       raz  = ( y(npnn) + y(npn) ) / y(npn)
       u(np) = ( ((1+raz)**2)*u(npn) - u(npnn) ) / ( raz*(2+raz) )
       v(np) = 0.0d0
       p(np) = ( ((1+raz)**2)*p(npn) - p(npnn) ) / ( raz*(2+raz) )
       T(np) = ( ((1+raz)**2)*T(npn) - T(npnn) ) / ( raz*(2+raz) )
       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
       npnn = np + 2*nx
       raz  = ( y(npnn-1) + y(npn-1) ) / y(npn-1)
       u(np) = ( ((1+raz)**2)*u(npn) - u(npnn) ) / ( raz*(2+raz) )
       v(np) = 0.0d0
       p(np) = ( ((1+raz)**2)*p(npn) - p(npnn) ) / ( raz*(2+raz) )
       T(np) = ( ((1+raz)**2)*T(npn) - T(npnn) ) / ( raz*(2+raz) )
       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(nps)
       v(np) = v(nps)
       p(np) = ( p(nps) + p(np) ) / 2
       T(np) = T(nps)
       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
	   npww = np - 2
       if ( i == 1 ) then
	      p(np)  = p(nps)
          xp(np) = x(nps)
          yp(np) = y(nps)
       end if
       if ( i == nx ) then
	      p(np) = 2.0*p(npw) - p(npww)
          xp(np) = x(nps-1)
          yp(np) = y(nps-1)
       end if
    end do

  end subroutine contornos_pos_processamento

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

  subroutine contornos_propriedades_pos_processamento

   	integer :: ii   ! auxiliar
	integer :: np   ! volume P 
	integer :: npe, npee  
	integer :: npw, npww  
	integer :: npn, npnn
	integer :: nps, npss
	real*8  :: dxf, dxr
	real*8  :: xw
	real*8  :: raz


	! Frações mássicas das espécies
	do ii = 1, Ne

	   ! Atualização na entrada do domínio
	   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)
		   Yi(np,ii) = Yi(npe,ii) + (dxf/dxr) * (Yi(npe,ii) - Yi(npee,ii))
	   end do

	   ! Atualização na saída do domínio
	   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)
		   Yi(np,ii) = Yi(npw,ii) + (dxf/dxr) * (Yi(npw,ii) - Yi(npww,ii))
	   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
          npnn = np + 2*nx
          raz  = ( yp(npnn) - yp(npn) ) / yp(npn)
          Yi(np,ii) = ( ((1+raz)**2)*Yi(npn,ii) - Yi(npnn,ii) ) / ( raz*(2+raz) )
       end do

       ! canto esquerdo inferior
	   j = 1
	   i = 1
	   np   = (j-1)*nx + i
	   npn  = np + nx
	   npnn = np + 2*nx
	   raz  = ( y(npnn) + y(npn) ) / y(npn)
	   Yi(np,ii) = ( ((1+raz)**2)*Yi(npn,ii) - Yi(npnn,ii) ) / ( raz*(2+raz) )

	   ! canto direito inferior
	   j = 1
	   i = nx
	   np   = (j-1)*nx + i
	   npw  = np - 1
	   npn  = np + nx
	   npnn = np + 2*nx
	   raz  = ( y(npnn-1) + y(npn-1) ) / y(npn-1)
	   Yi(np,ii) = ( ((1+raz)**2)*Yi(npn,ii) - Yi(npnn,ii) ) / ( raz*(2+raz) )

       ! contorno superior (parede)
       j = ny
       do i = 1, nx
          np  = (j-1)*nx + i
          nps = np - nx
          Yi(np,ii) = Yi(nps,ii)
       end do

      ! cantos superiores
      j = ny
      do i = 1, nx, nx-1
         np  = (j-1)*nx + i
         nps = np - nx
         Yi(np,ii) = Yi(nps,ii)
      end do


	end do

    ! Atualização na entrada do domínio (extrapolação linear)
    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)
	   cp(np) = cp(npe) + (dxf/dxr) * (cp(npe) - cp(npee))
	   Rg(np) = Rg(npe) + (dxf/dxr) * (Rg(npe) - Rg(npee))
	   gama(np) = gama(npe) + (dxf/dxr) * (gama(npe) - gama(npee))
	   gamaf(np) = gamaf(npe) + (dxf/dxr) * (gamaf(npe) - gamaf(npee))
	   ent(np) = ent(npe) + (dxf/dxr) * (ent(npe) - ent(npee))
    end do

    ! Atualização na saída do domínio (extrapolação linear)
    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)
	   cp(np) = cp(npw) + (dxf/dxr) * (cp(npw) - cp(npww))
	   Rg(np) = Rg(npw) + (dxf/dxr) * (Rg(npw) - Rg(npww))
	   gama(np) = gama(npw) + (dxf/dxr) * (gama(npw) - gama(npww))
	   gamaf(np) = gamaf(npw) + (dxf/dxr) * (gamaf(npw) - gamaf(npww))
	   ent(np) = ent(npw) + (dxf/dxr) * (ent(npw) - ent(npww))
    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
       npnn = np + 2*nx
       raz  = ( yp(npnn) - yp(npn) ) / yp(npn)
       cp(np) = ( ((1+raz)**2)*cp(npn) - cp(npnn) ) / ( raz*(2+raz) )
       Rg(np) = ( ((1+raz)**2)*Rg(npn) - Rg(npnn) ) / ( raz*(2+raz) )
       gama(np) = ( ((1+raz)**2)*gama(npn) - gama(npnn) ) / ( raz*(2+raz) )
       gamaf(np) = ( ((1+raz)**2)*gamaf(npn) - gamaf(npnn) ) / ( raz*(2+raz) )
       ent(np) = ( ((1+raz)**2)*ent(npn) - ent(npnn) ) / ( raz*(2+raz) )
    end do

   ! canto esquerdo inferior
   j = 1
   i = 1
   np   = (j-1)*nx + i
   npn  = np + nx
   npnn = np + 2*nx
   raz  = ( y(npnn) + y(npn) ) / y(npn)
   cp(np) = ( ((1+raz)**2)*cp(npn) - cp(npnn) ) / ( raz*(2+raz) )
   Rg(np) = ( ((1+raz)**2)*Rg(npn) - Rg(npnn) ) / ( raz*(2+raz) )
   gama(np) = ( ((1+raz)**2)*gama(npn) - gama(npnn) ) / ( raz*(2+raz) )
   gamaf(np) = ( ((1+raz)**2)*gamaf(npn) - gamaf(npnn) ) / ( raz*(2+raz) )
   ent(np) = ( ((1+raz)**2)*ent(npn) - ent(npnn) ) / ( raz*(2+raz) )

   ! canto direito inferior
   j = 1
   i = nx
   np   = (j-1)*nx + i
   npw  = np - 1
   npn  = np + nx
   npnn = np + 2*nx
   raz  = ( y(npnn-1) + y(npn-1) ) / y(npn-1)
   cp(np) = ( ((1+raz)**2)*cp(npn) - cp(npnn) ) / ( raz*(2+raz) )
   Rg(np) = ( ((1+raz)**2)*Rg(npn) - Rg(npnn) ) / ( raz*(2+raz) )
   gama(np) = ( ((1+raz)**2)*gama(npn) - gama(npnn) ) / ( raz*(2+raz) )
   gamaf(np) = ( ((1+raz)**2)*gamaf(npn) - gamaf(npnn) ) / ( raz*(2+raz) )
   ent(np) = ( ((1+raz)**2)*ent(npn) - ent(npnn) ) / ( raz*(2+raz) )

   ! contorno superior (parede)
   j = ny
   do i = 1, nx
      np  = (j-1)*nx + i
      nps = np - nx
      cp(np) = cp(nps)
      Rg(np) = Rg(nps)
      gama(np) = gama(nps)
      gamaf(np) = gamaf(nps)
      ent(np) = ent(nps)
   end do

   ! cantos superiores
   j = ny
   do i = 1, nx, nx-1
      np  = (j-1)*nx + i
      nps = np - nx
      cp(np) = cp(nps)
      Rg(np) = Rg(nps)
      gama(np) = gama(nps)
      gamaf(np) = gamaf(nps)
      ent(np) = ent(nps)
   end do


  end subroutine contornos_propriedades_pos_processamento

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

  subroutine grafico_x_velocidade_u

    integer :: k, nps, npn

    nome = trim(adjustl(caso)) // '.velocidade_u.dat'

    open(7,file=trim(adjustl(nome)))

    write(7,1)
    1 format('### VELOCIDADE u (m/s) ###', /, &
             '#', t4,'x (m)', t28,'1D', t52,'centro', t76,'parede')

    do i = 1, nx
       j = 1
       nps = (j-1)*nx + i
       j = ny
       npn = (j-1)*nx + i
       if ( i==1 .or. i==nx .or. mod(i,w_g)==0 ) &
          write(7,2) xp(nps), u1D(i), u(nps), u(npn)
       2 format(4(1pe24.15))
    end do

    close(7)

    open(7,file= 'Mach2D_6p0_velocidade_u.gnu' )

    do k = 1, 6
       read(7,*)
    end do

    write(7,3) trim(adjustl(nome))
    3 format("plot '", a, "' using 1:2 title '1D'")

    write(7,4) trim(adjustl(nome))
    4 format("replot '", a, "' using 1:3 title 'centro'")

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:4 title 'parede'")

    write(7,7) trim(adjustl(title))
    7 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_velocidade_u.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_x_velocidade_u

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

  subroutine grafico_x_velocidade_v

    integer :: k, nps, npn

    nome = trim(adjustl(caso)) // '.velocidade_v.dat'

    open(7,file=trim(adjustl(nome)))

    write(7,1)
    1 format('### VELOCIDADE v (m/s) ###', /, &
             '#', t4,'x (m)', t28,'1D', t52,'centro', t76,'parede')

    do i = 1, nx
       j = 1
       nps = (j-1)*nx + i
       j = ny
       npn = (j-1)*nx + i
       if ( i==1 .or. i==nx .or. mod(i,w_g)==0 ) &
          write(7,2) xp(nps), 0.0d0, v(nps), v(npn)
       2 format(4(1pe24.15))
    end do

    close(7)

    open(7,file= 'Mach2D_6p0_velocidade_v.gnu' )

    do k = 1, 6
       read(7,*)
    end do

    write(7,3) trim(adjustl(nome))
    3 format("plot '", a, "' using 1:2 title '1D'")

    write(7,4) trim(adjustl(nome))
    4 format("replot '", a, "' using 1:3 title 'centro'")

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:4 title 'parede'")

    write(7,6) trim(adjustl(title))
    6 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_velocidade_v.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_x_velocidade_v

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

  subroutine grafico_x_pressao

    integer :: k, nps, npn

    nome = trim(adjustl(caso)) // '.pressao.dat'

    open(7,file=trim(adjustl(nome)))

    write(7,1)
    1 format('### PRESSÃO (Pa) ###', /, &
             '#', t4,'x (m)', t28,'1D', t52,'centro', t76,'parede')

    do i = 1, nx
       j = 1
       nps = (j-1)*nx + i
       j = ny
       npn = (j-1)*nx + i
       if ( i==1 .or. i==nx .or. mod(i,w_g)==0 ) &
          write(7,2) xp(nps), p1D(i), p(nps), p(npn)
       2 format(4(1pe24.15))
    end do

    close(7)

    open(7,file= 'Mach2D_6p0_pressao.gnu' )

    do k = 1, 7
       read(7,*)
    end do

    write(7,3) trim(adjustl(nome))
    3 format("plot '", a, "' using 1:2 title '1D'")

    write(7,4) trim(adjustl(nome))
    4 format("replot '", a, "' using 1:3 title 'centro'")

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:4 title 'parede'")

    write(7,6) trim(adjustl(title))
    6 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_pressao.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_x_pressao

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

  subroutine grafico_x_temperatura

    integer :: k, nps, npn

    nome = trim(adjustl(caso)) // '.temperatura.dat'

    open(7,file=trim(adjustl(nome)))

    write(7,1)
    1 format('### TEMPERATURA (K) ###', /, &
             '#', t4,'x (m)', t28,'1D', t52,'centro', t76,'parede')

    do i = 1, nx
       j = 1
       nps = (j-1)*nx + i
       j = ny
       npn = (j-1)*nx + i
       if ( i==1 .or. i==nx .or. mod(i,w_g)==0 ) &
          write(7,2) xp(nps), T1D(i), T(nps), T(npn)
       2 format(4(1pe24.15))
    end do

    close(7)

    open(7,file= 'Mach2D_6p0_temperatura.gnu' )

    do k = 1, 6
       read(7,*)
    end do

    write(7,3) trim(adjustl(nome))
    3 format("plot '", a, "' using 1:2 title '1D'")

    write(7,4) trim(adjustl(nome))
    4 format("replot '", a, "' using 1:3 title 'centro'")

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:4 title 'parede'")

    write(7,6) trim(adjustl(title))
    6 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_temperatura.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_x_temperatura

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

  subroutine grafico_x_massa_especifica

    integer :: k, nps, npn

    nome = trim(adjustl(caso)) // '.massa_especifica.dat'

    open(7,file=trim(adjustl(nome)))

    write(7,1)
    1 format('### MASSA ESPECÍFICA (kg/m3) ###', /, &
             '#', t4,'x (m)', t28,'1D', t52,'centro', t76,'parede')

    do i = 1, nx
       j = 1
       nps = (j-1)*nx + i
       j = ny
       npn = (j-1)*nx + i
       if ( i==1 .or. i==nx .or. mod(i,w_g)==0 ) &
          write(7,2) xp(nps), p1D(i)/(Rg(i)*T1D(i)), ro(nps), ro(npn)
       2 format(4(1pe24.15))
    end do

    close(7)

    open(7,file= 'Mach2D_6p0_massa_especifica.gnu' )

    do k = 1, 6
       read(7,*)
    end do

    write(7,3) trim(adjustl(nome))
    3 format("plot '", a, "' using 1:2 title '1D'")

    write(7,4) trim(adjustl(nome))
    4 format("replot '", a, "' using 1:3 title 'centro'")

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:4 title 'parede'")

    write(7,6) trim(adjustl(title))
    6 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_massa_especifica.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_x_massa_especifica

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

  subroutine grafico_x_Mach

    integer :: k, nps, npn

    nome = trim(adjustl(caso)) // '.Mach.dat'

    open(7,file=trim(adjustl(nome)))

    write(7,1)
    1 format('### MACH ###', /, &
             '#', t4,'x (m)', t28,'1D', t52,'centro', t76,'parede')

    do i = 1, nx
       j = 1
       nps = (j-1)*nx + i
       j = ny
       npn = (j-1)*nx + i
       if ( i==1 .or. i==nx .or. mod(i,w_g)==0 ) &
          write(7,2) xp(nps), M1D(i), M(nps), M(npn)
       2 format(4(1pe24.15))
    end do

    close(7)

    open(7,file= 'Mach2D_6p0_Mach.gnu' )

    do k = 1, 6
       read(7,*)
    end do

    write(7,3) trim(adjustl(nome))
    3 format("plot '", a, "' using 1:2 title '1D'")

    write(7,4) trim(adjustl(nome))
    4 format("replot '", a, "' using 1:3 title 'centro'")

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:4 title 'parede'")

    write(7,6) trim(adjustl(title))
    6 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_Mach.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_x_Mach

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

  subroutine grafico_x_empuxo

    integer :: k,  nps
    real*8  :: Fd, Fp,  Fpv, Fv, F, rpp, dy

    nome = trim(adjustl(caso)) // '.empuxo.dat'

    open(7,file=trim(adjustl(nome)))

    write(7,1)
    1 format('### EMPUXO (N) ###', /, &
             '#', t4,'x (m)', t18,'dinâmico', t32,'pressão-mar', &
             t46,'pressão-vácuo', t60,'total-mar', t74,'total-vácuo' )

    do i = ig+1, nx
       Fd  = 0.0d0
       Fp  = 0.0d0
       Fpv = 0.0d0
       do j = 2, ny-1
          np  = (j-1)*nx + i
          rpp = rp(np)
          dy  = (ye(np-1)+ye(np))/2
		  if ( i == nx ) then
             rpp = re(np-1)
             dy  = ye(np-1)
          end if
          Fd  = Fd  + rpp * ro(np) * dy * (u(np)**2)
          Fp  = Fp  + (p(np) - pr) * rpp * dy
          Fpv = Fpv +  p(np)       * rpp * dy
       end do
       Fd  = Fd  * 2 * pi
       Fp  = Fp  * 2 * pi
       Fpv = Fpv * 2 * pi
       F   = Fd + Fp
       Fv  = Fd + Fpv
       j = 1
       nps = (j-1)*nx + i
       if ( i==ig+1 .or. i==nx .or. mod(i,w_g)==0 ) &
          write(7,2) xp(nps), Fd, Fp, Fpv, F, Fv 
       2 format(6(1pe14.5))
    end do

    close(7)

    open(7,file= 'Mach2D_6p0_empuxo.gnu' )

    do k = 1, 6
       read(7,*)
    end do

    write(7,3) trim(adjustl(nome))
    3 format("plot '", a, "' using 1:2 title 'dinâmico'")

    write(7,4) trim(adjustl(nome))
    4 format("replot '", a, "' using 1:3 title 'pressão-mar'")

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:4 title 'pressão-vácuo'")

    write(7,6) trim(adjustl(nome))
    6 format("replot '", a, "' using 1:5 title 'total-mar'")

    write(7,7) trim(adjustl(nome))
    7 format("replot '", a, "' using 1:6 title 'total-vácuo'")

    write(7,8) trim(adjustl(title))
    8 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot Mach2D_6p0_empuxo.gnu'

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_x_empuxo

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

  subroutine pos_processamento

    integer :: k, nps, npn
    real*8,dimension(nx) :: fis ! variável no contorno sul
    real*8,dimension(nx) :: fin ! variável no contorno norte

	character(50) propriedade

    ! escreve o valor da norma L1 dos resíduos na iteração 1
    write(10,1) norma1
    1 format(/,5x,'Norma L1 dos resíduos na iteração 1 =',1pe15.6)

    ! calcula a função de corrente
    call chain

    ! escreve campos

    call contornos_pos_processamento
	if ((modelo_f == 3).or.(modelo_f == 4)) call contornos_propriedades_pos_processamento

	M = dsqrt ( (u**2 + v**2) / (gama * Rg * T) )
    ro = p / ( Rg * T )

    if ( w_cam == 1 ) then

       write(10,09)
       09 format(//,10x,'*** campos ***', //, &
                    10x,'na horizontal (esquerda para direita): contorno oeste para leste', //, &
                    10x,'na vertical   (de cima para baixo):    contorno sul   para norte'  )

       write(10,10)
       10 format(//,30x,'M: número de Mach')
       call print1(M,ny,nx,nxy)

       write(10,11)
       11 format(//,30x,'u: velocidade cartesiana u nodal')
       call print1(u,ny,nx,nxy)

       write(10,12)
       12 format(//,30x,'v: velocidade cartesiana v nodal')
       call print1(v,ny,nx,nxy)

       write(10,13)
       13 format(//,30x,'p: pressão')
       call print1(p,ny,nx,nxy)

       write(10,14)
       14 format(//,30x,'T: Temperatura')
       call print1(T,ny,nx,nxy)

       write(10,15)
       15 format(//,30x,'ro: massa específica nodal')
       call print1(ro,ny,nx,nxy)

       write(10,55)
       55 format(//,30x,'cp: calor específico à pressão constante')
       call print1(cp,ny,nx,nxy)

       write(10,16)
       16 format(//,30x,'Uce: velocidade contravariante U na face leste')
       call print1(Uce,ny,nx,nxy)

       write(10,17)
       17 format(//,30x,'Vcn: velocidade contravariante V na face norte')
       call print1(Vcn,ny,nx,nxy)

       write(10,18)
       18 format(//,30x,'pl: variação de pressão')
       call print1(pl,ny,nx,nxy)

       write(10,19)
       19 format(//,30x,'b: fonte de pl')
       call print1(b,ny,nx,nxy)

       write(10,20)
       20 format(//,30x,'psi: função de corrente')
       call print1(psi,ny,nx,nxy)

       write(10,21)
       21 format(//,30x,'X: coordenada X no canto nordeste')
       call print1(x,ny,nx,nxy)

       write(10,22)
       22 format(//,30x,'Y: coordenada Y no canto nordeste')
       call print1(y,ny,nx,nxy)

       write(10,23)
       23 format(//,30x,'Jp: Jacobiano no centro do volume')
       call print1(Jp,ny,nx,nxy)

       write(10,24)
       24 format(//,30x,'xe: métrica x eta na face leste')
       call print1(xe,ny,nx,nxy)

       write(10,25)
       25 format(//,30x,'ye: métrica y eta na face leste')
       call print1(ye,ny,nx,nxy)

       write(10,26)
       26 format(//,30x,'xk: métrica x ksi na face norte')
       call print1(xk,ny,nx,nxy)

       write(10,27)
       27 format(//,30x,'yk: métrica y ksi na face norte')
       call print1(yk,ny,nx,nxy)

       write(10,28)
       28 format(//,30x,'de: coeficiente d do SIMPLEC na face leste')
       call print1(de,ny,nx,nxy)

       write(10,29)
       29 format(//,30x,'dn: coeficiente d do SIMPLEC na face norte')
       call print1(dn,ny,nx,nxy)

    end if

    if ((modelo_f == 2).or.(modelo_f == 3).or.(modelo_f == 4)) then
		call escreve_campos_especies
		propriedade = 'Constante do gás - Rg (J/kg.K)'
		call escreve_campos_propriedades1 (Rg, propriedade)
		propriedade = 'Calor específico - cp (J/kg.K)'
		call escreve_campos_propriedades1 (cp, propriedade)
		propriedade = 'Gama (adim.)'
		call escreve_campos_propriedades1 (gama, propriedade)
		propriedade = 'Entropia (J/kg.K)'
		if ((modelo_f == 3).or.(modelo_f == 4)) then
			propriedade = 'Gama congelado (adim.)'
			call escreve_campos_propriedades1 (gama, propriedade)
		end if
		propriedade = 'Entropia (J/kg.K)'
		call escreve_campos_propriedades2 (ent, propriedade)

	end if

    call resultados_principais

    write(10,31) tcpu
    31 format(//,f14.2,' = tcpu: tempo de CPU (s)')

    ! faz o gráfico dos resíduos ao longo das iterações

    call curva_residuo

    ! mostra listagem do arquivo das iterações

    nome = 'notepad '// trim(adjustl(caso)) // '.itn'

    if ( sem_a == 0 ) dos = system(nome)

    ! fecha o arquivo principal de resultados

    close (10)

    ! mostra listagem do arquivo principal de resultados

    nome = 'notepad '// trim(adjustl(caso)) // '.txt'

    if ( sem_a == 0 ) dos = system(nome)

    ! *** GRÁFICOS ***

    call grafico_x_velocidade_u

    call grafico_x_velocidade_v

    call grafico_x_pressao

    call grafico_x_temperatura

    call grafico_x_massa_especifica

    call grafico_x_Mach

    call grafico_x_empuxo

  end subroutine pos_processamento

!-------------------------------------------------
  subroutine escreve_campos_especies

	integer :: ii  ! auxiliar
	character(80) :: aux   ! auxiliar

	! Escreve campos de espécies no centro da tubeira


	do ii = 1, Ne

	   aux = ' *** CENTRO *** ' // ' *** espécie ' // especie(ii) // ' ***'
	   call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
							   1,                       & ! freqüência de escrita
							   aux                   ,  & ! 50 colunas
							   ' ponto',                & !  6 colunas
							   'x (m)                ', & ! 21 colunas
							   'Yi (adim.)           ', & ! 21 colunas
							   1,                       & ! número do ponto inicial
							   nx,                      & ! número de dados a escrever
							   xp,                      & ! coordenada espacial
							   Yi(:,ii)                 ) ! variável a escrever
	end do

	! Escreve campos de espécies na parede da tubeira
	do ii = 1, Ne

	   aux = ' *** PAREDE *** ' // ' *** espécie ' // especie(ii) // ' ***'
	   call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
							   1,                       & ! freqüência de escrita
							   aux                   ,  & ! 50 colunas
							   ' ponto',                & !  6 colunas
							   'x (m)                ', & ! 21 colunas
							   'Yi (adim.)           ', & ! 21 colunas
							   (ny-1)*nx+1,             & ! número do ponto inicial
							   nxy,                     & ! número de dados a escrever
							   xp,                      & ! coordenada espacial
							   Yi(:,ii)                 ) ! variável a escrever
	end do


  end subroutine escreve_campos_especies

!-------------------------------------------------
  
  subroutine escreve_campos_propriedades ( T, propriedade )

	integer :: ii  !auxiliar
	character(80) :: aux   ! auxiliar
	
	real*8, intent(in), dimension(:) :: T
	character(50), intent(in) :: propriedade

	! Escreve campos de espécies no centro da tubeira

    aux = ' *** CENTRO *** ' // trim(adjustl(propriedade)) // ' ***'
    call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
 						   1,                       & ! freqüência de escrita
						   aux                   ,  & ! 50 colunas
						   ' ponto',                & !  6 colunas
						   'x (m)                ', & ! 21 colunas
						   'sol. numérica        ', & ! 21 colunas
						   1,                       & ! número do ponto inicial
						   nx,                      & ! número de dados a escrever
						   xp,                      & ! coordenada espacial
						   T                        ) ! variável a escrever

	! Escreve campos de espécies na parede da tubeira
    aux = ' *** PAREDE *** ' // trim(adjustl(propriedade)) // ' ***'
    call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
						   1,                       & ! freqüência de escrita
						   aux                   ,  & ! 50 colunas
						   ' ponto',                & !  6 colunas
						   'x (m)                ', & ! 21 colunas
						   'sol. numérica        ', & ! 21 colunas
						   (ny-1)*nx+1,             & ! número do ponto inicial
						   nxy,                     & ! número de dados a escrever
						   xp,                      & ! coordenada espacial
						   T                        ) ! variável a escrever

  end subroutine escreve_campos_propriedades

!-------------------------------------------------
  
  subroutine escreve_campos_propriedades1 ( T, propriedade )

	integer :: ii  !auxiliar
	character(80) :: aux   ! auxiliar
	
	real*8, intent(in), dimension(:) :: T
	character(50), intent(in) :: propriedade

	! Escreve campos de espécies no centro da tubeira

   aux = ' *** CENTRO *** ' // trim(adjustl(propriedade)) // ' ***'
   call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
						   1,                       & ! freqüência de escrita
						   aux                   ,  & ! 50 colunas
						   ' ponto',                & !  6 colunas
						   'x (m)                ', & ! 21 colunas
						   'sol. numérica        ', & ! 21 colunas
						   1,                       & ! número do ponto inicial
						   nx,                      & ! número de dados a escrever
						   xp,                      & ! coordenada espacial
						   T                        ) ! variável a escrever

	! Escreve campos de espécies na parede da tubeira
   aux = ' *** PAREDE *** ' // trim(adjustl(propriedade)) // ' ***'
   call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
						   1,                       & ! freqüência de escrita
						   aux                   ,  & ! 50 colunas
						   ' ponto',                & !  6 colunas
						   'x (m)                ', & ! 21 colunas
						   'sol. numérica        ', & ! 21 colunas
						   (ny-1)*nx+1,             & ! número do ponto inicial
						   nxy,                     & ! número de dados a escrever
						   xp,                      & ! coordenada espacial
						   T                        ) ! variável a escrever


  end subroutine escreve_campos_propriedades1

!-------------------------------------------------
  
  subroutine escreve_campos_propriedades2 ( T, propriedade )

	integer :: ii  !auxiliar
	character(80) :: aux   ! auxiliar
	
	real*8, intent(in), dimension(:) :: T
	character(50), intent(in) :: propriedade

	! Escreve campos de espécies no centro da tubeira
   aux = ' *** CENTRO *** ' // trim(adjustl(propriedade)) // ' ***'
   call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
						   1,                       & ! freqüência de escrita
						   aux                   ,  & ! 50 colunas
						   ' ponto',                & !  6 colunas
						   'x (m)                ', & ! 21 colunas
						   'sol. numérica        ', & ! 21 colunas
						   2,                       & ! número do ponto inicial
						   nx-1,                    & ! número de dados a escrever
						   xp,                      & ! coordenada espacial
						   T                        ) ! variável a escrever

	! Escreve campos de espécies na parede da tubeira
   aux = ' *** PAREDE *** ' // trim(adjustl(propriedade)) // ' ***'
   call ARQUIVO_escreve_campo_2D ( 10,              & ! unidade de arquivo para escrever
						   1,                       & ! freqüência de escrita
						   aux                   ,  & ! 50 colunas
						   ' ponto',                & !  6 colunas
						   'x (m)                ', & ! 21 colunas
						   'sol. numérica        ', & ! 21 colunas
						   (ny-1)*nx+2,             & ! número do ponto inicial
						   nxy-1,                   & ! número de dados a escrever
						   xp,                      & ! coordenada espacial
						   T                        ) ! variável a escrever


  end subroutine escreve_campos_propriedades2


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

  subroutine ARQUIVO_escreve_campo_2D ( 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. mod(i-1,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 ARQUIVO_escreve_campo_2D

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

end module usuario
