 module Coef
 
 use Variaveis 
 
 implicit none

 contains

	subroutine Coeficientes

    real*8  :: x, y,u_x,u_xx,u_y,u_yy,v_x,v_xx,v_y,v_yy                ! coordenada espacial
	
	dx = 1.0d0/(Nx-1)
	dy = 1.0d0/(Ny-1)

	!i é x
	!j é y

    ! pontos internos

	if (vel==1) then

		select case (caso_ang)
					
		case(0)
			
			!Para o ângulo de 0 graus
			u=1.0d0
			vv=0.0d0

		case(30)

			!Para o angulo de 30 graus
			u=dsqrt(3.0d0)/2.0d0
			vv=1.0d0/2.0d0

		case(3)
				
			!Para velocidades nulas: EQUAÇÂO DE LAPLACE
			u=0.0d0
			vv=0.0d0

		case(45)

			select case(caso_mod)

				case(0)

				!Para velocidades constantes e iguais a 1
				u=1.0d0
				vv=1.0d0
				
				case(1)

				!Para o módulo 1
				
				u=dsqrt(2.0d0)/2.0d0
				vv=dsqrt(2.0d0)/2.0d0

				case(2)

				!Para o modulo 2

				u=3.0d0*dsqrt(2.0d0)/2.0d0
				vv=3.0d0*dsqrt(2.0d0)/2.0d0

				case(10)
				!Para o módulo 10
				
				u=dsqrt(5.0d1)
				vv=dsqrt(5.0d1)

				case(100)
				!Para o módulo 100

				u=dsqrt(5.0d3)
				vv=dsqrt(5.0d3)

				case(1000)
				!Para o módulo 1000

				u=dsqrt(5.0d5)
				vv=dsqrt(5.0d5)
		
		end select
						 

		case(60)

			!Para o angulo de 60 graus
			u=1.0d0/2.0d0
			vv=dsqrt(3.0d0)/2.0d0

		case(90)

			!Para o angulo de 90 graus
			u=0.0d0
			vv=1.0d0

		end select
			

		do i=2,Nx-1

			x=(i-1)*dx

			do j=2,Ny-1

				y=(j-1)*dx
				np = (j-1)*Nx+i

				AS(np) = - ( ( alphas / ( dy ** 2 ) ) + ( vv / dy ) )
				AW(np) = - ( ( alphas / ( dx ** 2 ) ) + ( u / dx ) )
				AP(np) = ( 2.0d0 * alphas / ( dx ** 2 ) ) + ( 2.0d0 * alphas / ( dy ** 2 ) ) + &
				( u / dx ) + ( vv / dy )
				AE(np) = - alphas / ( dx ** 2 )
				AN(np) = - alphas / ( dy ** 2 )
				BP(np) = (u**2)*vv*y+(vv**2)*u*x

			end do

		end do
	
	
	else ! Para as velocidades de Shih

		do i = 2,Nx-1

			x = (i-1)*dx

			do j = 2,Ny-1

				y = (j-1)*dy

				! Velocidades para o termo fonte

				u = 8.0d0 * ((x**4)-2.0d0*(x**3)+x**2) * (4.0d0*(y**3)-2.0d0*y)
				vv = -8.0d0 * (4.0d0*(x**3)-6.0d0*(x**2)+2.0d0*x)*(y**4-y**2)

				u_x = 32.0d0*(2.0d0*x**3-3.0d0*x**2+x) * (2.0d0*y**3-y)
				u_xx = 32.0d0 * (6.0d0*x**2-6.0d0*x+1.0d0)*(2.0d0*y**3-y)
				u_y = 16.0d0 * (x**4-2.0d0*x**3+x**2)*(6.0d0*y**2-1.0d0)
				u_yy = 192.0d0*y*(x**4-2.0d0*x**3+x**2)
				
				v_x = -16.0d0 * (6.0d0*x**2-6.0d0*x+1.0d0)*(y**4-y**2)
				v_xx = -96.0d0*(2.0d0*x-1.0d0)*(y**4-y**2)
				v_y = -32.0d0*(2.0d0*x**3-3.0d0*x**2+x)*(2.0d0*y**3-y)
				v_yy = -32.0d0*(2.0d0*x**3-3.0d0*x**2+x)*(6.0d0*y**2-1.0d0)
						
				np = (j-1)*Nx+i

				! Coeficientes
					
				AS(np) = - ( ( alphas / ( dy ** 2 ) ) + ( vv / dy ) )
				AW(np) = - ( ( alphas / ( dx ** 2 ) ) + ( u / dx ) )
				AP(np) = ( 2.0d0 * alphas / ( dx ** 2 ) ) + ( 2.0d0 * alphas / ( dy ** 2 ) ) + &
				( u / dx ) + ( vv / dy )
				AE(np) = - alphas / ( dx ** 2 )
				AN(np) = - alphas / ( dy ** 2 )
				
				BP(np) = u*y*(x*(u_x*vv+u*v_x)+u*vv)+vv*x*(y*(u_y*vv+u*v_y)+u*vv)-&
				alphas*y*(x*(u_xx*vv+2.0d0*u_x*v_x+u*v_xx)+2.0d0*(u_x*vv+v_x*u))-&
				alphas*x*(y*(u_yy*vv+2.0d0*u_y*v_y+u*v_yy)+2.0d0*(u_y*vv+u*v_y))
					
			end do

		end do

	end if

	
	! contorno norte

	pi=dacos(-1.0d0)

	j = Ny

	do i = 1,Nx

		x = (i-1)*dx
		y = (j-1)*dy
		
		np = (j-1)*Nx+i
				
		AS(np) = 0.0d0
		AW(np) = 0.0d0
		AP(np) = 1.0d0
		AE(np) = 0.0d0
		AN(np) = 0.0d0		
		
		if(vel==1) then  ! Velocidades constantes


			if (caso_ang==3) then !Equação de Laplace

				BP(np)=dsin(pi*x)

			else


				BP(np) = u*vv*x*y 

			end if

		else !Velocidades de Shih

			BP(np)=0.0d0

		end if
		
			
	end do

	
	! contorno sul

    j = 1

	do i = 1,Nx

		x = (i-1)*dx

		np = (j-1)*Nx + i
				
		AS(np) = 0.0d0
		AW(np) = 0.0d0
		AP(np) = 1.0d0
		AE(np) = 0.0d0
		AN(np) = 0.0d0
		BP(np)=0.0d0
		
	end do

	! contorno oeste

    i = 1

	do j = 1,Ny


		x = (i-1)*dx

		np = (j-1)*Nx + i

		AS(np) = 0.0d0
		AW(np) = 0.0d0
		AP(np) = 1.0d0
		AE(np) = 0.0d0
		AN(np) = 0.0d0
		BP(np)=0.0d0
					
	end do
    
	! contorno leste

	i = Nx

	do j = 1,Ny

		x = (i-1)*dx
		y = (j-1)*dy
	
		np = (j-1)*Nx+i

		AS(np) = 0.0d0
		AW(np) = 0.0d0
		AP(np) = 1.0d0
		AE(np) = 0.0d0
		AN(np) = 0.0d0

		if(vel==1) then

			BP(np) = u*vv*x*y 

		else

			BP(np)=0.0d0
		
		end if
		
	end do

  end subroutine Coeficientes

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

