! UNIVERSIDADE FEDERAL DO PARANÁ
! DEPARTAMENTO DE ENGENHARIA MECÂNICA
! PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA (PG-MEC)
! LAPLACE 2D EM MALHAS TRIANGULARES	IMPARES(linear)
! PROFESSOR ORIENTADOR: CARLOS HENRIQUE MARCHI
! PROFESSOR CO-ORIENTADOR: LUCIANO KIYOSHI ARAKI
! ALUNA: ARILEIDE CRISTINA ALVES


program Laplace_Malha2D_triangular

use portlib
use msfwin
use msflib

!********************************************************************
! Declaração de variáveis *******************************************
!********************************************************************

implicit none

integer :: Np	    ! Número de pontos
integer :: Ne		! Número de elementos
integer :: Ner      ! Número de elementos reais
integer :: Nef      ! Número de elementos fictícios
integer :: qer      ! Quantidade de elementos reais contiguos ao contorno leste
integer :: passoe   ! Passo para numero do elemento real contiguo ao contorno leste
integer :: passon   ! Passo para numero do elemento real contiguo ao contorno norte
integer :: passow   ! Passo para numero do elemento real contiguo ao contorno oeste
integer :: passos   ! Passo para numero do elemento real contiguo ao contorno sul
integer :: e        ! Número do elemento
integer :: f        ! Número do elemento
integer :: iti      ! Contador	de iterações do Gauss Seidel
integer :: cont     ! Contador	para temperatura no meio do dominio (saídas Richardson)
integer :: numiti   ! Numero de iteraçoes no Gauss Seidel
integer :: p        ! Número do ponto
integer :: aux_fict	! Auxiliar para contagem dos volumes fictícios
integer :: malha    ! Nível da malha
integer :: drf      ! Número de elementos reais da malha/ número de elementos fictícios
integer :: necce  	! Número do elemento contíguo do real ao leste na malha 1
integer :: neccs  	! Número do elemento contíguo do real ao sul na malha 1
integer :: neccn  	! Número do elemento contíguo do real ao norte na malha 1
integer :: neccw  	! Número do elemento contíguo do real ao oeste na malha 1

real*8  :: T1, T2, TP       !Auxiliares na contagem do tempo de cpu
real*8  :: areatotal		! Area do dominio de calculo
real*8  :: h				! Métrica da malha
real*8  :: somanorm			! Auxiliar no calculo da norma do erro
real*8  :: somatemp			! Auxiliar no calculo da temperatura media
real*8  :: Tmedia			! Variavel de interesse: Temperatura Média
real*8  :: somatte          ! Auxiliar no calculo da taxa de transferencia de calor (contorno leste)
real*8  :: somattn          ! Auxiliar no calculo da taxa de transferencia de calor (contorno norte)
real*8  :: somattw          ! Auxiliar no calculo da taxa de transferencia de calor (contorno oeste)
real*8  :: somatts			! Auxiliar no calculo da taxa de transferencia de calor (contorno sul)
real*8  :: tte  			! Variavel de interesse: Taxa de Transferencia de Calor (contorno leste)
real*8  :: ttn  			! Variavel de interesse: Taxa de Transferencia de Calor (contorno norte)
real*8  :: ttw  			! Variavel de interesse: Taxa de Transferencia de Calor (contorno oeste)
real*8  :: tts  			! Variavel de interesse: Taxa de Transferencia de Calor (contorno sul)
real*8  :: T    			! Variavel de interesse: Temperatura no meio do dominio (Richardson)
real*8  :: area				! Area de cada elemento triangular
real*8  :: mx, my			! Pontos médios das arestas
real*8  :: norx, nory		! Normais às arestas
real*8  :: dbaric			! Distância entre baricentros
real*8  :: vunitx, vunity	! Vetores unitários
real*8  :: Jx, Jy			! Vetores auxiliares no calculo do fluxo de temperatura 
real*8  :: soma_fict		! Auxiliar para Gauss-Seidel (fictícios)
real*8  :: l0               ! Auxiliar no calculo da norma inicial do residuo 
real*8  :: lnr              ! Auxiliar no calculo da norma do residuo por iteraçao
integer :: quociente  
integer :: resto
logical :: ver

