module congelado

! 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 solucao_congelado

    integer :: unit_it = 81 ! identificador do arquivo das iterações
    integer :: unit_au = 82 ! identificador do arquivo dos campos auxiliares

	integer :: ii          ! auxiliar
	character(80) :: aux   ! auxiliar

! *** campos auxiliares ***

    real*8,dimension(:),allocatable :: C_rec  ! coeficiente de recuperação de Taw
    real*8,dimension(:),allocatable :: Taw    ! temperatura de parede adiabática (K)
    real*8,dimension(:),allocatable :: qlt    ! fluxo de calor total para a parede (W/m2)
    real*8,dimension(:),allocatable :: qlr    ! fluxo de calor radiativo para a parede (W/m2)

    real*8,dimension(:),allocatable :: f      ! fator de atrito de Darcy (adim)
    real*8,dimension(:),allocatable :: h      ! coeficiente de convecção (W/m2.K)
    real*8,dimension(:),allocatable :: mi     ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),allocatable :: k      ! condutividade térmica (W/m.K)

    ! alocação de memória  
    allocate ( f(nxy),   h(nxy),   mi(nxy),  k(nxy)  )

    if ( tipo_cp == 2 ) then
       call GIBBS_EQUILIBRIO_dados_recebe ( unit_au, 0, modelo_q, itimax, itemax,  &
                                            tol_e, tol_n, T0, p0/101325.0d0, OF )
       call GIBBS_EQUILIBRIO_inicio1
       call GIBBS_EQUILIBRIO_calculos1
       call GIBBS_MACH1D_envia2 ( Ne )
       allocate ( Yi(nxy,Ne), especie(Ne) )
       call GIBBS_MACH1D_envia3 ( especie )
       call GIBBS_MACH1D_envia6 ( cp_cte, gama_cte, Rg_cte, Yi(1,:) )
       do i = 2, nxy
          Yi(i,:) = Yi(1,:)
       end do
       call NUM_CONG_cp_gama_congelados ( nxy, T, Yi, cp, gama )
    else
       allocate ( Yi(nxy,1), especie(1) )
       Yi = 0.0d0
       Ne = 1
       especie = 'nada'
       cp   = cp_cte
       gama = gama_cte
    end if
	Rg = Rg_cte

    ! inicialização das propriedades
    call GIBBS_mi_kc_inicio
    call NUM_CONG_Mach ( Rg, gama, u, T, M )

    ! zeragem do tempo de CPU

    tcpu = timef()

    ! ciclo iterativo externo no tempo

    do it = 1, itmax

       ! zeragem da norma a cada iteração
       norma = 0.0d0

       ! atualização das incógnitas para novo incremento de tempo
       ua   = u
       va   = v
       pa   = p
       Ta   = T
       roa  = ro
       Ucea = Uce
       Vcna = Vcn

       ! atualização das condições na entrada da tubeira
       pina = pin
       call Tin_pin_entrada
       call plin_entrada_p_ficticio

       ! cálculo das propriedades

       if ( it >= it_q .and. ( it == itmax .or. mod(it,w_q)==0 ) ) then

          if ( tipo_cp == 2 ) call NUM_CONG_cp_gama_congelados ( nxy, T, Yi, cp, gama )

          call NUM_CONG_mi_k_congelados ( nxy, T, mi, k )

       end if

       ! solução de u
       call coeficientes_gerais
       call fonte_velocidade_u
       call ccu
       call lu2d5(a,dl,du,nxy,nx,ny)
       call fb2d5(a,dl,du,ny,nx,nxy,u,tolu,nitm_u,b)
!	   call gauss_seidel(a,u,b,nx,ny,nxy)
       call norma_L1 (nx,ny,nxy,a,b,u,norma)

       ! solução de v
       call fonte_velocidade_v
       call ccv
       call lu2d5(a,dl,du,nxy,nx,ny)
       call fb2d5(a,dl,du,ny,nx,nxy,v,tolu,nitm_u,b)
!	   call gauss_seidel(a,v,b,nx,ny,nxy)
       call norma_L1 (nx,ny,nxy,a,b,v,norma)

       ! solução de T
       call fonte_temperatura
       call ccT
       call lu2d5(a,dl,du,nxy,nx,ny)
       call fb2d5(a,dl,du,ny,nx,nxy,T,tolu,nitm_u,b)