subroutine Estimativa_inicial

		real*8  :: x,y                ! coordenada espacial
		
		
		dx = 1.0d0/(Nx-1)
		dy = 1.0d0/(Ny-1)

		select case (ci)
						
			case (1) ! Estimativa inicial nula

		
			!i é x
			!j é y

			! pontos internos

			
			TNum = 0.0d0

			! contorno norte
		
								
				j = Ny

				do i = 1,Nx
					
					y = (j-1)*dy
					x = (i-1)*dx
					np = (j-1)*Nx+i
						
					if(vel==1) then  !Velocidades Constantes

						if (caso_ang==3) then

							TNum(np)=dsin(pi*x)

						else

						TNum(np) = u*vv*x*y

						end if

					else
						
						TNum(np)=0.0d0

					end if
																
				end do


			! contorno leste

				pi=dacos(-1.0d0)
				
				i = Nx

				do j = 1,Ny

					x = (i-1)*dx
			
					np = (j-1)*Nx+i
				
					y = (j-1)*dy
					
					if(vel==1) then

						TNum(np) = u*vv*y*x

					else

						TNum = 0.0d0
			
		
					end if
															
				end do


			case (2)  ! Estimativa inicial é a solução analítica

			!Pontos internos

			
			if (vel==1) then   !Velocidades constantes

				if (caso_ang==3) then !Problema de Laplace

					do i=1,Nx
					
						x = (i-1)*dx
					
						do j=1,Ny

							y = (j-1)*dy
							
							np = (j-1)*Nx + i
											
							TA(np)=dsin(pi*x)*(dsinh(pi*y)/dsinh(pi))

						end do
									
					end do

				else !Velocidades Constantes

					do i=1,Nx
					
						x = (i-1)*dx
					
						do j=1,Ny

							y = (j-1)*dy
							
							np = (j-1)*Nx + i
											
							TA(np)=u*vv*x*y					

						end do
									
					end do

				end if

			else !Velocidades de Shih

				do i=1,Nx
				
					x = (i-1)*dx
				
					do j=1,Ny

						y = (j-1)*dy
						
						np = (j-1)*Nx + i
										
						u = 8.0d0 * ((x**4) - 2.0d0*(x**3) + x**2) * (4.0d0*(y**3) - 2.0d0*y)
						vv = -8.0d0 * (4.0d0*(x**3) - 6.0d0*(x**2) + 2.0d0*x) * (y**4 - y**2)
						
						TA(np)=u*vv*x*y
					

					end do
								
				end do

			end if
			
			TNum = TA

			!contorno norte		
								
			if (vel==1) then

				if(caso_ang==3) then

					do i=1,Nx
					
						x = (i-1)*dx
						
						do j=1,Ny

							y = (j-1)*dy
								
							np = (j-1)*Nx + i
												
							TNum(np)=dsin(pi*x)

						end do
									
					end do

				else ! Velocidades Constantes
						
					j = Ny

					do i = 1,Nx

						y = (j-1)*dy
						x = (i-1)*dx
						np = (j-1)*Nx+i
								
						TNum(np) = u*vv*x*y
												
					end do

				end if

				else ! Velocidades de Shih

					j = Ny

					do i = 1,Nx

						y = (j-1)*dy
						x = (i-1)*dx
						np = (j-1)*Nx+i
								
						u = 8.0d0 * ((x**4) - 2.0d0*(x**3) + x**2) * (4.0d0*(y**3) - 2.0d0*y)
						vv = -8.0d0 * (4.0d0*(x**3) - 6.0d0*(x**2) + 2.0d0*x) * (y**4 - y**2)
						
						TNum(np)=u*vv*x*y
																										
					end do

				end if


				! contorno leste
	
				pi=dacos(-1.0d0)
				
				if(vel==1) then ! Velocidades Constantes

					i = Nx

					do j = 1,Ny

						x = (i-1)*dx
				
						np = (j-1)*Nx+i
					
						y = (j-1)*dy
									
						TNum(np) = u*vv*y*x

					end do

				else !Velocidades de Shih
				
					i = Nx

					do j = 1,Ny

						x = (i-1)*dx
				
						np = (j-1)*Nx+i
					
						y = (j-1)*dy
						
						u = 8.0d0 * ((x**4) - 2.0d0*(x**3) + x**2) * (4.0d0*(y**3) - 2.0d0*y)
						vv = -8.0d0 * (4.0d0*(x**3) - 6.0d0*(x**2) + 2.0d0*x) * (y**4 - y**2)
					
						TNum(np)=u*vv*x*y
																			
					end do

				end if


			case (3) ! Estimativa inicial igual a 0.5


			!i é x
			!j é y

			! pontos internos

			
			TNum = 0.5d0

			! contorno norte

			pi=dacos(-1.0d0)
										
				if (vel==1) then

					if (caso_ang==3) then ! Para a Equação de Laplace

						j = Ny
						do i = 1,Nx
						
							y = (j-1)*dy
							x = (i-1)*dx
							np = (j-1)*Nx+i
								
							Tnum(np)=dsin(pi*x)

						end do
					
					else ! Para velocidades constantes

						j = Ny
						do i = 1,Nx
							
							y = (j-1)*dy
							x = (i-1)*dx
							np = (j-1)*Nx+i
								
							TNum(np) = u*vv*x*y

						end do

					end if

				else ! Para velocidades de Shih

					TNum = 0.5d0
						
				end if
						
		
			! contorno leste

				pi=dacos(-1.0d0)
				
				i = Nx

				do j = 1,Ny

					x = (i-1)*dx
			
					np = (j-1)*Nx+i
				
					y = (j-1)*dy
					
					if(vel==1) then

						TNum(np) = u*vv*y*x

					else

						TNum = 0.5d0
			
		
					end if
															
				end do

		
		end select

		