integer, dimension(:),  allocatable :: i1, i2, i3           ! Numero do vertice do elemento triangular
integer, dimension(:),  allocatable :: necce1                ! Numero do elemento real contiguo ao contorno leste
integer, dimension(:),  allocatable :: neccn1               ! Numero do elemento real contiguo ao contorno norte
integer, dimension(:),  allocatable :: neccw1                ! Numero do elemento real contiguo ao contorno oeste
integer, dimension(:),  allocatable :: neccs1                ! Numero do elemento real contiguo ao contorno sul
integer, dimension(:),  allocatable :: necce2                ! Numero do elemento real contiguo ao contorno leste
integer, dimension(:),  allocatable :: neccn2                ! Numero do elemento real contiguo ao contorno norte
integer, dimension(:),  allocatable :: neccw2                ! Numero do elemento real contiguo ao contorno oeste
integer, dimension(:),  allocatable :: neccs2                ! Numero do elemento real contiguo ao contorno sul

integer, dimension(:,:),allocatable :: MCC                  ! Matriz de conectividade dos centróides		   
real*8,  dimension(:),  allocatable :: cp, v1, v2, v3, b	! Coeficientes
real*8,  dimension(:),  allocatable :: x, y				    ! Coordenadas dos pontos que geram elementos
real*8,  dimension(:),  allocatable :: xc, yc			    ! Coordenadas dos centroides do elemento
real*8,  dimension(:),  allocatable :: soma3			    ! Variavel auxiliar
real*8,  dimension(:),  allocatable :: Tpn				    ! Temperatura Numerica
real*8,  dimension(:),  allocatable :: Tpa				    ! Temperatura Analitica
real*8,  dimension(:),  allocatable :: R				    ! Norma residuo iteração inicial
real*8,  dimension(:),  allocatable :: Rr				    ! Residuo por iteração de cada elemento
real*8,  dimension(:),  allocatable :: normaerro		    ! Norma do erro de discretizaçao
real*8,  dimension(:,:),allocatable :: jotanormal		    ! Vetor fluxo de temperatura
real*8,  dimension(:),  allocatable :: nrr       		    ! Norma do residuo por iteração de elemento real


!********************************************************************
! Leitura da malha **************************************************
!********************************************************************

open(10,file = 'dados.ent')

read(10,*) Np
read(10,*) Ne
read(10,*) Ner
read(10,*) Nef

allocate( x(Np), y(Np) )
allocate( i1(Ne), i2(Ne), i3(Ne) )

do p = 1, Np
	read(10,*) x(p), y(p)
end do

do e = 1, Ne
	read(10,*) i1(e), i2(e), i3(e)
end do

close(10)

!********************************************************************
! Leitura da conectividade ******************************************
!********************************************************************

allocate( MCC(Ne,3) )

open(20,file = 'conectividades.ent')

do e = 1, Ne
	read(20,*) MCC(e,1), MCC(e,2), MCC(e,3)
end do

close(20)

!********************************************************************
! Leitura dos baricentros *******************************************
!********************************************************************

allocate ( xc(Ne), yc(Ne) )

open(30,file = 'baricentros.ent')

do e = 1, Ne
	read(30,*) xc(e), yc(e)
end do

!********************************************************************
! Area dos elementos ************************************************
!********************************************************************

areatotal = 1.0d0
area = 1.0d0 / Ner		 
!Calculo da métrica da malha  
h = (areatotal / Ner) ** (1.0d0 / 2.0d0)

!********************************************************************
! Cálculo dos coeficientes dos volumes de controle ******************
!********************************************************************

! Alocação da matriz jotanormal
allocate ( jotanormal(Ne,3) )

! Alocação do auxiliar soma
allocate ( soma3(Ne) )

! Alocação de vetores de coeficientes e termos-fontes
allocate ( cp(Ne), v1(Ne), v2(Ne), v3(Ne), b(Ne) )

! Volumes reais -----------------------------------------------------
! Inicialização de variáveis
jotanormal = 0.0d0