!	   call gauss_seidel(a,T,b,nx,ny,nxy)
       call norma_L1 (nx,ny,nxy,a,b,T,norma)

       ! cálculo da massa especifica
       call massa_especifica_nos
       call massa_especifica_faces

       ! cálculo dos coeficientes do SIMPLEC e das velocidades contravariantes nas faces
       call coeficientes_SIMPLEC
       call velocidades_contravariantes_faces
       call Uce_Vcn_de_dn_faces_contornos

       ! ciclo iterativo da massa
       
	   do im = 1, imax

          ! solução de pl
          call coeficientes_fonte_pressao
          call ccp
          call lu2d5(a,dl,du,nxy,nx,ny)
          call fb2d5(a,dl,du,ny,nx,nxy,pl,tolp,nitm_p,b)
!		  call gauss_seidel(a,pl,b,nx,ny,nxy)
       
          ! correção de p, u, v, Uce, Vcn, ro, roe e ron com pl
          call pressao_e_massa_especifica_com_pl
          call massa_especifica_faces
          call u_v_nos_reais_com_pl
          call u_v_nos_ficticios_com_pl
          call Uce_Vcn_faces_internas_com_pl
          call Uce_Vcn_faces_contornos_com_pl

       end do ! fim do ciclo iterativo da massa

       call norma_L1 (nx,ny,nxy,a,b,pl,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

    ! finalização do tempo de CPU

    tcpu = timef()

    ! cria o arquivo dos campos auxiliares

       if ( tipo_cp > 1 ) then
          call GIBBS_EQUILIBRIO_dados_recebe ( unit_au, 0, modelo_q, itimax, itemax,  &
                                               tol_e, tol_n, T0, p0/101325.0d0, OF )
          call GIBBS_EQUILIBRIO_calculos1
          call GIBBS_si_aloca_memoria
          call GIBBS_MACH1D_camara_entropia ( entin, ny )
          call GIBBS_EQUILIBRIO_dados_escreve ( unit_au )
          call GIBBS_EQUILIBRIO_escreve1
          call GIBBS_CONGELADO_calculos2
          call GIBBS_CONGELADO_escreve2 ( unit_au )
          call NUM_CONG_cp_gama_congelados ( nxy, T, Yi, cp, gama )
          do i = 1, nxy
             call GIBBS_MACH1D_entropia ( T(i), p(i)/101325.0d0, Rg(i), Yi(i,:), ent(i) )
          end do

       end if

    close(8) ! fecha arquivo das iterações

    close(11) ! fecha arquivo dos resíduos

  end subroutine solucao_congelado

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

  subroutine NUM_CONG_cp_gama_congelados ( nxy, T, Yi, cp, gama )

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

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

    real*8,dimension(:),  intent(inout) :: cp    ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),  intent(inout) :: gama  ! razão de calores específicos (adimensional)

    do i = 1, nxy

	    call GIBBS_EQUILIBRIO_dados_recebe ( 1, 0, modelo_q, itimax, itemax,  &
											tol_e, tol_n, T(i), p0/101325.0d0, OF )
	    call GIBBS_MACH1D_recebe1 ( Yi(i,:) )

	    call GIBBS_cp_calculo

	    call GIBBS_gama_calculo

	    call GIBBS_MACH1D_envia1 ( cp(i), gama(i) )

    end do
    
  end subroutine NUM_CONG_cp_gama_congelados

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

  subroutine NUM_CONG_fr_C_rec ( gama, Pr, M, fr, C_rec )

    real*8,dimension(:),intent(in)  :: gama  ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in)  :: Pr    ! número de Prandtl (adim)
    real*8,dimension(:),intent(in)  :: M     ! número de Mach (adimensional)

    real*8,dimension(:),intent(out) :: fr    ! fator de recuperação (adim)
    real*8,dimension(:),intent(out) :: C_rec ! coeficiente de recuperação de Taw

    select case ( tipo_fr )

       case ( 1 )
          fr = fr_cte

       case ( 2 )
          fr = Pr**(1/3.0d0)

    end select

    C_rec = 1.0d0 + fr*(gama-1.0d0)*(M**2)/2.0d0
    
  end subroutine NUM_CONG_fr_C_rec

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

  subroutine NUM_CONG_h ( n, fm, rc, Sg, sp, M, gama, cp, mi, T_wall, Pr, h )

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

    real*8, intent(in) :: fm ! fluxo de massa numérico (kg/s)
    real*8, intent(in) :: rc ! raio de curvatura na garganta da tubeira (m)
    real*8, intent(in) :: Sg ! área da garganta (m2)

    real*8,dimension(:),intent(in) :: sp     ! área no centro do volume de controle
    real*8,dimension(:),intent(in) :: M      ! número de Mach (adimensional)
	real*8,dimension(:),intent(in) :: gama   ! razão de calores específicos (adimensional)
    real*8,dimension(:),intent(in) :: cp     ! calor específico à pressão constante (J/kg.K)
    real*8,dimension(:),intent(in) :: mi     ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),intent(in) :: T_wall ! temperatura da parede da tubeira (K)
    real*8,dimension(:),intent(in) :: Pr     ! número de Prandtl (adim)

    real*8,dimension(:),intent(inout) :: h   ! coeficiente de convecção (W/m2.K)

    real*8 :: Dg ! diâmetro da garganta da tubeira (m)
    real*8 :: pi ! constante matemática
    real*8 :: cte1, cte2, cte3, cte4, cte5, cte6 ! auxiliares

    select case ( tipo_h )

       case ( 1 )
          h = h_cte

       case ( 2 )
          pi = dacos(-1.0d0)
          Dg = 4.0d0 * Sg / pi
          cte1 = 0.026d0*((mi(1)**0.2d0)*cp(1)/(Pr(1)**0.6d0))/(Dg**0.2d0)
          cte2 = ((fm/Sg)**0.8d0)*((Dg/rc)**0.1d0)
          h(1) = 0.0d0
          do i = 2, n-1
             cte3 = 1.0d0 + 0.5d0*(gama(i)-1.0d0)*(M(i)**2)
             cte4 = (0.5d0*(1.0d0+T_wall(i)*cte3/T0))**0.68d0
             cte5 = cte3**0.12d0
             cte6 = 1.0d0 / ( cte4 * cte5 )
             h(i) = cte1*cte2*((Sg/sp(i))**0.9d0)*cte6
          end do
          h(n) = 0.0d0

    end select
    
  end subroutine NUM_CONG_h

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

  subroutine NUM_CONG_Mach ( Rg, gama, u, T, M )

    real*8,dimension(:),intent(in)  :: Rg   ! constante do gás (J/kg.K)
    real*8,dimension(:),intent(in)  :: gama ! 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 ( gama * Rg * T )
  
  end subroutine NUM_CONG_Mach

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

  subroutine NUM_CONG_mi_k_congelados ( nxy, T, mi, k )

    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(inout) :: mi ! viscosidade absoluta (Pa.s)
    real*8,dimension(:),intent(inout) :: k  ! condutividade térmica (W/m.K)

    if ( tipo_mi == 2 .or. tipo_k  == 2 ) then
          
       do i = 1, nxy

          call GIBBS_MACH1D_recebe2 ( T(i) )

          call GIBBS_CONGELADO_calculos2

          call GIBBS_CONGELADO_envia2 ( mi(i), k(i) )

       end do

    end if

    if ( tipo_mi == 1 ) mi = mi_cte

    if ( tipo_k  == 1 ) k  = k_cte
   
  end subroutine NUM_CONG_mi_k_congelados

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

  subroutine NUM_CONG_pressao_entrada ( gama_cte, Rg_cte, p_cam, u_in, T_in, p_in )

    real*8,intent(in)  :: gama_cte ! razão de calores específicos (adimensional)
    real*8,intent(in)  :: Rg_cte   ! constante do gás (J/kg.K)
    real*8,intent(in)  :: p_cam    ! pressão na câmara de combustão (Pa)
    real*8,intent(in)  :: u_in     ! velocidade na entrada da tubeira (m/s)
    real*8,intent(in)  :: T_in     ! temperatura na entrada da tubeira (K)

    real*8,intent(out) :: p_in  ! pressão na entrada da tubeira (Pa)

    real*8 :: k    ! auxiliar
    real*8 :: M_in ! número de Mach na entrada

    M_in = u_in / dsqrt ( gama_cte * Rg_cte * T_in )

    k = 1.0d0 + 0.5d0 * ( gama_cte - 1.0d0 ) * ( M_in ** 2 )

    p_in = p_cam / ( k ** ( gama_cte / ( gama_cte - 1.0d0 ) ) )

  end subroutine NUM_CONG_pressao_entrada

!-------------------------------------------------
end module congelado