module resultados

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

! Objetivos: processar o algoritmo principal do código,
!            calcular a função de corrente,
!            escrever resultados,
!            apresentar gráficos e
!            listar resultados.

! Versão atual: 27 Jan 09

! Relação das sub-rotinas em ordem de apresentação:
!    1) solucao_numerica
!    2) pos_processamento
!    3) print1
!    4) grafico_x_velocidade_u
!    5) grafico_x_velocidade_v
!    6) grafico_x_pressao
!    7) grafico_x_temperatura
!    8) grafico_x_massa_especifica
!    9) grafico_x_Mach
!   10) grafico_x_empuxo
!   11) chain 
!   12) globais
!   13) resultados_principais
!   14) curva_residuo

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

use usuario
use solvers

implicit none

contains

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

  subroutine solucao_numerica

    ! leitura dos dados principais
    
	call le_dados

    ! criação de arquivos de resultados
    
    ! arquivo das iterações
  	open( 8,file=trim(adjustl(caso))//'.itn')

    ! arquivo principal
    open(10,file=trim(adjustl(caso))//'.txt')
    write(10,7) 
    7 format(/,'Listagem gerada pelo programa Mach2D 5.5.08, de 29 Jan 2009')

    ! arquivo dos resíduos
    open(11,file=trim(adjustl(caso))//'.residuo.dat')

    ! escrita dos dados principais

    call data_hora
    call mostra_dados

    ! alocacao de memória

    call aloca_memoria

    ! geração da malha e cálculo de parâmetros geométricos

    call define_dominio
    call gera_malha
    call pre_processamento
    call x_y_nodal
    call metricas
    call raios

    ! estimativa inicial das incógnitas

    call estimativa_inicial

    ! texto inicial nos arquivos de resultados

    write(8,1)
    1 format(t1,'Iteração', t14,'fmi / fme', t38,'fme / fm1D', t62,'Fd / Fd1D' )
    write(10,2) 
    2 format(//,t3,'Iteração',t14,'L1(R)n / L1(R)1', t38,'dt (s)', /)

    write(11,3) 
    3 format('# Iteração',t14,'L1(R)n / L1(R)1')

    ! zeragem do tempo de CPU

!    call cpu_time ( tcpu1 )
    tcpu1 = timef()

    ! ciclo iterativo externo no tempo

    do it = ito, itmax

       ! cálculo do dt de cada iteração
       if ( it <= It1 ) dt = dt1 
       if ( it >  It1 .and. it < It2 ) &
          dt = dt1 + (dt2-dt1)*(it-It1)/(It2-It1)
       if ( it >= It2 ) dt = dt2 

       ! zeragem da norma a cada iteração
       norma = 0.0d0

       ! atualização das incógnitas para novo incremento de tempo
       ua   = u
       va   = v
       uea  = ue
       vea  = ve
       una  = un
       vna  = vn
       pa   = p
       Ta   = T
       roa  = ro

       ! atualização das condições na entrada da tubeira
       pina = pin
       call Tin_pin_entrada
       call plin_entrada_p_ficticio

       ! cálculo de propriedades
       call calor_especifico

       ! solução de u
       call coeficientes_gerais
       call fonte_velocidade_u
       call ccu
       call lu2d5(a,dl,du,nxy,nx,ny)
       call fb2d5(a,dl,du,ny,nx,nxy,u,tolu,nitm_u,b)
       call norma_L1 (nx,ny,nxy,a,b,u,norma)

       ! solução de v
       call fonte_velocidade_v
       call ccv
       call lu2d5(a,dl,du,nxy,nx,ny)
       call fb2d5(a,dl,du,ny,nx,nxy,v,tolu,nitm_u,b)
       call norma_L1 (nx,ny,nxy,a,b,v,norma)

       ! solução de T
       call fonte_temperatura
       call ccT
       call lu2d5(a,dl,du,nxy,nx,ny)
       call fb2d5(a,dl,du,ny,nx,nxy,T,tolu,nitm_u,b)
       call norma_L1 (nx,ny,nxy,a,b,T,norma)

       ! cálculo da massa especifica
       call massa_especifica_nos
       call massa_especifica_faces

       ! cálculo dos coeficientes do SIMPLEC e das velocidades contravariantes nas faces
       call coeficientes_SIMPLEC
       call velocidades_contravariantes_faces
       call Uce_Vcn_de_dn_faces_contornos

       ! ciclo iterativo da massa
       
	   do im = 1, imax

          ! solução de pl
          call coeficientes_fonte_pressao
          call ccp
          call lu2d5(a,dl,du,nxy,nx,ny)
          call fb2d5(a,dl,du,ny,nx,nxy,pl,tolp,nitm_p,b)
       
          ! correção de p, u, v, Uce, Vcn, ro, roe e ron com pl
          call pressao_e_massa_especifica_com_pl
          call massa_especifica_faces
          call u_v_nos_reais_com_pl
          call u_v_nos_ficticios_com_pl
          call Uce_Vcn_faces_internas_com_pl
          call Uce_Vcn_faces_contornos_com_pl

       end do ! fim do ciclo iterativo da massa

       call norma_L1 (nx,ny,nxy,a,b,pl,norma)

       if ( it == 1 ) norma1 = norma

       if ( it==1 .or. it==itmax .or. norma/norma1<=tole .or. mod(it,w_r)==0 ) then
          write(10,4) it, norma/norma1, dt
          write(11,4) it, norma/norma1
          4 format(2x,i8,1pe15.6,9x,1pe15.6)
          call globais
       end if

       if ( norma/norma1 <= tole ) exit

    end do ! fim do ciclo iterativo externo no tempo

    if ( it > itmax ) it = itmax

    ! finalização do tempo de CPU

!    call cpu_time ( tcpu2 )
    tcpu2 = timef()

    tcpu = tcpu2-tcpu1 + tcpuo

    open(99,file=trim(adjustl(caso))//'.bin',form='unformatted')
    write(99) it,norma1,tcpu,pin,u,v,p,pl,ue,vn,T,ro,ve,un,Uce,Vcn,roe,ron
    close(99)

    close(8) ! fecha arquivo das iterações

    close(11) ! fecha arquivo dos resíduos

    ! executa o pós-processamento
  
    call pos_processamento

  end subroutine solucao_numerica

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

  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

    ! 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

    ! calcula e armazena valores dos contornos nos fictícios
    call contornos_pos_processamento

    ! calcula o campo do número de Mach
    M = dsqrt ( (u**2 + v**2) / (gama * Rg * T) )

    ! calcula o campo de massa especíica
    ro = p / ( Rg * T )

    ! escreve campos

    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 norte para sul' )

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    end if

    call resultados_principais

    write(10,31) tcpuo, tcpu2-tcpu1, tcpu
    31 format(//,f14.2,' = tcpuo: tempo de CPU (s) anterior acumulado deste caso', &
               /,f14.2,' = dtcpu: tempo de CPU (s) da simulação atual', &
               /,f14.2,' = tcpu:  tempo de CPU (s) total acumulado deste caso')

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

    ! mostra listagem do arquivo das iterações
    if ( sem_a == 0 ) dos = system('notepad '// trim(adjustl(caso)) // '.itn')

    ! fecha o arquivo principal de resultados
    close (10)

    ! mostra listagem do arquivo principal de resultados
    if ( sem_a == 0 ) dos = system('notepad '// trim(adjustl(caso)) // '.txt')

    ! lista arquivo das variáveis de interesse
    if ( sem_a == 0 ) dos = system('notepad '// trim(adjustl(caso)) // '.Richardson_3p0')

    ! *** gera e apresenta GRÁFICOS dos resultados ***

    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 print1 ( b, nj, ni, nij )

    integer,intent(in) :: nij, ni, nj
    real*8,intent(in),dimension(nij) :: b ! incógnita a imprimir
    real*8  :: vmin, vmax, xmin, xmax, ymin, ymax
    integer :: lmin(1), lmax(1)
    integer :: k, ll, npi, nps, l, ka

    k=0
    ll=10

    2 continue

    if ( k+10 <= ni ) then
       write(10,51) (ka,ka=k+1,k+10)
       51 format(/,'i=', 10(i10,5x))
    else
       write(10,51) (ka,ka=k+1,ni)
    end if

    do j=nj,1,-1
       np=(j-1)*ni + k
       npi=np+1
       nps=np+ll
       if(ni.lt.10) nps=np+ni
       write(10,10) j,(b(l),l=npi,nps)
    end do

    if(ni.lt.10) goto 3
    if(nps.eq.ni) goto 3

    k=k+10
    ll=ni-k
    if(ll.gt.10)ll=10

    goto 2

    10 format('j=',i4,10(1pe15.7))
  
    3 continue

    vmin = minval(b)
    lmin = minloc(b)
    xmin = xp(lmin(1))
    ymin = yp(lmin(1))

    vmax = maxval(b)
    lmax = maxloc(b)
    xmax = xp(lmax(1))
    ymax = yp(lmax(1))

    write(10,61) vmin, xmin, ymin
    61 format(/, 'Mínimo do campo:', 5x, 1pe15.7, ' = valor', 5x, &
                  1pe15.7, ' = x', 5x, 1pe15.7, ' = y') 

    write(10,62) vmax, xmax, ymax
    62 format(/, 'Máximo do campo:', 5x, 1pe15.7, ' = valor', 5x, &
                  1pe15.7, ' = x', 5x, 1pe15.7, ' = y') 

  end subroutine print1

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

  subroutine grafico_x_velocidade_u

    integer :: k, nps, npn

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

!    open(7,file=trim(adjustl(nome)))
    open(7,file=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_5p5_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_5p5_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)))
    open(7,file=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_5p5_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_5p5_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)))
    open(7,file=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_5p5_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_5p5_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)))
    open(7,file=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_5p5_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_5p5_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)))
    open(7,file=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*T1D(i)), ro(nps), ro(npn)
       2 format(4(1pe24.15))
    end do

    close(7)

    open(7,file= 'Mach2D_5p5_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_5p5_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)))
    open(7,file=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_5p5_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_5p5_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)))
    open(7,file=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_5p5_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_5p5_empuxo.gnu'

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

  end subroutine grafico_x_empuxo

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

  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, uex

    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
       uex = (u(np) + u(np+1)) / 2
       Fd  = Fd  + dfm * uex
    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, hmedio
    integer :: npe

    ! 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    = po * 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 + 1 ) / ( ( 1 + rcag ) ** 2 )
    KL2 = ( 8 * gama - 27 ) / ( 2304 * ( 1 + rcag ) )
    KL3 = 754 * ( gama ** 2 ) - ( 757 * gama ) + 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')

    open(12,file=trim(adjustl(caso))//'.Richardson_3p0')

    i   = 1
    j   = ny-1
    np  = (j-1)*nx + i
    hmedio = ( x(nx-1) - x(1) ) / ( nx - 2 )

    write(12,5) hmedio,  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
 
    5 format(  2x, 1pe25.15, ' = h médio em X (m)',/, &
             ' 1', 1pe25.15, ' = coeficiente de descarga',/, &
             ' 2', 1pe25.15, ' = empuxo dinâmico',/, &
             ' 3', 1pe25.15, ' = empuxo de pressão ao nível do mar (p = 101325 Pa)',/, &
             ' 4', 1pe25.15, ' = empuxo de pressão no vácuo',/, &
             ' 5', 1pe25.15, ' = empuxo total ao nível do mar',/, &
             ' 6', 1pe25.15, ' = empuxo total no vácuo',/, &
             ' 7', 1pe25.15, ' = coeficiente de empuxo ao nível do mar',/, &
             ' 8', 1pe25.15, ' = coeficiente de empuxo no vácuo',/, &
             ' 9', 1pe25.15, ' = velocidade característica',/, &
             '10', 1pe25.15, ' = velocidade de ejeção efetiva ao nível do mar',/, &
             '11', 1pe25.15, ' = velocidade de ejeção efetiva no vácuo',/, &
             '12', 1pe25.15, ' = impulso específico ao nível do mar',/, &
             '13', 1pe25.15, ' = impulso específico no vácuo')

    close(12)

  end subroutine resultados_principais

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

  subroutine curva_residuo

    integer :: k ! auxiliar

    nome = 'Mach2D_5p5_residuo.gnu'

!    open(7,file=trim(adjustl(nome)))
    open(7,file=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_5p5_residuo.gnu'

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

  end subroutine curva_residuo

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

end module resultados