! Cálculo do produto jotanormal -------------------------------------
do e = 1, Ne         ! Varredura de todos os volumes
	do f = 1, 3		 ! Verificação apenas para os volumes conectados
		if ( MCC(e,f) /= 0 ) then
			! Caso 1
			if ( (i1(e) == i1(MCC(e,f))).and.(i2(e)==i3(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i1(e)) + x(i2(e)) )
				my = 0.5d0*( y(i1(e)) + y(i2(e)) )
				! Normal à aresta
				norx =   y(i2(e)) - y(i1(e))
				nory = -(x(i2(e)) - x(i1(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
			! Caso 2
			if ( (i1(e) == i1(MCC(e,f))).and.(i3(e)==i2(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i1(e)) + x(i3(e)) )
				my = 0.5d0*( y(i1(e)) + y(i3(e)) )
				! Normal à aresta
				norx =   y(i3(e)) - y(i1(e))
				nory = -(x(i3(e)) - x(i1(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
			! Caso 3
			if ( (i2(e) == i2(MCC(e,f))).and.(i3(e)==i1(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i2(e)) + x(i3(e)) )
				my = 0.5d0*( y(i2(e)) + y(i3(e)) )
				! Normal à aresta
				norx =   y(i3(e)) - y(i2(e))
				nory = -(x(i3(e)) - x(i2(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
			! Caso 4
			if ( (i2(e) == i2(MCC(e,f))).and.(i1(e)==i3(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i2(e)) + x(i1(e)) )
				my = 0.5d0*( y(i2(e)) + y(i1(e)) )
				! Normal à aresta
				norx =   y(i1(e)) - y(i2(e))
				nory = -(x(i1(e)) - x(i2(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
 			! Caso 5
			if ( (i3(e) == i3(MCC(e,f))).and.(i1(e)==i2(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i3(e)) + x(i1(e)) )
				my = 0.5d0*( y(i3(e)) + y(i1(e)) )
				! Normal à aresta
				norx =   y(i1(e)) - y(i3(e))
				nory = -(x(i1(e)) - x(i3(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
   			! Caso 6
			if ( (i3(e) == i3(MCC(e,f))).and.(i2(e)==i1(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i3(e)) + x(i2(e)) )
				my = 0.5d0*( y(i3(e)) + y(i2(e)) )
				! Normal à aresta
				norx =   y(i2(e)) - y(i3(e))
				nory = -(x(i2(e)) - x(i3(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
   			! Caso 7
			if ( (i1(e) == i2(MCC(e,f))).and.(i2(e)==i1(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i1(e)) + x(i2(e)) )
				my = 0.5d0*( y(i1(e)) + y(i2(e)) )
				! Normal à aresta
				norx =   y(i2(e)) - y(i1(e))
				nory = -(x(i2(e)) - x(i1(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
   			! Caso 8
			if ( (i1(e) == i3(MCC(e,f))).and.(i3(e)==i1(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i1(e)) + x(i3(e)) )
				my = 0.5d0*( y(i1(e)) + y(i3(e)) )
				! Normal à aresta
				norx =   y(i3(e)) - y(i1(e))
				nory = -(x(i3(e)) - x(i1(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
   			! Caso 9
			if ( (i2(e) == i3(MCC(e,f))).and.(i3(e)==i2(MCC(e,f))) ) then
				! Ponto médio da aresta
				mx = 0.5d0*( x(i2(e)) + x(i3(e)) )
				my = 0.5d0*( y(i2(e)) + y(i3(e)) )
				! Normal à aresta
				norx =   y(i3(e)) - y(i2(e))
				nory = -(x(i3(e)) - x(i2(e)))
				! Distância entre baricentros
				dbaric = dsqrt((xc(e) - xc(MCC(e,f)))**2 + (yc(e) - yc(MCC(e,f)))**2)
				! Vetor unitário entre baricentros
				vunitx = ( xc(e) - xc(MCC(e,f)) ) / dbaric
				vunity = ( yc(e) - yc(MCC(e,f)) ) / dbaric
				! Jota unitário
				Jx = vunitx / dbaric
				Jy = vunity / dbaric
				! Jota-normal
				Jotanormal(e,f) = -dabs(Jx*norx + Jy*nory)
			end if
		end if
	end do
end do

!********************************************************************
! Solução analítica *************************************************
!********************************************************************


allocate ( Tpa(Ne) )

do e = 1, Ne
	Tpa(e) = xc(e) * yc(e)
end do

! Volumes fictícios -------------------------------------------------
aux_fict = Nef / 4

! Contorno sul
do e = 1, aux_fict
	cp(Ner + e) =  1.0d0
	v1(Ner + e) = -1.0d0
	v2(Ner + e) =  0.0d0
	v3(Ner + e) =  0.0d0
	b(Ner + e)  =  0.0d0
end do

! Contorno leste
do e = 1, aux_fict
	cp(Ner + aux_fict + e) =  1.0d0
	v1(Ner + aux_fict + e) = -1.0d0
	v2(Ner + aux_fict + e) =  0.0d0
	v3(Ner + aux_fict + e) =  0.0d0 
	b(Ner + aux_fict + e)  =  2.0d0*yc(Ner + aux_fict + e)
end do

! Contorno norte
do e = 1, aux_fict
	cp(Ner + 2*aux_fict + e) =  1.0d0
	v1(Ner + 2*aux_fict + e) = -1.0d0
	v2(Ner + 2*aux_fict + e) =  0.0d0
	v3(Ner + 2*aux_fict + e) =  0.0d0
	b(Ner + 2*aux_fict + e)  =  2.0d0*xc(Ner + 2*aux_fict + e)
end do

! Contorno oeste
do e = 1, aux_fict
	cp(Ner + 3*aux_fict + e) =  1.0d0
	v1(Ner + 3*aux_fict + e) = -1.0d0
	v2(Ner + 3*aux_fict + e) =  0.0d0
	v3(Ner + 3*aux_fict + e) =  0.0d0
	b(Ner + 3*aux_fict + e)  =  0.0d0
end do

! Volumes reais -----------------------------------------------------
! Inicialização da soma
soma3 = 0.0d0
do e = 1, Ner
	do f = 1, 3
		soma3(e) = soma3(e) + jotanormal(e,f)
	end do	
end do

! Coeficientes e termos fontes --------------------------------------
do e = 1, Ner
	cp(e) = 1.0d0
	v1(e) = jotanormal(e,1) / soma3(e)
	v2(e) = jotanormal(e,2) / soma3(e)
	v3(e) = jotanormal(e,3) / soma3(e)
	b(e)  = 0.0d0
end do

!********************************************************************
! Solução numérica **************************************************
!********************************************************************


T1 = TIMEF()



open(10, file = "saidanormaerro.txt")
write(*,*) "Qual o nível da malha?"
read(*,*) malha
write(*,*) "Qual o numero de iteracoes desejadas para o método de Gauss Seidel?"
read(*,*) numiti
write(10,8)
write(10,9) malha
write(10,10)
allocate ( Tpn(Ne) )
allocate ( R(Ner))
allocate ( Rr(Ner))

allocate(normaerro(numiti))	
allocate(nrr(numiti))  				

Tpn = 0.0d0
!Calculo da norma da temperatura antes do processo iterativo

l0 = 0.0d0

do e = 1, Ner
	R(e) = v1(e)*Tpn(MCC(e,1)) + v2(e)*Tpn(MCC(e,2)) + v3(e)*Tpn(MCC(e,3)) + b(e) - cp(e)*Tpn(e)
	l0 = l0 + dabs(R(e))
end do 

if (l0 == 0.0d0) l0 = 1.0d0

do iti = 1, numiti
	somanorm = 0.0d0

	! Volumes reais

	lnr=0.0d0

	do e = 1, Ner
		Tpn(e) = ( v1(e)*Tpn(MCC(e,1)) + v2(e)*Tpn(MCC(e,2)) + v3(e)*Tpn(MCC(e,3)) + b(e) ) / cp(e)
	end do
	
	! Volumes fictícios

	do e = Ner+1, Ne
		soma_fict = 0.0d0
		if (MCC(e,1) /= 0) soma_fict = soma_fict + v1(e)*Tpn(MCC(e,1))
		if (MCC(e,2) /= 0) soma_fict = soma_fict + v2(e)*Tpn(MCC(e,2))
		if (MCC(e,3) /= 0) soma_fict = soma_fict + v3(e)*Tpn(MCC(e,3))
		Tpn(e) = (soma_fict + b(e)) / cp(e)
	end do

	do e = 1, Ner
		Rr(e) = v1(e)*Tpn(MCC(e,1)) + v2(e)*Tpn(MCC(e,2)) + v3(e)*Tpn(MCC(e,3)) + b(e) - cp(e)*Tpn(e)
		lnr = lnr + dabs(Rr(e))
		somanorm = somanorm + dabs(Tpa(e) - Tpn(e))
	end do
	
	nrr(iti) = lnr / l0
	normaerro(iti)= somanorm / Ner

	if(numiti >= 20000) then

		quociente = iti / 100
		resto = iti - (quociente * 100)
			
		if(resto == 0) then
			write(10,20) iti, normaerro(iti),nrr(iti)
		end if
    else

		write(10,20) iti, normaerro(iti),nrr(iti)

	end if



end do 
   
somatemp = 0.0d0

do e = 1, Ner
	somatemp = somatemp + (Tpn(e) * area)
end do 
 8 format(3/,"EQUAÇÃO DE LAPLACE LINEAR *** MALHAS TRIANGULARES VOLUMES FINITOS *** FICTÍCIOS ANTISSIMÉTRICOS")
 9 format(1/,"NÍVEL DA MALHA:",i3)
10 format (2/,"Iteração", 2x, "Norma Erro Discretização",11x,"Norma L1 do Residuo")
20 format (i8,2X, 1pe25.15,10x,1pe25.15)

!********************************************************************
!Temperatura Media **************************************************
!********************************************************************

Tmedia = somatemp / areatotal
qer = sqrt(Ner / 4.0d0)  !numero de volumes contiguos aos contornos
allocate(necce1(qer),neccn1(qer),neccw1(qer),neccs1(qer))
allocate(necce2(qer),neccn2(qer),neccw2(qer),neccs2(qer))

!********************************************************************
!Taxa de Transferencia de Calor nos contornos ***********************
!********************************************************************

drf = Ner / Nef

if (malha == 1) then

	necce = 2
	neccs = 1
	neccn = necce
	neccw = neccs

end if

if (malha /= 1)	then

necce1(1) = ((drf - 1.0d0) * 8.0d0) + 3.0d0
necce2(1) = necce1(1) + 1.0d0
neccn1(1) = (drf * (drf - 1.0d0) * 8.0d0) + 5.0d0
neccn2(1) = neccn1(1) + 1.0d0
neccw1(1) = 7.0d0
neccw2(1) = neccw1(1) + 1.0d0
neccs1(1) = 1.0d0
neccs2(1) = neccs1(1) + 1.0d0
 
passos = 8.0d0
passoe = Nef
passon = passos
passow = passoe


do e = 2, qer
	necce1(e) = necce1(e - 1) + passoe
	necce2(e) = necce1(e) + 1.0d0
	neccn1(e) = neccn1(e - 1) + passon
	neccn2(e) = neccn1(e) + 1.0d0
	neccw1(e) = neccw1(e - 1) + passow
	neccw2(e) = neccw1(e) + 1.0d0
	neccs1(e) = neccs1(e - 1) + passos
	neccs2(e) = neccs1(e) + 1.0d0
end do 

end if

if (malha == 1) then
	tte = yc(2) - Tpn(2)
	ttn = xc(2) - Tpn(2)
	ttw	= Tpn(1)
	tts	= Tpn(1) 		
end if


if (malha /= 1)	then

	somatte = 0.0d0 !calculo da taxa de transferencia ao leste
	somattn = 0.0d0	!calculo da taxa de transferencia ao norte
	somatts = 0.0d0	!calculo da taxa de transferencia ao sul
	somattw = 0.0d0	!calculo da taxa de transferencia ao oeste

	do e = 1, qer	
		somatte = somatte + (yc(necce1(e)) - Tpn(necce1(e)))
		somatte = somatte + (yc(necce2(e)) - Tpn(necce2(e))) 
		somattn = somattn + (xc(neccn1(e)) - Tpn(neccn1(e)))
		somattn = somattn + (xc(neccn2(e)) - Tpn(neccn2(e)))
		somattw = somattw + (Tpn(neccw1(e)) - yc(neccw1(e)))
		somattw = somattw + (Tpn(neccw2(e)) - yc(neccw2(e)))
		somatts = somatts + (Tpn(neccs1(e)) - xc(neccs1(e)))
		somatts = somatts + (Tpn(neccs2(e)) - xc(neccs2(e))) 	 
	end do 

	tte = - somatte / (1.0d0 - xc(necce1(1))) * ((1.0d0 / qer)/2.0d0)
	ttn = - somattn / (1.0d0 - yc(neccn1(1))) * ((1.0d0 / qer)/2.0d0)
	ttw = - somattw / xc(necce1(1)) * ((1.0d0 / qer)/2.0d0)
	tts = - somatts / yc(neccn1(1)) * ((1.0d0 / qer)/2.0d0)

end if



T2 = TIMEF()

TP = T2 - T1
  
!********************************************************************
! Saída de dados ****************************************************
!********************************************************************
		  
write(10,30)
do e = 1, Ne
	write(10,40) e, xc(e), yc(e), Tpa(e), Tpn(e), Tpa(e)-Tpn(e)
end do

write(10,50) Tmedia
write(10,60) tte,ttn,ttw,tts
write(10,70) numiti
write(10,80) TP

30 format (3/,4x,"Volume", t13, "Posição x", t38, "Posição y", t63, "Sol. analítica", t88, "Sol. numérica", t113, "Erro")
40 format ( i8,<5>(1pe25.15) )
50 format (3/,"Temperatura Média:",1x, 1pe25.15)
60 format (1//,"Taxa de Transferencia de Calor no Contorno Leste:", 1pe25.15,3x,&
           1/, "Taxa de Transferencia de Calor no Contorno Norte:", 1pe25.15,3x,&
		   1/, "Taxa de Transferencia de Calor no Contorno Oeste:", 1pe25.15,3x,&
		   1/, "Taxa de Transferencia de Calor no Contorno Sul:",2x,1pe25.15)

70 format(4/,"Numero de iterações =",i10)

80 format(1/,"Tempo de Processamento =",1pe25.6,2x, "segundos")

close(10)
ver = system('Notepad saidanormaerro.txt')

!ARQUIVO PARA SAIDAS RICHARDSON	
!Temperatura em x=y=1/2

T = 0.0d0
cont = 0.0d0

if (Ner == 2.0d0) then			
		T = (Tpn(1) + Tpn(2))/2.0d0

	else
	do e = 1, Ner
		if((x(i1(e)) == 0.5d0 .and. y(i1(e)) == 0.5d0) .or. (x(i2(e)) == 0.5d0 .and. y(i2(e)) == 0.5d0) .or. &
		   (x(i3(e)) == 0.5d0 .and. y(i3(e)) == 0.5d0))  then
			cont = cont + 1
			T = T + Tpn(e)
		end if
	end do 
		T = T / cont
end if

open(30, file = "saidasrichardson.txt")
write(30,71) h
write(30,72) T
write(30,73) Tmedia
write(30,74) tte
write(30,75) ttn
write(30,76) ttw
write(30,77) tts
write(30,78) normaerro(numiti)
close(30)

71 format (2x,1pe24.15,2x, "(Métrica da Malha m2)")
72 format (" 1",1pe24.15,2x, "(Temperatura em x=y=1/2 em °C)")
73 format (" 2",1pe24.15,2x, "(Temperatura Média)")
74 format (" 3",1pe24.15,2x, "(Taxa de Transferencia de Calor no Contorno Leste em °C/m)")
75 format (" 4",1pe24.15,2x, "(Taxa de Transferencia de Calor no Contorno Norte em °C/m)")
76 format (" 5",1pe24.15,2x, "(Taxa de Transferencia de Calor no Contorno Oeste em °C/m)")
77 format (" 6",1pe24.15,2x, "(Taxa de Transferencia de Calor no Contorno Sul em °C/m)")
78 format (" 7",1pe24.15,2x, "(Media da Norma L1 do Erro Numerico)")  

ver = system('Notepad saidasrichardson.txt')

end program


