module equilibrio

! objetivo: processar o algoritmo principal do código.

! Versão atual:  28 Out 06

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

use usuario
use solvers
use Gibbs_1p4

implicit none

contains

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

  subroutine NUM_EQ_cp_gama_congelados ( nx, ny, T, p, cp, gcp, Rgp, Yi, hi )

    integer,intent(in) :: nx ! número de volumes de controle com fictícios (direção axial)
    integer,intent(in) :: ny ! número de volumes de controle com fictícios (direção radial)

    real*8,dimension(:),  intent(in)    :: T     ! temperatura (K)
    real*8,dimension(:),  intent(in)    :: P     ! pressao
    
    real*8,dimension(:),  intent(inout) :: cp    ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),  intent(inout) :: gcp   ! razão de calores específicos (adimensional)
    real*8,dimension(:),  intent(inout) :: Rgp   ! constante da mistura de gases (J/kg)
    
    real*8,dimension(:,:),intent(inout) :: Yi    ! fração mássica (adim.)
    real*8,dimension(:,:),intent(inout) :: hi    ! entalpia (J/kg.K)
    
    integer :: nps
    
    do i = 2, nx-1
		do j = 1, ny

	       np = (j-1)*nx + i
		   
	       call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax,  &
											    tol_ej, tol_n, T(np), p(np)/101325.0d0, OF )
	       call GIBBS_EQUILIBRIO_calculos1
	       call GIBBS_hi_massa_calculo
    	   call GIBBS_MACH1D_envia4 ( cp(np), gcp(np), Rgp(np), Yi(np,:), hi(np,:) )

		end do
    end do
    
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + 1
       nps = np - nx
       
       cp(np)   = cp(nps)
       gcp(np)  = gcp(nps)
       Rgp(np)  = Rgp(nps)
       hi(np,:) = hi(nps,:)
       
    end do
    

	i = 1
	do j = 2, ny-1

       np = (j-1)*nx + i

       call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax,  &
										    tol_ej, tol_n, T(np), p(np)/101325.0d0, OF )
       call GIBBS_EQUILIBRIO_calculos1
       call GIBBS_hi_massa_calculo
	   call GIBBS_MACH1D_envia4 ( cp(np), gcp(np), Rgp(np), Yi(np,:), hi(np,:) )

	end do

	i = nx
	do j = 2, ny-1

       np = (j-1)*nx + i

       call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax,  &
										    tol_ej, tol_n, T(np), p(np)/101325.0d0, OF )
       call GIBBS_EQUILIBRIO_calculos1
       call GIBBS_hi_massa_calculo
	   call GIBBS_MACH1D_envia4 ( cp(np), gcp(np), Rgp(np), Yi(np,:), hi(np,:) )
	   
	end do
	
  end subroutine NUM_EQ_cp_gama_congelados

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

  subroutine NUM_EQ_gama_efetivo ( nx, ny, T, p, gcpf, Rgp, gcp )

    integer,intent(in) :: nx ! número de volumes de controle com dois fictícios (direção x)
    integer,intent(in) :: ny ! número de volumes de controle com dois fictícios (direção y)

    real*8,dimension(:),intent(in)    :: T     ! temperatura (K)
    real*8,dimension(:),intent(in)    :: p     ! pressão (Pa)
    real*8,dimension(:),intent(in)    :: gcpf  ! gama congelado (adimensional)
    real*8,dimension(:),intent(in)    :: Rgp   ! constante da mistura de gases (J/kg.K)

    real*8,dimension(:),intent(inout) :: gcp   ! razão de calores específicos (adimensional)

    integer :: nps

    do i = 2, nx-1
		do j = 1, ny

		   np = (j-1)*nx + i
           
		   call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax, tol_ej, &
												tol_n, T(np), (p(np)+dpY*p(np)/2.0d0)/101325.0d0, OF )
		   call GIBBS_EQUILIBRIO_calculos1
		   call GIBBS_MACH1D_reativo5
		   call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax, tol_ej, &
												tol_n, T(np), (p(np)-dpY*p(np)/2.0d0)/101325.0d0, OF )
		   call GIBBS_EQUILIBRIO_calculos1
		   call GIBBS_MACH1D_reativo6
		   call GIBBS_MACH1D_reativo4 ( dpY*p(np), gcpf(np), 8314.51d0/Rgp(np), p(np), gcp(np) )

		end do
    end do

    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + 1
       nps = np - nx
       
       gcp(np)  = gcp(nps)
       
    end do

	i = 1
	do j = 2, ny-1

	   np = (j-1)*nx + i

	   call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax, tol_ej, &
											tol_n, T(np), (p(np)+dpY*p(np)/2.0d0)/101325.0d0, OF )
	   call GIBBS_EQUILIBRIO_calculos1
	   call GIBBS_MACH1D_reativo5
	   call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax, tol_ej, &
											tol_n, T(np), (p(np)-dpY*p(np)/2.0d0)/101325.0d0, OF )
	   call GIBBS_EQUILIBRIO_calculos1
	   call GIBBS_MACH1D_reativo6
	   call GIBBS_MACH1D_reativo4 ( dpY*p(np), gcpf(np), 8314.51d0/Rgp(np), p(np), gcp(np) )
	
	end do

	i = nx
	do j = 2, ny-1

	   np = (j-1)*nx + i

	   call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax, tol_ej, &
											tol_n, T(np), (p(np)+dpY*p(np)/2.0d0)/101325.0d0, OF )
	   call GIBBS_EQUILIBRIO_calculos1
	   call GIBBS_MACH1D_reativo5
	   call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modq, itimax, itemax, tol_ej, &
											tol_n, T(np), (p(np)-dpY*p(np)/2.0d0)/101325.0d0, OF )
	   call GIBBS_EQUILIBRIO_calculos1
	   call GIBBS_MACH1D_reativo6
	   call GIBBS_MACH1D_reativo4 ( dpY*p(np), gcpf(np), 8314.51d0/Rgp(np), p(np), gcp(np) )
	
	end do

    
  end subroutine NUM_EQ_gama_efetivo

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

  subroutine NUM_EQ_Mach ( Rgp, gcp, u, T, M )

    real*8,dimension(:),intent(in)  :: Rgp  ! constante do gás (J/kg.K)
    real*8,dimension(:),intent(in)  :: gcp  ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in)  :: u    ! velocidade (m/s)
    real*8,dimension(:),intent(in)  :: T    ! temperatura (K)

    real*8,dimension(:),intent(out) :: M    ! número de Mach (adimensional)

    M = u / dsqrt ( gcp * Rgp * T )
  
  end subroutine NUM_EQ_Mach

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

  subroutine NUM_EQ_mi_k_congelados ( nxy, T, Yi, vlp, kp )

    integer,intent(in) :: nxy ! número de volumes de controle com dois fictícios

    real*8,dimension(:),intent(in)    :: T   ! temperatura (K)
    
    real*8,dimension(:,:),intent(in)  :: Yi  ! fração mássica (adim.)

    real*8,dimension(:),intent(inout) :: vlp ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),intent(inout) :: kp  ! condutividade térmica (W/m.K)
    
    if ( tipo_vl == 2 .or. tipo_k  == 2 ) then
       vlp = 0.0d0
       kp  = 0.0d0
          
       do i = 1, nxy
          

          call GIBBS_MACH1D_recebe2 ( T(i) )

          call GIBBS_CONGELADO_calculos2
          
          call GIBBS_CONGELADO_envia2 ( vlp(i), kp(i) )

       end do

    end if

    if ( tipo_vl == 1 ) vlp = visc_cte

    if ( tipo_k  == 1 ) kp  = k_cte
   
  end subroutine NUM_EQ_mi_k_congelados

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

end module equilibrio