module resultados

! objetivo: escrever resultados e apresentar gráficos

use calculos

implicit none

contains

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

  subroutine mostra_resultados

    call escreve_ordens_h
    call escreve_ordens_i
    call escreve_ordens_bi
    call escreve_ordens_c12
    call escreve_ordens_c13

    call escreve_Th_e_extrapolacoes
    call escreve_dTh_e_deltas_extrapolacoes
    if ( com_solucao_analitica == 1 ) &
       call escreve_erros_Th_e_extrapolacoes

    call escreve_incertezas
    if ( com_solucao_analitica == 1 ) &
       call escreve_Us_por_Eh
    call escreve_GCI_e_outros
    if ( com_solucao_analitica == 1 ) &
       call escreve_razoes_Es
    call escreve_razao_entre_incertezas

    call escreve_multiextrapolacoes
    call escreve_deltas_multiextrapolacoes
    if ( com_solucao_analitica == 1 ) &
       call escreve_erros_multiextrapolacoes
    call escreve_incertezas_multiextrapolacoes
    if ( com_solucao_analitica == 1 ) then
       call escreve_razao_incertezas_erros_multiextrapolacoes
       call escreve_Eh_por_Em
    end if
    call escreve_multiextrapolacoes_Uri_por_Um
	
    call grafico_ordens

    if ( com_solucao_analitica == 1 ) then
       call escreve_erros_para_graficos
       call grafico_erros
    end if
    call grafico_incertezas_1
    call grafico_incertezas_2

  end subroutine mostra_resultados

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

  subroutine escreve_ordens_h

    write(out,1) trim(adjustl(pVt(1))), trim(adjustl(pVt(1)))
    1 format (/, t3,'h', t15,'pE(Th)', t37,'pU(Th)', t59,'pE(Th) - pL(',a,')', &
                t81,'pU(Th) - pL(',a,')',  t104,'it_pU(Th)' )

    write(out,2) h(1)

    if ( pE_h(2) < limite ) then
       write(out,2) h(2), pE_h(2), pE_h(2)-pL
    else
       write(out,2) h(2)
    end if

    2 format(1pe12.4,1pe22.14,22x,1pe22.14)
    3 format(1pe12.4,4(1pe22.14),i10)
    4 format(1pe12.4,22x,1pe22.14,22x,1pe22.14,i10)

    if ( nm >= 3 ) then
       do i = 3, nm
          if ( pE_h(i) < limite .and. pU_h(i) < limite ) then
             write(out,3) h(i), pE_h(i), pU_h(i), pE_h(i)-pL, pU_h(i)-pL, it_pU_h(i)
          end if
          if ( pE_h(i) < limite .and. pU_h(i) > limite ) then
             write(out,2) h(i), pE_h(i), pE_h(i)-pL
          end if
          if ( pE_h(i) > limite .and. pU_h(i) < limite ) then
             write(out,4) h(i), pU_h(i), pU_h(i)-pL, it_pU_h(i)
          end if
          if ( pE_h(i) > limite .and. pU_h(i) > limite ) then
             write(out,2) h(i)
          end if
       end do
    end if

    open(out2,file='pL.dat')
    write(out2,7)
    7 format ('#', t4,'log(h)', t27,'pL' )
    do i = 1, nm
       write(out2,8) dlog10(h(i)), pL
    end do
    8 format(1pe12.4,I8)
    close(out2)

    open(out2,file='pE_h.dat')
    write(out2,5)
    5 format ('#', t4,'log(h)', t27,'pE_h' )
    do i = 2, nm
       if ( pE_h(i) < limite ) then
          write(out2,2) dlog10(h(i)), pE_h(i)
       else
          write(out2,2)
       end if
    end do
    close(out2)

    open(out2,file='pU_h.dat')
    write(out2,6)
    6 format ('#', t4,'log(h)', t27,'pU_h' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( pU_h(i) < limite ) then
             write(out2,2) dlog10(h(i)), pU_h(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

  end subroutine escreve_ordens_h

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

  subroutine escreve_ordens_i

    write(out,1) trim(adjustl(pVt(2))), trim(adjustl(pVt(2)))
    1 format (/, t3,'h', t15,'pE(Ti_pU)', t37,'pU(Ti_pU)', t59,'pE(Ti_pU) - p2(',a,')', &
                t81,'pU(Ti_pU) - p2(',a,')',  t104,'it_pU(Ti_pU)' )

    do i = 1, 3
       if ( i > nm ) then
          exit
       else
          write(out,2) h(i)
       end if
    end do

    if ( nm >= 4 ) then
       i = 4
       if ( pE_i(i) < limite ) then
          write(out,2) h(i), pE_i(i), pE_i(i)-pV(2)
       else
          write(out,2) h(i)
       end if
    end if

    2 format(1pe12.4,1pe22.14,22x,1pe22.14)
    3 format(1pe12.4,4(1pe22.14),i10)
    4 format(1pe12.4,22x,1pe22.14,22x,1pe22.14,i10)

    if ( nm >= 5 ) then
       do i = 5, nm
          if ( pE_i(i) < limite .and. pU_i(i) < limite ) then
             write(out,3) h(i), pE_i(i), pU_i(i), pE_i(i)-pV(2), pU_i(i)-pV(2), it_pU_i(i)
          end if
          if ( pE_i(i) < limite .and. pU_i(i) > limite ) then
             write(out,2) h(i), pE_i(i), pE_i(i)-pV(2)
          end if
          if ( pE_i(i) > limite .and. pU_i(i) < limite ) then
             write(out,4) h(i), pU_i(i), pU_i(i)-pV(2), it_pU_i(i)
          end if
          if ( pE_i(i) > limite .and. pU_i(i) > limite ) then
             write(out,2) h(i)
          end if
       end do
    end if

    open(out2,file='p2.dat')
    write(out2,7)
    7 format ('#', t4,'log(h)', t27,'pV(2)' )
    do i = 1, nm
       write(out2,8) dlog10(h(i)), pV(2)
    end do
    8 format(1pe12.4,I8)
    close(out2)

    open(out2,file='pE_i.dat')
    write(out2,5)
    5 format ('#', t4,'log(h)', t27,'pE_i' )
    if ( nm >= 4 ) then
       do i = 4, nm
          if ( pE_i(i) < limite ) then
             write(out2,2) dlog10(h(i)), pE_i(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='pU_i.dat')
    write(out2,6)
    6 format ('#', t4,'log(h)', t27,'pU_i' )
    if ( nm >= 5 ) then
       do i = 5, nm
          if ( pU_i(i) < limite ) then
             write(out2,2) dlog10(h(i)), pU_i(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

  end subroutine escreve_ordens_i

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

  subroutine escreve_ordens_bi

    write(out,1) trim(adjustl(pVt(3))), trim(adjustl(pVt(3)))
    1 format (/, t3,'h', t15,'pE(Tbi_pU)', t37,'pU(Tbi_pu)', t59,'pE(Tbi_pU) - p3(',a,')', &
                t81,'pU(Tbi_pU) - p3(',a,')',  t104,'it_pU(Tbi_pU)' )

    do i = 1, 5
       if ( i > nm ) then
          exit
       else
          write(out,2) h(i)
       end if
    end do

    if ( nm >= 6 ) then
       i = 6
       if ( pE_bi(i) < limite ) then
          write(out,2) h(i), pE_bi(i), pE_bi(i)-pV(3)
       else
          write(out,2) h(i)
       end if
    end if

    2 format(1pe12.4,1pe22.14,22x,1pe22.14)
    3 format(1pe12.4,4(1pe22.14),i10)
    4 format(1pe12.4,22x,1pe22.14,22x,1pe22.14,i10)

    if ( nm >= 7 ) then
       do i = 7, nm
          if ( pE_bi(i) < limite .and. pU_bi(i) < limite ) then
             write(out,3) h(i), pE_bi(i), pU_bi(i), pE_bi(i)-pV(3), pU_bi(i)-pV(3), it_pU_bi(i)
          end if
          if ( pE_bi(i) < limite .and. pU_bi(i) > limite ) then
             write(out,2) h(i), pE_bi(i), pE_bi(i)-pV(3)
          end if
          if ( pE_bi(i) > limite .and. pU_bi(i) < limite ) then
             write(out,4) h(i), pU_bi(i), pU_bi(i)-pV(3), it_pU_bi(i)
          end if
          if ( pE_bi(i) > limite .and. pU_bi(i) > limite ) then
             write(out,2) h(i)
          end if
       end do
    end if

    open(out2,file='p3.dat')
    write(out2,7)
    7 format ('#', t4,'log(h)', t27,'pV(3)' )
    do i = 1, nm
       write(out2,8) dlog10(h(i)), pV(3)
    end do
    8 format(1pe12.4,I8)
    close(out2)

    open(out2,file='pE_bi.dat')
    write(out2,5)
    5 format ('#', t4,'log(h)', t27,'pE_bi' )
    if ( nm >= 6 ) then
       do i = 6, nm
          if ( pE_bi(i) < limite ) then
             write(out2,2) dlog10(h(i)), pE_bi(i)
          else
             write(out2,2)
          end if
       end do
       close(out2)
    end if
    close(out2)

    open(out2,file='pU_bi.dat')
    write(out2,6)
    6 format ('#', t4,'log(h)', t27,'pU_bi' )
    if ( nm >= 7 ) then
       do i = 7, nm
          if ( pU_bi(i) < limite ) then
             write(out2,2) dlog10(h(i)), pU_bi(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

  end subroutine escreve_ordens_bi

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

  subroutine escreve_ordens_c12

    write(out,1) trim(adjustl(pVt(2))), trim(adjustl(pVt(2)))
    1 format (/, t3,'h', t15,'pE(Tc_12)', t37,'pU(Tc_12)', t59,'pE(Tc_12) - p2(',a,')', &
                t81,'pU(Tc_12) - p2(',a,')',  t104,'it_pU(Tc_12)' )

    do i = 1, 3
       if ( i > nm ) then
          exit
       else
          write(out,2) h(i)
       end if
    end do

    if ( nm >= 4 ) then
       i = 4
       if ( pE_c12(i) < limite ) then
          write(out,2) h(i), pE_c12(i), pE_c12(i)-pV(2)
       else
          write(out,2) h(i)
       end if
    end if

    2 format(1pe12.4,1pe22.14,22x,1pe22.14)
    3 format(1pe12.4,4(1pe22.14),i10)
    4 format(1pe12.4,22x,1pe22.14,22x,1pe22.14,i10)

    if ( nm >= 5 ) then
       do i = 5, nm
          if ( pE_c12(i) < limite .and. pU_c12(i) < limite ) then
             write(out,3) h(i), pE_c12(i), pU_c12(i), pE_c12(i)-pV(2), pU_c12(i)-pV(2), it_pU_c12(i)
          end if
          if ( pE_c12(i) < limite .and. pU_c12(i) > limite ) then
             write(out,2) h(i), pE_c12(i), pE_c12(i)-pV(2)
          end if
          if ( pE_c12(i) > limite .and. pU_c12(i) < limite ) then
             write(out,4) h(i), pU_c12(i), pU_c12(i)-pV(2), it_pU_c12(i)
          end if
          if ( pE_c12(i) > limite .and. pU_c12(i) > limite ) then
             write(out,2) h(i)
          end if
       end do
    end if

    open(out2,file='pE_c12.dat')
    write(out2,5)
    5 format ('#', t4,'log(h)', t27,'pE_c12' )
    if ( nm >= 4 ) then
       do i = 4, nm
          if ( pE_c12(i) < limite ) then
             write(out2,2) dlog10(h(i)), pE_c12(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='pU_c12.dat')
    write(out2,6)
    6 format ('#', t4,'log(h)', t27,'pU_c12' )
    if ( nm >= 5 ) then
       do i = 5, nm
          if ( pU_c12(i) < limite ) then
             write(out2,2) dlog10(h(i)), pU_c12(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

  end subroutine escreve_ordens_c12

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

  subroutine escreve_ordens_c13

    write(out,1) trim(adjustl(pVt(2))), trim(adjustl(pVt(2)))
    1 format (/, t3,'h', t15,'pE(Tc_13)', t37,'pU(Tc_13)', t59,'pE(Tc_13) - p2(',a,')', &
                t81,'pU(Tc_13) - p2(',a,')',  t104,'it_pU(Tc_13)' )

    do i = 1, 3
       if ( i > nm ) then
          exit
       else
          write(out,2) h(i)
       end if
    end do

    if ( nm >= 4 ) then
       i = 4
       if ( pE_c13(i) < limite ) then
          write(out,2) h(i), pE_c13(i), pE_c13(i)-pV(2)
       else
          write(out,2) h(i)
       end if
    end if

    2 format(1pe12.4,1pe22.14,22x,1pe22.14)
    3 format(1pe12.4,4(1pe22.14),i10)
    4 format(1pe12.4,22x,1pe22.14,22x,1pe22.14,i10)

    if ( nm >= 5 ) then
       do i = 5, nm
          if ( pE_c13(i) < limite .and. pU_c13(i) < limite ) then
             write(out,3) h(i), pE_c13(i), pU_c13(i), pE_c13(i)-pV(2), pU_c13(i)-pV(2), it_pU_c13(i)
          end if
          if ( pE_c13(i) < limite .and. pU_c13(i) > limite ) then
             write(out,2) h(i), pE_c13(i), pE_c13(i)-pV(2)
          end if
          if ( pE_c13(i) > limite .and. pU_c13(i) < limite ) then
             write(out,4) h(i), pU_c13(i), pU_c13(i)-pV(2), it_pU_c13(i)
          end if
          if ( pE_c13(i) > limite .and. pU_c13(i) > limite ) then
             write(out,2) h(i)
          end if
       end do
    end if

    open(out2,file='pE_c13.dat')
    write(out2,5)
    5 format ('#', t4,'log(h)', t27,'pE_c13' )
    if ( nm >= 4 ) then
       do i = 4, nm
          if ( pE_c13(i) < limite ) then
             write(out2,2) dlog10(h(i)), pE_c13(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='pU_c13.dat')
    write(out2,6)
    6 format ('#', t4,'log(h)', t27,'pU_c13' )
    if ( nm >= 5 ) then
       do i = 5, nm
          if ( pU_c13(i) < limite ) then
             write(out2,2) dlog10(h(i)), pU_c13(i)
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

  end subroutine escreve_ordens_c13

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

  subroutine escreve_Th_e_extrapolacoes

    write(out,1)
    1 format (/,  t3,'h',      t15,'Th',     t37,'Ti_12', t59,'Ti_13',  &
                 t81,'Ti_pU', t103,'Ti_23', t125,'Ti_bi', 17x,'Ti_tri', &
                t169,'Tc_12', t191,'Tc_13', t213,'Tbi_pU' )

    write(out,2) h(1), Th(1)

    write(out,2) h(2), Th(2), Ti_12(2)

    2 format(1pe12.4,10(1pe22.14))
    3 format(1pe12.4,3(1pe22.14),22x,3(1pe22.14),44x,1pe22.14)
    4 format(1pe12.4,6(1pe22.14),22x,3(1pe22.14))

    if ( nm >= 3 ) then
       i = 3
          if ( Ti_pU(i) < limite ) &
             write(out,4) h(i), Th(i),    Ti_12(i),  Ti_13(i), Ti_pU(i), Ti_23(i), &
                                Ti_bi(i),            Tc_12(i), Tc_13(i)
          if ( Ti_pU(i) > limite ) &
             write(out,3) h(i), Th(i),    Ti_12(i),  Ti_13(i),           Ti_23(i), &
                                Ti_bi(i)
    end if

    if ( nm >= 4 ) then
       i = 4
          if ( Ti_pU(i) < limite ) &
             write(out,2) h(i), Th(i),    Ti_12(i),  Ti_13(i), Ti_pU(i), Ti_23(i), &
                                Ti_bi(i), Ti_tri(i), Tc_12(i), Tc_13(i)
          if ( Ti_pU(i) > limite ) &
             write(out,3) h(i), Th(i),    Ti_12(i),  Ti_13(i),           Ti_23(i), &
                                Ti_bi(i), Ti_tri(i)
    end if

    if ( nm >= 5 ) then
       do i = 5, nm
          if ( Ti_pU(i) < limite .and. Tbi_pU(i) < limite ) &
             write(out,2) h(i), Th(i),    Ti_12(i),  Ti_13(i), Ti_pU(i), Ti_23(i), &
                                Ti_bi(i), Ti_tri(i), Tc_12(i), Tc_13(i), Tbi_pU(i)
             if ( Ti_pU(i) < limite .and. Tbi_pU(i) > limite) &
             write(out,2) h(i), Th(i),    Ti_12(i),  Ti_13(i), Ti_pU(i), Ti_23(i), &
                                Ti_bi(i), Ti_tri(i), Tc_12(i), Tc_13(i)
          if ( Ti_pU(i) > limite .and. Tbi_pU(i) < limite) &
             write(out,3) h(i), Th(i),    Ti_12(i),  Ti_13(i),           Ti_23(i), &
                                Ti_bi(i), Ti_tri(i),                     Tbi_pU(i)
          if ( Ti_pU(i) > limite .and. Tbi_pU(i) > limite ) &
             write(out,3) h(i), Th(i),    Ti_12(i),  Ti_13(i),           Ti_23(i), &
                                Ti_bi(i), Ti_tri(i)                               
       end do
    end if

  end subroutine escreve_Th_e_extrapolacoes

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

  subroutine escreve_dTh_e_deltas_extrapolacoes

    real*8 :: dTh, dTi_12, dTi_13, dTi_pU, dTi_23, dTi_bi, dTi_tri, &
                   dTc_12, dTc_13, dTbi_pU

    write(out,1)
    1 format (/,  t3,'h',       t15,'dTh',     t37,'dTi_12', t59,'dTi_13',  &
                 t81,'dTi_pU', t103,'dTi_23', t125,'dTi_bi', 16x,'dTi_tri', &
                t169,'dTc_12', t191,'dTc_13', t213,'dTbi_pU' )

    write(out,2) h(1)

    write(out,2) h(2), Th(2)-Th(1)

    if ( nm >= 3 ) write(out,2) h(3), Th(3)-Th(2), Ti_12(3)-Ti_12(2)

    2 format(1pe12.4,10(1pe22.14))
    3 format(1pe12.4,3(1pe22.14),22x,3(1pe22.14),44x,1pe22.14)
    4 format(1pe12.4,6(1pe22.14),22x,3(1pe22.14))

    if ( nm >= 4 ) then
       i = 4
          dTh = Th(i) - Th(i-1)
          dTi_12 = Ti_12(i) - Ti_12(i-1)
          dTi_13 = Ti_13(i) - Ti_13(i-1)
          dTi_pU = Ti_pU(i) - Ti_pU(i-1)
          dTi_23 = Ti_23(i) - Ti_23(i-1)
          dTi_bi = Ti_bi(i) - Ti_bi(i-1)
          dTc_12 = Tc_12(i) - Tc_12(i-1)
          dTc_13 = Tc_13(i) - Tc_13(i-1)
          if ( Ti_pU(i) < limite ) &
             write(out,4) h(i), dTh,    dTi_12,  dTi_13, dTi_pU, dTi_23, &
                                dTi_bi,          dTc_12, dTc_13
          if ( Ti_pU(i) > limite ) &
             write(out,3) h(i), dTh,    dTi_12,  dTi_13,         dTi_23, &
                                dTi_bi
    end if

    if ( nm >= 5 ) then
       i = 5
          dTh = Th(i) - Th(i-1)
          dTi_12 = Ti_12(i) - Ti_12(i-1)
          dTi_13 = Ti_13(i) - Ti_13(i-1)
          dTi_pU = Ti_pU(i) - Ti_pU(i-1)
          dTi_23 = Ti_23(i) - Ti_23(i-1)
          dTi_bi = Ti_bi(i) - Ti_bi(i-1)
          dTi_tri = Ti_tri(i) - Ti_tri(i-1)
          dTc_12 = Tc_12(i) - Tc_12(i-1)
          dTc_13 = Tc_13(i) - Tc_13(i-1)
          if ( Ti_pU(i) < limite ) &
             write(out,2) h(i), dTh,    dTi_12,  dTi_13, dTi_pU, dTi_23, &
                                dTi_bi, dTi_tri, dTc_12, dTc_13
          if ( Ti_pU(i) > limite ) &
             write(out,3) h(i), dTh,    dTi_12,  dTi_13,         dTi_23, &
                                dTi_bi, dTi_tri
    end if

    if ( nm >= 6 ) then
       do i = 6, nm
          dTh = Th(i) - Th(i-1)
          dTi_12 = Ti_12(i) - Ti_12(i-1)
          dTi_13 = Ti_13(i) - Ti_13(i-1)
          dTi_pU = Ti_pU(i) - Ti_pU(i-1)
          dTi_23 = Ti_23(i) - Ti_23(i-1)
          dTi_bi = Ti_bi(i) - Ti_bi(i-1)
          dTi_tri = Ti_tri(i) - Ti_tri(i-1)
          dTc_12 = Tc_12(i) - Tc_12(i-1)
          dTc_13 = Tc_13(i) - Tc_13(i-1)
          dTbi_pU = Tbi_pU(i) - Tbi_pU(i-1)
          if ( Ti_pU(i) < limite .and. Tbi_pU(i) < limite ) &
             write(out,2) h(i), dTh,    dTi_12,  dTi_13, dTi_pU, dTi_23, &
                                dTi_bi, dTi_tri, dTc_12, dTc_13, dTbi_pU
             if ( Ti_pU(i) < limite .and. Tbi_pU(i) > limite) &
             write(out,2) h(i), dTh,    dTi_12,  dTi_13, dTi_pU, dTi_23, &
                                dTi_bi, dTi_tri, dTc_12, dTc_13
          if ( Ti_pU(i) > limite .and. Tbi_pU(i) < limite) &
             write(out,3) h(i), dTh,    dTi_12,  dTi_13,         dTi_23, &
                                dTi_bi, dTi_tri,                 dTbi_pU
          if ( Ti_pU(i) > limite .and. Tbi_pU(i) > limite ) &
             write(out,3) h(i), dTh,    dTi_12,  dTi_13,         dTi_23, &
                                dTi_bi, dTi_tri                              
       end do
    end if

  end subroutine escreve_dTh_e_deltas_extrapolacoes

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

  subroutine escreve_erros_Th_e_extrapolacoes

    write(out,1)
    1 format (/,  t3,'h',         t15,'E(Th)',     t37,'E(Ti_12)', t59,'E(Ti_13)',  &
                 t81,'E(Ti_pU)', t103,'E(Ti_23)', t125,'E(Ti_bi)', 14x,'E(Ti_tri)', &
                t169,'E(Tc_12)', t191,'E(Tc_13)', t213,'E(Tbi_pU)' )

    write(out,2) h(1), Eh(1)

    write(out,2) h(2), Eh(2), Ei_12(2)

    2 format(1pe12.4,10(1pe22.14))
    3 format(1pe12.4,3(1pe22.14),22x,3(1pe22.14),44x,1pe22.14)
    4 format(1pe12.4,6(1pe22.14),22x,3(1pe22.14))

    if ( nm >= 3 ) then
       i = 3
          if ( Ti_pU(i) < limite ) &
             write(out,4) h(i), Eh(i),    Ei_12(i),  Ei_13(i), Ei_pU(i), Ei_23(i), &
                                Ei_bi(i),            Ec_12(i), Ec_13(i)
          if ( Ti_pU(i) > limite ) &
             write(out,3) h(i), Eh(i),    Ei_12(i),  Ei_13(i),           Ei_23(i), &
                                Ei_bi(i)
    end if

    if ( nm >= 4 ) then
       i = 4
          if ( Ti_pU(i) < limite ) &
             write(out,2) h(i), Eh(i),    Ei_12(i),  Ei_13(i), Ei_pU(i), Ei_23(i), &
                                Ei_bi(i), Ei_tri(i), Ec_12(i), Ec_13(i)
          if ( Ti_pU(i) > limite ) &
             write(out,3) h(i), Eh(i),    Ei_12(i),  Ei_13(i),           Ei_23(i), &
                                Ei_bi(i), Ei_tri(i)
    end if

    if ( nm >= 5 ) then
       do i = 5, nm
          if ( Ti_pU(i) < limite .and. Tbi_pU(i) < limite ) &
             write(out,2) h(i), Eh(i),    Ei_12(i),  Ei_13(i), Ei_pU(i), Ei_23(i), &
                                Ei_bi(i), Ei_tri(i), Ec_12(i), Ec_13(i), Ebi_pU(i)
             if ( Ti_pU(i) < limite .and. Tbi_pU(i) > limite) &
             write(out,2) h(i), Eh(i),    Ei_12(i),  Ei_13(i), Ei_pU(i), Ei_23(i), &
                                Ei_bi(i), Ei_tri(i), Ec_12(i), Ec_13(i)
          if ( Ti_pU(i) > limite .and. Tbi_pU(i) < limite) &
             write(out,3) h(i), Eh(i),    Ei_12(i),  Ei_13(i),           Ei_23(i), &
                                Ei_bi(i), Ei_tri(i),                     Ebi_pU(i)
          if ( Ti_pU(i) > limite .and. Tbi_pU(i) > limite ) &
             write(out,3) h(i), Eh(i),    Ei_12(i),  Ei_13(i),           Ei_23(i), &
                                Ei_bi(i), Ei_tri(i)                               
       end do
    end if

  end subroutine escreve_erros_Th_e_extrapolacoes

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

  subroutine escreve_incertezas

    real*8 :: aux1, aux2, aux3, aux4, aux5 ! variáveis auxiliares

    if ( com_solucao_analitica == 1 ) then

       write(out,1)
       1 format (/,  t3,'h',           t15,'E(Th)',       t37,'Uri_12(Th)', t59,'Uri_13(Th)',  &
                    t81,'Uri_pU(Th)', t103,'Uri_23(Th)', t125,'Uri_bi(Th)', 12x,'Uri_tri(Th)', &
                   t169,'Uc(Tc_12)',  t191,'Uc(Tc_13)',  t213,'E(Tc_12)',  t235,'E(Tc_13)'     )

       write(out,2) h(1), Eh(1)

       write(out,2) h(2), Eh(2), Uri_12(2)

       2 format(1pe12.4,11(1pe22.14))
       3 format(1pe12.4,3(1pe22.14),22x,3(1pe22.14),44x,1pe22.14)
       4 format(1pe12.4,6(1pe22.14),22x,4(1pe22.14))

       if ( nm >= 3 ) then
          i = 3
          if ( Ti_pU(i) < limite ) &
             write(out,4) h(i), Eh(i),     Uri_12(i),  Uri_13(i), Uri_pU(i), Uri_23(i), &
                                Uri_bi(i),             Uc_12(i),  Uc_13(i),  Ec_12(i),  &
                                Ec_13(i)
          if ( Ti_pU(i) > limite ) &
             write(out,3) h(i), Eh(i),     Uri_12(i),  Uri_13(i),            Uri_23(i), &
                                Uri_bi(i)
       end if

       if ( nm >= 4 ) then
          do i = 4, nm
             if ( Ti_pU(i) < limite ) &
                write(out,2) h(i), Eh(i),     Uri_12(i),  Uri_13(i), Uri_pU(i), Uri_23(i), &
                                   Uri_bi(i), Uri_tri(i), Uc_12(i),  Uc_13(i),  Ec_12(i),  &
                                   Ec_13(i)
             if ( Ti_pU(i) > limite ) &
                write(out,3) h(i), Eh(i),     Uri_12(i),  Uri_13(i),            Uri_23(i), &
                                   Uri_bi(i), Uri_tri(i)
          end do
       end if

    else

       write(out,11)
       11 format (/,  t3,'h',           t15,'Uri_12(Th)',  t37,'Uri_13(Th)',  t59,'Uri_pU(Th)', t81,'Uri_23(Th)', &
                    t103,'Uri_bi(Th)', t125,'Uri_tri(Th)', 11x,'Uc(Tc_12)',  t169,'Uc(Tc_13)' )

       write(out,12) h(1)
 
       write(out,12) h(2), Uri_12(2)

       12 format(1pe12.4,8(1pe22.14))
       13 format(1pe12.4,2(1pe22.14),22x,3(1pe22.14),44x,1pe22.14)
       14 format(1pe12.4,5(1pe22.14),22x,2(1pe22.14))

       if ( nm >= 3 ) then
          i = 3
          if ( Ti_pU(i) < limite ) &
             write(out,14) h(i), Uri_12(i),  Uri_13(i), Uri_pU(i), Uri_23(i), &
                                 Uri_bi(i),             Uc_12(i),  Uc_13(i)
          if ( Ti_pU(i) > limite ) &
             write(out,13) h(i), Uri_12(i),  Uri_13(i),            Uri_23(i), &
                                 Uri_bi(i)
       end if

       if ( nm >= 4 ) then
          do i = 4, nm
             if ( Ti_pU(i) < limite ) &
                write(out,12) h(i), Uri_12(i), Uri_13(i),  Uri_pU(i), Uri_23(i), &
                                    Uri_bi(i), Uri_tri(i), Uc_12(i),  Uc_13(i)
             if ( Ti_pU(i) > limite ) &
                write(out,13) h(i), Uri_12(i), Uri_13(i),             Uri_23(i), &
                                    Uri_bi(i), Uri_tri(i)
          end do
       end if

    end if

    open(out2,file='U_12.dat')
    write(out2,7)
    7 format ('#', t3,'log(h)', t15,'log|Ud|', t37,'log|Uri_12|', t59,'log(GCI_pL)' )
    do i = 2, nm
       if ( Ud(i) /= 0.0d0 ) then
          aux1 = dlog10(dabs(Ud(i)))
          aux2 = dlog10(dabs(Uri_12(i)))
          aux3 = dlog10(dabs(GCI_pL(i)))
          write(out2,2) dlog10(h(i)), aux1, aux2, aux3
       else
          write(out2,2)
       end if
    end do
    close(out2)

    open(out2,file='Uri_13.dat')
    write(out2,31)
    31 format ('#', t3,'log(h)', t15,'log|Uri_13|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Uri_13(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Uri_13(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Uri_23.dat')
    write(out2,32)
    32 format ('#', t3,'log(h)', t15,'log|Uri_23|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Uri_23(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Uri_23(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Uri_bi.dat')
    write(out2,33)
    33 format ('#', t3,'log(h)', t15,'log|Uri_bi|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Uri_bi(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Uri_bi(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Uri_pU.dat')
    write(out2,6)
    6 format ('#', t3,'log(h)', t15,'log|Uri_pU|',  t37,'log|Uri_p|', t59,'log(GCI_p)' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Uri_pU(i) /= 0.0d0 .and. Uri_pU(i) < limite ) then
             aux1 = dlog10(dabs(Uri_pU(i)))
             aux2 = dlog10(dabs(Uri_p(i)))
             aux3 = dlog10(dabs(GCI_p(i)))
             write(out2,2) dlog10(h(i)), aux1, aux2, aux3
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Uc_12.dat')
    write(out2,16)
    16 format ('#', t3,'log(h)', t15,'log(Uc_12)' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Uc_12(i) /= 0.0d0 .and. Uc_12(i) < limite ) then
             aux1 = dlog10(dabs(Uc_12(i)))
             write(out2,2) dlog10(h(i)), aux1
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Uc_13.dat')
    write(out2,17)
    17 format ('#', t3,'log(h)', t15,'log(Uc_13)' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Uc_13(i) /= 0.0d0 .and. Uc_13(i) < limite ) then
             aux1 = dlog10(dabs(Uc_13(i)))
             write(out2,2) dlog10(h(i)), aux1
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Utri.dat')
    write(out2,8)
    8 format ('#', t3,'log(h)', t15,'log|Uri_tri|' )
    if ( nm >= 4 ) then
       do i = 4, nm
          if ( Uri_tri(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Uri_tri(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

  end subroutine escreve_incertezas

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

  subroutine escreve_Us_por_Eh

    write(out,1)
    1 format (/, t3,'h',           t15,'Ud/|Eh|-1', t27,'U12/Eh-1', t39,'U13/Eh-1',  &
                t51,'UpU/Eh-1',    t63,'U23/Eh-1',  t75,'Ubi/Eh-1', t87,'Utri/Eh-1', &
                t99,'U/|Ec12|-1', t111,'U/|Ec13|-1' )

    write(out,2) h(1)

    i = 2
    if ( Eh(i) /= 0.0d0 ) write(out,2) h(i), Ud(i)/dabs(Eh(i))-1, Uri_12(i)/Eh(i)-1

    2 format(10(1pe12.4))
    3 format( 4(1pe12.4), 12x, 3(1pe12.4) )
    4 format( 7(1pe12.4), 12x, 2(1pe12.4) )

    if ( nm >= 3 ) then
       i = 3
       if ( Ti_pU(i) < limite .and. Eh(i) /= 0.0d0 .and. Ec_12(i) /= 0.0d0 .and. Ec_13(i) /= 0.0d0 ) &
          write(out,4) h(i), Ud(i)/dabs(Eh(i))-1, Uri_12(i)/Eh(i)-1, Uri_13(i)/Eh(i)-1, &
                             Uri_pU(i)/Eh(i)-1,   Uri_23(i)/Eh(i)-1, Uri_bi(i)/Eh(i)-1, &
                                                  Uc_12(i)/dabs(Ec_12(i))-1, &
                             Uc_13(i)/dabs(Ec_13(i))-1
       if ( Ti_pU(i) > limite .and. Eh(i) /= 0.0d0 ) &
          write(out,3) h(i), Ud(i)/dabs(Eh(i))-1, Uri_12(i)/Eh(i)-1, Uri_13(i)/Eh(i)-1, &
                                                  Uri_23(i)/Eh(i)-1, Uri_bi(i)/Eh(i)-1
    end if

    if ( nm >= 4 ) then
       do i = 4, nm
          if ( Ti_pU(i) < limite .and. Eh(i) /= 0.0d0 .and. Ec_12(i) /= 0.0d0 .and. Ec_13(i) /= 0.0d0 ) &
             write(out,2) h(i), Ud(i)/dabs(Eh(i))-1, Uri_12(i)/Eh(i)-1, Uri_13(i)/Eh(i)-1, &
                                Uri_pU(i)/Eh(i)-1,   Uri_23(i)/Eh(i)-1, Uri_bi(i)/Eh(i)-1, &
                                Uri_tri(i)/Eh(i)-1,  Uc_12(i)/dabs(Ec_12(i))-1, &
                                Uc_13(i)/dabs(Ec_13(i))-1
          if ( Ti_pU(i) > limite .and. Eh(i) /= 0.0d0 ) &
             write(out,3) h(i), Ud(i)/dabs(Eh(i))-1, Uri_12(i)/Eh(i)-1, Uri_13(i)/Eh(i)-1, &
                                                     Uri_23(i)/Eh(i)-1, Uri_bi(i)/Eh(i)-1, &
                                Uri_tri(i)/Eh(i)-1
       end do
    end if

  end subroutine escreve_Us_por_Eh

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

  subroutine escreve_GCI_e_outros

    if ( com_solucao_analitica == 1 ) then

       write(out,1)
       1 format (/, t3,'h',           &
                   t15,'Uri_p(Th)',   t37,'GCI_pL(Th)',    t59,'GCI_p(Th)',   &
                   t81,'Urip/U12-1',  t93,'GpL/|U12|-Fs', t105,'Gp/|U12|-Fs', &
                  t117,'Urip/Eh-1',  t129,'GpL/|Eh|-Fs',  t141,'Gp/|Eh|-Fs'   )

       write(out,2) h(1)

       i = 2
       if ( Eh(i) /= 0.0d0 .and. Uri_12(i) /= 0.0d0 ) &
          write(out,3) h(i),                       GCI_pL(i),                                   &
                                                   GCI_pL(i)/dabs(Uri_12(i))-Fs,                       &
                                                   GCI_pL(i)/dabs(Eh(i))-Fs

       if ( nm >= 3 ) then
          do i = 3, nm
             if ( Ti_pU(i) < limite .and. Eh(i) /= 0.0d0 .and. Uri_12(i) /= 0.0d0 ) &
                write(out,2) h(i), Uri_p(i),             GCI_pL(i),                    GCI_p(i),                    &
                                   Uri_p(i)/Uri_12(i)-1, GCI_pL(i)/dabs(Uri_12(i))-Fs, GCI_p(i)/dabs(Uri_12(i))-Fs, &
                                   Uri_p(i)/Eh(i)-1,     GCI_pL(i)/dabs(Eh(i))-Fs,     GCI_p(i)/dabs(Eh(i))-Fs
             if ( Ti_pU(i) > limite .and. Eh(i) /= 0.0d0 .and. Uri_12(i) /= 0.0d0 ) &
                write(out,3) h(i),                       GCI_pL(i),                                          &
                                                         GCI_pL(i)/dabs(Uri_12(i))-Fs,                       &
                                                         GCI_pL(i)/dabs(Eh(i))-Fs
          end do
       end if

    else

       write(out,6)
       6 format (/, t3,'h',           &
                   t15,'Uri_p(Th)',   t37,'GCI_pL(Th)',    t59,'GCI_p(Th)',   &
                   t81,'Urip/U12-1',  t93,'GpL/|U12|-Fs', t105,'Gp/|U12|-Fs'  )

       write(out,2) h(1)

       i = 2
       if ( Uri_12(i) /= 0.0d0 ) &
          write(out,3) h(i),                       GCI_pL(i),                                   &
                                                   GCI_pL(i)/dabs(Uri_12(i))-Fs

       if ( nm >= 3 ) then
          do i = 3, nm
             if ( Ti_pU(i) < limite .and. Uri_12(i) /= 0.0d0 ) &
                write(out,2) h(i), Uri_p(i),             GCI_pL(i),                    GCI_p(i),                    &
                                   Uri_p(i)/Uri_12(i)-1, GCI_pL(i)/dabs(Uri_12(i))-Fs, GCI_p(i)/dabs(Uri_12(i))-Fs
             if ( Ti_pU(i) > limite .and. Uri_12(i) /= 0.0d0 ) &
                write(out,3) h(i),                       GCI_pL(i),                                          &
                                                         GCI_pL(i)/dabs(Uri_12(i))-Fs
          end do
       end if

    end if

    2 format( 1pe12.4, 3(1pe22.14), 6(1pe12.4) )
    3 format( 1pe12.4, 22x, 1pe22.14, 34x, 1pe12.4, 24x, 1pe12.4 )

  end subroutine escreve_GCI_e_outros

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

  subroutine escreve_razoes_Es

    write(out,1)
    1 format (/, t3,'h',          t15,'Eh/Ei12',    t27,'Ei13/Ei12',   t39,'EipU/Ei12', &
                t51,'Ei23/Ei12',  t63,'E_bi/Ei12',  t75,'E_tri/Ei12',  t87,'Ec12/Ei12', &
                t99,'Ec13/Ei12', t111,'Eh/Ec12',   t123,'Eh/Ec13',    t135,'Ec13/Ec12'  )

    write(out,2) h(1)

    i = 2
    if ( Ei_12(i) /= 0.0d0 ) write(out,2) h(i), Eh(i)/Ei_12(i)

    2 format(12(1pe12.4))
    3 format( 3(1pe12.4), 12x, 3(1pe12.4) )
    4 format( 6(1pe12.4), 12x, 5(1pe12.4) )

    if ( nm >= 3 ) then
       i = 3
          if ( Ti_pU(i) < limite .and. Ei_12(i) /= 0.0d0 .and. Ec_12(i) /= 0.0d0 .and. Ec_13(i) /= 0.0d0 ) &
             write(out,4) h(i),              Eh(i)/Ei_12(i),    Ei_13(i)/Ei_12(i), Ei_pU(i)/Ei_12(i), &
                          Ei_23(i)/Ei_12(i), Ei_bi(i)/Ei_12(i),                    Ec_12(i)/Ei_12(i), &
                          Ec_13(i)/Ei_12(i), Eh(i)/Ec_12(i),    Eh(i)/Ec_13(i),    Ec_13(i)/Ec_12(i) 
          if ( Ti_pU(i) > limite .and. Ei_12(i) /= 0.0d0 ) &
             write(out,3) h(i),                    Eh(i)/Ei_12(i),       Ei_13(i)/Ei_12(i),    &
                          Ei_23(i)/Ei_12(i),       Ei_bi(i)/Ei_12(i)
    end if

    if ( nm >= 4 ) then
       do i = 4, nm
          if ( Ti_pU(i) < limite .and. Ei_12(i) /= 0.0d0 .and. Ec_12(i) /= 0.0d0 .and. Ec_13(i) /= 0.0d0 ) &
             write(out,2) h(i),              Eh(i)/Ei_12(i),    Ei_13(i)/Ei_12(i),  Ei_pU(i)/Ei_12(i), &
                          Ei_23(i)/Ei_12(i), Ei_bi(i)/Ei_12(i), Ei_tri(i)/Ei_12(i), Ec_12(i)/Ei_12(i), &
                          Ec_13(i)/Ei_12(i), Eh(i)/Ec_12(i),    Eh(i)/Ec_13(i),     Ec_13(i)/Ec_12(i) 
          if ( Ti_pU(i) > limite .and. Ei_12(i) /= 0.0d0 ) &
             write(out,3) h(i),              Eh(i)/Ei_12(i),    Ei_13(i)/Ei_12(i),  &
                          Ei_23(i)/Ei_12(i), Ei_bi(i)/Ei_12(i), Ei_tri(i)/Ei_12(i)
       end do
    end if

  end subroutine escreve_razoes_Es

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

  subroutine escreve_razao_entre_incertezas

    write(out,1)
    1 format (/, t3,'h',           t15,'Ud/U12-1',  t27,'U13/U12-1',  t39,'UpU/U12-1',  &
                t51,'U23/U12-1',   t63,'Ubi/U12-1', t75,'Utri/U12-1', t87,'|U12|/Uc12', &
                t99,'|U12|/Uc13', t111,'Uc13/Uc12'  )

    write(out,2) h(1)

    i = 2
    if ( Uri_12(i) /= 0.0d0 ) write(out,2) h(i), Ud(i)/Uri_12(i)-1

    2 format(10(1pe12.4))
    3 format( 3(1pe12.4), 12x, 3(1pe12.4) )
    4 format( 6(1pe12.4), 12x, 3(1pe12.4) )

    if ( nm >= 3 ) then
       i = 3
       if ( Ti_pU(i) < limite .and. Uri_12(i) /= 0.0d0 .and. Uc_12(i) /= 0.0d0 .and. Uc_13(i) /= 0.0d0 ) &
          write(out,4) h(i), Ud(i)/Uri_12(i)-1,        Uri_13(i)/Uri_12(i)-1,    Uri_pU(i)/Uri_12(i)-1,  &
                             Uri_23(i)/Uri_12(i)-1,    Uri_bi(i)/Uri_12(i)-1,                            &
                             dabs(Uri_12(i))/Uc_12(i), dabs(Uri_12(i))/Uc_13(i), Uc_13(i)/Uc_12(i)
       if ( Ti_pU(i) > limite .and. Uri_12(i) /= 0.0d0 ) &
          write(out,3) h(i), Ud(i)/Uri_12(i)-1,        Uri_13(i)/Uri_12(i)-1,                            &
                             Uri_23(i)/Uri_12(i)-1,    Uri_bi(i)/Uri_12(i)-1
    end if

    if ( nm >= 4 ) then
       do i = 4, nm
          if ( Ti_pU(i) < limite .and. Uri_12(i) /= 0.0d0 .and. Uc_12(i) /= 0.0d0 .and. Uc_13(i) /= 0.0d0 ) &
             write(out,2) h(i), Ud(i)/Uri_12(i)-1,        Uri_13(i)/Uri_12(i)-1,    Uri_pU(i)/Uri_12(i)-1,  &
                                Uri_23(i)/Uri_12(i)-1,    Uri_bi(i)/Uri_12(i)-1,    Uri_tri(i)/Uri_12(i)-1, &
                                dabs(Uri_12(i))/Uc_12(i), dabs(Uri_12(i))/Uc_13(i), Uc_13(i)/Uc_12(i)
          if ( Ti_pU(i) > limite .and. Uri_12(i) /= 0.0d0 ) &
             write(out,3) h(i), Ud(i)/Uri_12(i)-1,        Uri_13(i)/Uri_12(i)-1,                            &
                                Uri_23(i)/Uri_12(i)-1,    Uri_bi(i)/Uri_12(i)-1,    Uri_tri(i)/Uri_12(i)-1
       end do
    end if

  end subroutine escreve_razao_entre_incertezas

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

  subroutine escreve_multiextrapolacoes

    character*2   :: num
    character*500 :: texto

    texto = 'h           Th(p=' // trim(adjustl(pVt(1))) // ')'
    do m = 1, nm-1
       write(num,2) m
       2 format(i2)
       if ( m == 1 ) &
          texto = trim(adjustl(texto)) // '  ' // '             T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // ')'
       if ( m >= 2 ) &
          texto = trim(adjustl(texto)) // '             T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // ')'
    end do
    write(out,1) trim(adjustl(texto))
    1 format (/, 2x, a )

    write(out,3) h(1), Th(1)

    3 format(1pe12.4,30(1pe22.14))

    do i = 2, nm
       write(out,3) h(i), Th(i), ( Tm_12(i,m), m=1,i-1 )
    end do

  end subroutine escreve_multiextrapolacoes

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

  subroutine escreve_deltas_multiextrapolacoes

    character*2   :: num
    character*500 :: texto

    texto = 'h           dTh(p=' // trim(adjustl(pVt(1))) // ')'
    do m = 1, nm-1
       write(num,2) m
       2 format(i2)
       if ( m == 1 ) &
          texto = trim(adjustl(texto)) // '  ' // '            dT_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // ')'
       if ( m >= 2 ) &
          texto = trim(adjustl(texto)) // '            dT_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // ')'
    end do
    write(out,1) trim(adjustl(texto))
    1 format (/, 2x, a )

    write(out,3) h(1)

    write(out,3) h(2), Th(2)-Th(1)

    3 format(1pe12.4,30(1pe22.14))

    if ( nm >= 3 ) then
       do i = 3, nm
          write(out,3) h(i), Th(i)-Th(i-1), ( Tm_12(i,m)-Tm_12(i-1,m), m=1,i-2 )
       end do
    end if

  end subroutine escreve_deltas_multiextrapolacoes

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

  subroutine escreve_erros_multiextrapolacoes

    character*2   :: num
    character*500 :: texto

    texto = 'h           E(Th(p=' // trim(adjustl(pVt(1))) // '))'
    do m = 1, nm-1
       write(num,2) m
       2 format(i2)
       if ( m == 1 ) &
          texto = trim(adjustl(texto)) // '  ' // '          E(T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // '))'
       if ( m >= 2 ) &
          texto = trim(adjustl(texto)) // '          E(T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // '))'
    end do
    write(out,1) trim(adjustl(texto))
    1 format (/, 2x, a )

    write(out,3) h(1), Eh(1)

    3 format(1pe12.4,30(1pe22.14))

    do i = 2, nm
       write(out,3) h(i), Eh(i), ( T_exato-Tm_12(i,m), m=1,i-1 )
    end do

  end subroutine escreve_erros_multiextrapolacoes

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

  subroutine escreve_Eh_por_Em

    real*8        :: del
    character*2   :: num
    character*500 :: texto

    texto = 'h           E(Th)/E(Th)'
    do m = 1, nm-1
       write(num,2) m
       2 format(i2)
       if ( m == 1 ) &
          texto = trim(adjustl(texto)) // '  ' // '         E(Th)/E(T_' // &
                  trim(adjustl(num)) // 'i)'
       if ( m >= 2 ) &
          texto = trim(adjustl(texto)) // '         E(Th)/E(T_' // &
                  trim(adjustl(num)) // 'i)'
    end do
    write(out,1) trim(adjustl(texto))
    1 format (/, 2x, a )

    if ( Eh(1) /= 0.0d0 ) write(out,3) h(1), Eh(1)/Eh(1)

    3 format(1pe12.4,30(1pe22.14))

    do i = 2, nm
       del = 1.0d0
       do m = 1, i-1
          if ( T_exato-Tm_12(i,m) == 0.0d0 ) del = 0.0d0
       end do
       if ( Eh(i) /= 0.0d0 .and. del /= 0.0d0 ) &
          write(out,3) h(i), Eh(i)/Eh(i), ( Eh(i)/(T_exato-Tm_12(i,m)), m=1,i-1 )
    end do

  end subroutine escreve_Eh_por_Em

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

  subroutine escreve_incertezas_multiextrapolacoes

    character*2   :: num
    character*500 :: texto

    texto = 'h           Uri_12(Th(p=' // trim(adjustl(pVt(1))) // '))'
    do m = 1, nm-1
       write(num,2) m
       2 format(i2)
       if ( m == 1 ) &
          texto = trim(adjustl(texto)) // '  ' // '     Uri_12(T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // '))'
       if ( m >= 2 ) &
          texto = trim(adjustl(texto)) // '     Uri_12(T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // '))'
    end do
    write(out,1) trim(adjustl(texto))
    1 format (/, 2x, a )

    write(out,3) h(1)

    write(out,3) h(2), Uri_12(2)

    3 format(1pe12.4,30(1pe22.14))

    if ( nm >= 3 ) then
       do i = 3, nm
          write(out,3) h(i), Uri_12(i), ( Um_12(i,m), m=1,i-2 )
       end do
    end if

  end subroutine escreve_incertezas_multiextrapolacoes

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

  subroutine escreve_razao_incertezas_erros_multiextrapolacoes

    real*8        :: del
    character*2   :: num
    character*500 :: texto

    texto = 'h           (U/E-1)(Th(p=' // trim(adjustl(pVt(1))) // '))'
    do m = 1, nm-1
       write(num,2) m
       2 format(i2)
       if ( m == 1 ) &
          texto = trim(adjustl(texto)) // '  ' // '    (U/E-1)(T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // '))'
       if ( m >= 2 ) &
          texto = trim(adjustl(texto)) // '    (U/E-1)(T_' // &
                  trim(adjustl(num)) // 'i(p=' // trim(adjustl(pVt(m+1))) // '))'
    end do
    write(out,1) trim(adjustl(texto))
    1 format (/, 2x, a )

    write(out,3) h(1)

    if ( Eh(2) /= 0.0d0 ) write(out,3) h(2), Uri_12(2)/Eh(2)-1

    3 format(1pe12.4,30(1pe22.14))

    if ( nm >= 3 ) then
       do i = 3, nm
          del = 1.0d0
          do m = 1, i-2
             if ( T_exato-Tm_12(i,m) == 0.0d0 ) del = 0.0d0
          end do
          if ( Eh(i) /= 0.0d0 .and. del /= 0.0d0 ) &
             write(out,3) h(i), Uri_12(i)/Eh(i)-1, ( Um_12(i,m)/(T_exato-Tm_12(i,m))-1, m=1,i-2 )
       end do
    end if

  end subroutine escreve_razao_incertezas_erros_multiextrapolacoes

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

  subroutine escreve_multiextrapolacoes_Uri_por_Um

    real*8        :: del
    character*2   :: num
    character*500 :: texto

    texto = 'h           U12(Th)/U12(Th)'
    do m = 1, nm-1
       write(num,2) m
       2 format(i2)
       if ( m == 1 ) &
          texto = trim(adjustl(texto)) // '  ' // '     U12(Th)/U12(T_' // &
                  trim(adjustl(num)) // 'i)'
       if ( m >= 2 ) &
          texto = trim(adjustl(texto)) // '     U12(Th)/U12(T_' // &
                  trim(adjustl(num)) // 'i)'
    end do
    write(out,1) trim(adjustl(texto))
    1 format (/, 2x, a )

    write(out,3) h(1)

    if ( Uri_12(2) /= 0.0d0 ) write(out,3) h(2), Uri_12(2)/Uri_12(2)

    3 format(1pe12.4,30(1pe22.14))

    if ( nm >= 3 ) then
       do i = 3, nm
          del = 1.0d0
          do m = 1, i-2
             if ( Um_12(i,m) == 0.0d0 ) del = 0.0d0
          end do
          if ( Uri_12(i) /= 0.0d0 .and. del /= 0.0d0 ) &
             write(out,3) h(i), Uri_12(i)/Uri_12(i), ( Uri_12(i)/Um_12(i,m), m=1,i-2 )
       end do
    end if

  end subroutine escreve_multiextrapolacoes_Uri_por_Um

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

  subroutine grafico_ordens

    if ( com_solucao_analitica == 1 ) then

       ! gráfico de pE e pU

       open(out2,file='Richardson_3p0_pE_pU.gnu')

       do k = 1, 20
          read(out2,*)
       end do

       write(out2,1) trim(adjustl(head))
       1 format("set title '", a, /, "replot")

       close(out2)

       dos = system('wgnuplot Richardson_3p0_pE_pU.gnu')

    else

       ! gráfico de pU

       open(out2,file='Richardson_3p0_pU.gnu')

       do k = 1, 15
          read(out2,*)
       end do

       write(out2,1) trim(adjustl(head))

       close(out2)

       dos = system('wgnuplot Richardson_3p0_pU.gnu')

    end if

  end subroutine grafico_ordens

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

  subroutine escreve_erros_para_graficos

    real*8 :: aux1, aux2, aux3, aux4, Emin ! variáveis auxiliares

    open(out2,file='Eh.dat')
    write(out2,6)
    6 format ('#', t4,'log(h)', t27,'log|Eh|' )
    do i = 1, nm
       if ( Eh(i) /= 0.0d0 ) then
          write(out2,2) dlog10(h(i)), dlog10(dabs(Eh(i)))
       else
          write(out2,2)
       end if
    end do
    close(out2)

    open(out2,file='Ei_12.dat')
    write(out2,1)
    1 format ('#', t4,'log(h)', t27,'log|Ei_12|' )
    do i = 2, nm
       if ( Ei_12(i) /= 0.0d0 ) then
          write(out2,2) dlog10(h(i)), dlog10(dabs(Ei_12(i)))
       else
          write(out2,2)
       end if
    end do
    2 format(6(1pe23.14))
    close(out2)

    open(out2,file='Ei_13.dat')
    write(out2,3)
    3 format ('#', t4,'log(h)', t27,'log|Ei_13|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Ei_13(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Ei_13(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Ei_23.dat')
    write(out2,9)
    9 format ('#', t4,'log(h)', t27,'log|Ei_23|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Ei_23(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Ei_23(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Ei_bi.dat')
    write(out2,19)
    19 format ('#', t4,'log(h)', t27,'log|Ei_bi|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Ei_bi(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Ei_bi(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Ei_tri.dat')
    write(out2,18)
    18 format ('#', t4,'log(h)', t27,'log|Ei_tri|' )
    if ( nm >= 4 ) then
       do i = 4, nm
          if ( Ei_tri(i) /= 0.0d0 ) then
             write(out2,2) dlog10(h(i)), dlog10(dabs(Ei_tri(i)))
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Ei_pU.dat')
    write(out2,4)
    4 format ('#', t4,'log(h)', t27,'log|Ei_pU|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if ( Ei_pU(i) /= 0.0d0 .and. Ti_pU(i) < limite ) then
             aux1 = dlog10(dabs(Ei_pU(i)))
             write(out2,2) dlog10(h(i)), aux1
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Ec_12.dat')
    write(out2,14)
    14 format ('#', t4,'log(h)', t27,'log|Ec_12|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if (  Ec_12(i) /= 0.0d0 .and. Tc_12(i) < limite ) then
             aux1 = dlog10(dabs(Ec_12(i)))
             write(out2,2) dlog10(h(i)), aux1
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Ec_13.dat')
    write(out2,15)
    15 format ('#', t4,'log(h)', t27,'log|Ec_13|' )
    if ( nm >= 3 ) then
       do i = 3, nm
          if (  Ec_13(i) /= 0.0d0 .and. Tc_13(i) < limite ) then
             aux1 = dlog10(dabs(Ec_13(i)))
             write(out2,2) dlog10(h(i)), aux1
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='Ebi.dat')
    write(out2,5)
    5 format ('#', t4,'log(h)', t27,'log|Ebi_pU|' )
    if ( nm >= 5 ) then
       do i = 5, nm
          if ( Ebi_pU(i) /= 0.0d0 .and. Ebi_pU(i) < limite ) then
             aux1 = dlog10(dabs(Ebi_pU(i)))
             write(out2,2) dlog10(h(i)), aux1
          else
             write(out2,2)
          end if
       end do
    end if
    close(out2)

    open(out2,file='E_pL_p2.dat')
    write(out2,7)
    7 format ('#', t4,'log(h)', t27,'log|E_pL|', t50,'log|E_p2|')
    if ( Eh(1) /= 0.0d0 ) then
       aux1 = dlog10(dabs(Eh(1)))
       write(out2,2) dlog10(h(1)), aux1, aux1
       aux1 = dlog10(dabs(Eh(1)*((h(nm)/h(1))**pL)))
       aux2 = dlog10(dabs(Eh(1)*((h(nm)/h(1))**pV(2))))
       write(out2,2) dlog10(h(nm)), aux1, aux2
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1, aux1
       aux1 = dlog10(1.0d0*((h(nm)/h(1))**pL))
       aux2 = dlog10(1.0d0*((h(nm)/h(1))**pV(2)))
       write(out2,2) dlog10(h(nm)), aux1, aux2
    end if
    close(out2)

    Emin = dmin1( dabs(Eh(1)),     dabs(Ei_12(nm)), dabs(Ei_13(nm)),  dabs(Ei_23(nm)), &
                  dabs(Ei_pU(nm)), dabs(Ei_bi(nm)), dabs(Ei_tri(nm)), dabs(Ec_12(nm)), &
                  dabs(Ec_13(nm)), dabs(Ebi_pU(nm)) )

    open(out2,file='E_p3.dat')
    write(out2,8)
    8 format ('#', t4,'log(h)', t27,'log|E_p3|')
    if ( Eh(1) /= 0.0d0 .and. Emin /= 0.0d0 ) then
       aux1 = dlog10(dabs(Eh(1)))
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(dabs(Emin))
       aux2 = (dabs(Emin/Eh(1))*(h(1)**pV(3)))**(1.0d0/pV(3)) ! h de Emin
       write(out2,2) dlog10(aux2), aux1
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(1.0d-3)
       aux2 = (1.0d-3*(h(1)**pV(3)))**(1.0d0/pV(3)) ! h de Emin
       write(out2,2) dlog10(aux2), aux1
    end if
    close(out2)

  end subroutine escreve_erros_para_graficos

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

  subroutine grafico_erros

    open(out2,file='Richardson_3p0_E.gnu')

    do k = 1, 19
       read(out2,*)
    end do

    write(out2,1) trim(adjustl(head))
    1 format("set title '", a, /, "replot")

    close(out2)

    dos = system('wgnuplot Richardson_3p0_E.gnu')

  end subroutine grafico_erros

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

  subroutine grafico_incertezas_1

    real*8 :: aux1, aux2, aux3, Umin ! variáveis auxiliares

    open(out2,file='U_pL.dat')
    write(out2,9)
    9 format ('#', t3,'log(h)', t15,'log|E_pL|')
    if ( Uri_12(2) /= 0.0d0 ) then
       aux1 = dlog10(dabs(Uri_12(2)))
       write(out2,2) dlog10(h(2)), aux1
       aux1 = dlog10(dabs(Uri_12(2)*((h(nm)/h(2))**pL)))
       write(out2,2) dlog10(h(nm)), aux1
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(1.0d0*((h(nm)/h(1))**pL))
       write(out2,2) dlog10(h(nm)), aux1
    end if
    close(out2)

    Umin = dmin1 ( dabs(Uri_12(2)),  dabs(Uri_12(nm)), dabs(Ud(nm)),     dabs(Uri_13(nm)), &
                   dabs(Uri_23(nm)), dabs(Uri_pU(nm)), dabs(Uri_bi(nm)), dabs(Uri_tri(nm)) )

    open(out2,file='U_p2.dat')
    write(out2,10)
    10 format ('#', t3,'log(h)', t15,'log|E_p2|')
    if ( Uri_12(2) /= 0.0d0 .and. Umin /= 0.0d0 ) then
       aux1 = dlog10(dabs(Uri_12(2)))
       write(out2,2) dlog10(h(2)), aux1
       aux1 = dlog10(dabs(Umin))
       aux2 = (dabs(Umin/Uri_12(2))*(h(2)**pV(2)))**(1.0d0/pV(2)) ! h de Umin para p2
       write(out2,2) dlog10(aux2), aux1
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(1.0d-3)
       aux2 = (1.0d-3*(h(1)**pV(2)))**(1.0d0/pV(2)) ! h de Umin
       write(out2,2) dlog10(aux2), aux1
    end if
    close(out2)

    open(out2,file='U_p3.dat')
    write(out2,20)
    20 format ('#', t3,'log(h)', t15,'log|E_p3|')
    if ( Uri_12(2) /= 0.0d0 .and. Umin /= 0.0d0 ) then
       aux1 = dlog10(dabs(Uri_12(2)))
       write(out2,2) dlog10(h(2)), aux1
       aux1 = dlog10(dabs(Umin))
       aux3 = (dabs(Umin/Uri_12(2))*(h(2)**pV(3)))**(1.0d0/pV(3)) ! h de Umin para p3
       write(out2,2) dlog10(aux3), aux1
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(1.0d-3)
       aux2 = (1.0d-3*(h(1)**pV(3)))**(1.0d0/pV(3)) ! h de Umin
       write(out2,2) dlog10(aux2), aux1
    end if
    close(out2)

    2 format(1pe12.4,11(1pe22.14))

    if ( com_solucao_analitica == 1 ) then

       ! gráfico de erro e incertezas

       open(out2,file='Richardson_3p0_E_U.gnu')

       do k = 1, 19
          read(out2,*)
       end do

       write(out2,1) trim(adjustl(head))
       1 format("set title '", a, /, "replot")

       close(out2)

       dos = system('wgnuplot Richardson_3p0_E_U.gnu')

    else

       ! gráfico de incertezas

       open(out2,file='Richardson_3p0_U.gnu')

       do k = 1, 16
          read(out2,*)
       end do

       write(out2,1) trim(adjustl(head))

       close(out2)

       dos = system('wgnuplot Richardson_3p0_U.gnu')

    end if

  end subroutine grafico_incertezas_1

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

  subroutine grafico_incertezas_2

    real*8 :: aux1, aux2, aux3, Umin ! variáveis auxiliares

    open(out2,file='U_pL.dat')
    write(out2,9)
    9 format ('#', t3,'log(h)', t15,'log|E_pL|')
    if ( Uri_12(2) /= 0.0d0 ) then
       aux1 = dlog10(dabs(Uri_12(2)))
       write(out2,2) dlog10(h(2)), aux1
       aux1 = dlog10(dabs(Uri_12(2)*((h(nm)/h(2))**pL)))
       write(out2,2) dlog10(h(nm)), aux1
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(1.0d0*((h(nm)/h(1))**pL))
       write(out2,2) dlog10(h(nm)), aux1
    end if
    close(out2)

    Umin = dmin1 ( dabs(Uri_12(2)),  dabs(Uri_12(nm)), dabs(Uri_13(nm)), &
                   dabs(Uri_pU(nm)), dabs(Uc_12(nm)),  dabs(Uc_13(nm))   )

    open(out2,file='U_p2.dat')
    write(out2,10)
    10 format ('#', t3,'log(h)', t15,'log|E_p2|')
    if ( Uri_12(2) /= 0.0d0 .and. Umin /= 0.0d0 ) then
       aux1 = dlog10(dabs(Uri_12(2)))
       write(out2,2) dlog10(h(2)), aux1
       aux1 = dlog10(dabs(Umin))
       aux2 = (dabs(Umin/Uri_12(2))*(h(2)**pV(2)))**(1.0d0/pV(2)) ! h de Umin para p2
       write(out2,2) dlog10(aux2), aux1
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(1.0d-3)
       aux2 = (1.0d-3*(h(1)**pV(2)))**(1.0d0/pV(2)) ! h de Umin
       write(out2,2) dlog10(aux2), aux1
    end if
    close(out2)

    open(out2,file='U_p3.dat')
    write(out2,20)
    20 format ('#', t3,'log(h)', t15,'log|E_p3|')
    if ( Uri_12(2) /= 0.0d0 .and. Umin /= 0.0d0 ) then
       aux1 = dlog10(dabs(Uri_12(2)))
       write(out2,2) dlog10(h(2)), aux1
       aux1 = dlog10(dabs(Umin))
       aux3 = (dabs(Umin/Uri_12(2))*(h(2)**pV(3)))**(1.0d0/pV(3)) ! h de Umin para p3
       write(out2,2) dlog10(aux3), aux1
    else
       aux1 = dlog10(1.0d0)
       write(out2,2) dlog10(h(1)), aux1
       aux1 = dlog10(1.0d-3)
       aux2 = (1.0d-3*(h(1)**pV(3)))**(1.0d0/pV(3)) ! h de Umin
       write(out2,2) dlog10(aux2), aux1
    end if
    close(out2)

    2 format(1pe12.4,11(1pe22.14))

    if ( com_solucao_analitica == 1 ) then

       ! gráfico de erro e incertezas

       open(out2,file='Richardson_3p0_E_U_2.gnu')

       do k = 1, 20
          read(out2,*)
       end do

       write(out2,1) trim(adjustl(head))
       1 format("set title '", a, /, "replot")

       close(out2)

       dos = system('wgnuplot Richardson_3p0_E_U_2.gnu')

    else

       ! gráfico de incertezas

       open(out2,file='Richardson_3p0_U_2.gnu')

       do k = 1, 15
          read(out2,*)
       end do

       write(out2,1) trim(adjustl(head))

       close(out2)

       dos = system('wgnuplot Richardson_3p0_U_2.gnu')

    end if

  end subroutine grafico_incertezas_2

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

end module resultados