end subroutine Estimativa_inicial

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

	subroutine Condicao_contorno

		real*8  :: x,y                ! coordenada espacial
				
		dx = 1.0d0/(Nx-1)
		dy = 1.0d0/(Ny-1)

		!i é x
		!j é y

		! contorno norte

		pi=dacos(-1.0d0)

		if (vel==1) then

			if (caso_ang==3) then ! Para a Equação de Laplace
				
				j = Ny

				do i = 1,Nx
				
					x = (i-1)*dx
					y = (j-1)*dy
							
					np = (j-1)*Nx+i

					T_r(np)=dsin(pi*x)

				end do
			
			else ! Para as velocidades constantes
				
				j = Ny

				do i = 1,Nx
						
					x = (i-1)*dx
					y = (j-1)*dy
					np = (j-1)*Nx+i

					T_r(np) = u*vv*x*y 

				end do

			end if

		else ! Para as velocidades de Shih

			j = Ny

				do i = 1,Nx
				
					x = (i-1)*dx
					y = (j-1)*dy
							
					np = (j-1)*Nx+i

					T_r(np)=0.0d0

				end do
			
		
		end if
			
		
		! contorno sul

		j = 1

		do i = 1,Nx

			np = (j-1)*Nx + i

			T_r(np)=0.0d0
		
		end do

		! contorno oeste

		i = 1

		do j = 1,Ny
	
			np = (j-1)*Nx + i

			T_r(np)=0.0d0
		
		end do
    
		! contorno leste

		i = Nx

		do j = 1,Ny

			x = (i-1)*dx
			y = (j-1)*dy
	
			np = (j-1)*Nx+i

			if(vel==1) then

				T_r(np) = u*vv*x*y 

			else

				T_r(np)=0.0d0
		
			end if

					
		end do

										
	end subroutine Condicao_contorno

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

end module Coef