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,
!            escrever resultados,
!            apresentar gráficos e
!            listar resultados.

! Versão atual: 2 Ago 09

! Relação das sub-rotinas em ordem de apresentação:
!    1) solucao_numerica
!    2) pos_processamento
!    3) print1
!    7) grafico_x_temperatura
!   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 Parede2D 1.1, de 2 Ago 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 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,'Tmax2D', t38,'q2D' )
    write(10,2) 
    2 format(//,t3,'Iteração',t14,'L1(R)n / L1(R)1', /)
    write(11,3) 
    3 format('# Iteração',t14,'L1(R)n / L1(R)1')

    ! zeragem do tempo de CPU
    call cpu_time ( tcpu1 )

    ! ciclo iterativo externo no tempo ----------------------------------------

    do it = ito, itmax

       ! cálculo de propriedades
       call condutividade_termica_nos
       call condutividade_termica_faces

       ! solução de T
       call coeficientes_fonte_temperatura
       call ccT
       call lu2d9(a,dl9,du9,nxy,nx,ny)
       call fb2d9(a,dl9,du9,ny,nx,nxy,T,toli,nitm,b)
       call norma_L1_9d (nx,ny,nxy,a,b,T,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
          write(11,4) it, norma/norma1
          4 format(2x,i8,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 )

    tcpu = tcpu2-tcpu1 + tcpuo

    open(99,file=trim(adjustl(caso))//'.bin',form='unformatted')
    write(99) it,norma1,tcpu,T
    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 :: 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 e armazena valores dos contornos nos fictícios
    call contornos_pos_processamento

    ! escreve campos

    if ( w_cam == 1 ) then

       write(10,09)
       09 format(//,10x,'*** CAMPOS ***', //, &
                    10x,'Os nós fictícios contêm os valores dos contornos.', //, &
                    10x,'Na horizontal (esquerda para direita): contorno oeste para leste.', //, &
                    10x,'Na vertical   (de cima para baixo):    contorno norte para sul.' )

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

       write(10,57)
       57 format(//,30x,'kp: condutividade térmica nodal (W/m.K)')
       call print1(kp,ny,nx,nxy)

       write(10,58)
       58 format(//,30x,'ke: condutividade térmica na face leste (W/m.K)')
       call print1(ke,ny,nx,nxy)

       write(10,59)
       59 format(//,30x,'kn: condutividade térmica na face norte (W/m.K)')
       call print1(kn,ny,nx,nxy)

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

       write(10,61)
       61 format(//,30x,'XP: coordenada X no centro de cada volume (m)')
       call print1(xp,ny,nx,nxy)

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

       write(10,62)
       62 format(//,30x,'YP: coordenada Y no centro de cada volume (m)')
       call print1(yp,ny,nx,nxy)

       write(10,71)
       71 format(//,30x,'R: raio no canto nordeste (m)')
       call print1(raio,ny,nx,nxy)

       write(10,72)
       72 format(//,30x,'RP: raio no centro de cada volume (m)')
       call print1(rp,ny,nx,nxy)

       write(10,73)
       73 format(//,30x,'re: raio na face leste (m)')
       call print1(re,ny,nx,nxy)

       write(10,74)
       74 format(//,30x,'rn: raio na face norte (m)')
       call print1(rn,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,23)
       23 format(//,30x,'Jp: Jacobiano no centro do volume (1/m2)')
       call print1(Jp,ny,nx,nxy)

       write(10,33)
       33 format(//,30x,'Je: Jacobiano na face leste (1/m2)')
       call print1(Je,ny,nx,nxy)

       write(10,34)
       34 format(//,30x,'Jn: Jacobiano na face norte (1/m2)')
       call print1(Jn,ny,nx,nxy)

       write(10,81)
       81 format(//,30x,'Alfae: Alfa na face leste (m2)')
       call print1(Alfae,ny,nx,nxy)

       write(10,82)
       82 format(//,30x,'Betae: Beta na face leste (m2)')
       call print1(Betae,ny,nx,nxy)

       write(10,83)
       83 format(//,30x,'Gaman: Gama na face norte (m2)')
       call print1(Gaman,ny,nx,nxy)

       write(10,84)
       84 format(//,30x,'Betan: Beta na face norte (m2)')
       call print1(Betan,ny,nx,nxy)

    end if

    ! 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')

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

    call grafico_x_temperatura

    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')

    ! 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')

  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 :: 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_temperatura

    integer :: nps, npn

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

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

    write(7,1)
    1 format('### TEMPERATURAS NAS PAREDES (K) ###', /, &
             '#', t4,'x (m)', t28,'wh-1D', t52,'wh-2D', t76,'wc-1D', t100,'wc-2D')

    write(10,2)
    2 format(/,'### TEMPERATURAS NAS PAREDES (K) ###', /, &
             t4,'x (m)', t28,'wh-1D', t52,'wh-2D', t76,'wc-1D', t100,'wc-2D')

    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 ) then
          write( 7,3) xp(nps), Twh1D(i), T(nps), Twc1D(i), T(npn)
          write(10,3) xp(nps), Twh1D(i), T(nps), Twc1D(i), T(npn)
          3 format(5(1pe24.15))
       end if
    end do

    close(7)

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

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

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

    write(7,5) trim(adjustl(nome))
    5 format("replot '", a, "' using 1:3 title 'wh-2D'")

    write(7,6) trim(adjustl(nome))
    6 format("replot '", a, "' using 1:4 title 'wc-1D'")

    write(7,7) trim(adjustl(nome))
    7 format("replot '", a, "' using 1:5 title 'wc-2D'")

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

    close(7)

    nome = 'wgnuplot Parede2D_1p1_temperatura.gnu'

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

  end subroutine grafico_x_temperatura

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

  subroutine globais

    real*8 :: dfm, uex

    write(8,1) it, Tmax2D, q2D

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

  end subroutine globais

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

  subroutine resultados_principais

    write(10,1) minval(Twc1D), maxval(Twh1D), maxval(Twh1D)-minval(Twc1D), q1D
    1 format(/,3x,'MODELO UNIDIMENSIONAL (1D)', &
             /,5x,'Twcmin: temperatura mínima da parede fria   (K) = ',1pe24.15, &
             /,5x,'Twhmax: temperatura máxima da parede quente (K) = ',1pe24.15, &
             /,5x,'Twhmax - Twcmin: diferença máxima de temperatura (K) = ',1pe24.15, &
             /,5x,'q1D: taxa de transferência de calor (W) = ',1pe24.15)

    write(10,2) Tmin2D, Tmax2D, Tmax2D-Tmin2D, q2D
    2 format(/,3x,'MODELO BIDIMENSIONAL (2D)', &
             /,5x,'Twcmin: temperatura mínima da parede fria   (K) = ',1pe24.15, &
             /,5x,'Twhmax: temperatura máxima da parede quente (K) = ',1pe24.15, &
             /,5x,'Twhmax - Twcmin: diferença máxima de temperatura (K) = ',1pe24.15, &
             /,5x,'q2D: taxa de transferência de calor (W) = ',1pe24.15)

    write(10,3) Tmax2D/maxval(Twh1D), q2D/q1D
    3 format(//,5x,'razão Tmax-2D/1D = ',1pe24.15, &
             //,5x,'razão q2D / q1D  = ',1pe24.15)

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

    write(12,4) hmedio, Tmax2D/maxval(Twh1D), q2D/q1D
    4 format(  2x, 1pe25.15, ' = h médio (m)',/, &
             ' 1', 1pe25.15, ' = razão Tmax-2D/1D',/, &
             ' 2', 1pe25.15, ' = razão q2D / q1D')

    close(12)

  end subroutine resultados_principais

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

  subroutine curva_residuo

    nome = 'Parede2D_1p1_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 Parede2D_1p1_residuo.gnu'

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

  end subroutine curva_residuo

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

end module resultados
