module EXATO_MACH

! sigla = EXATO_MACH_

! objetivo: calcular a solução analítica do escoamento (Mach)
!           isentrópico Q1D

! versão original 1.0 (21 Mai 99)

! programador: Carlos Henrique Marchi
! local:       SINMEC/UFSC - Florianópolis - SC

! versão atual    1.1 ( 9 Dez 99)

! última alteração =   14 Fev 00

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

! RELAÇÃO DAS ROTINAS DESTE MÓDULO:

  ! EXATO_MACH_abre_Newton_M
  ! EXATO_MACH_fluxo_massa
  ! EXATO_MACH_le_dados
  ! EXATO_MACH_le_inicio
  ! EXATO_MACH_Mach
  ! EXATO_MACH_Mach_u
  ! EXATO_MACH_massa_especifica
  ! EXATO_MACH_newton
  ! EXATO_MACH_p_tot
  ! EXATO_MACH_passa_dados
  ! EXATO_MACH_pressao
  ! EXATO_MACH_T_tot
  ! EXATO_MACH_temperatura
  ! EXATO_MACH_u_Mach

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

! rotinas onde são lidos ARQUIVOS de dados neste MÓDULO:

  ! EXATO_MACH_le_dados
  ! EXATO_MACH_le_inicio

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

! rotinas onde são gerados ARQUIVOS neste MÓDULO:

  ! EXATO_MACH_Mach

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

use ARQUIVO ! versão 1.2 (19 Jan 00)

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

implicit none

real*8,private  :: gama  ! razão de calores específicos (adimensional)
real*8,private  :: Rg    ! constante do gás (J/kg.K)
real*8,private  :: p_cam ! pressão na câmara de combustão (Pa)
real*8,private  :: T_cam ! temperatura na câmara de combustão (K)
real*8,private  :: g     ! aceleração gravitacional (m/s2)
real*8,private  :: tol   ! tolerância no método de Newton-Raphson

integer,private :: nitm  ! número máximo de iterações no método de Newton-Raphson
integer,private :: w_it  ! freqüência de escrita das iterações de M; 0=não escreve

real*8, private :: M_inicio_sub      ! Mach inicial subsônico
real*8, private :: M_inicio_super    ! Mach inicial supersônico
real*8, private :: M_inicio_negativo ! Mach inicial se negativo

integer,private :: lista_itera ! lista conteúdo da solução iterativa do Mach
integer,private :: lista       ! lista conteúdo do arquivo de dados do módulo EXATO_MACH
integer,private :: le          ! lê conteúdo do arquivo de dados do módulo EXATO_MACH
integer,private :: copia       ! gera back up do arquivo de dados do módulo EXATO_MACH

character*50,private :: nome_exato_mach ! nome do arquivo de dados do módulo EXATO_MACH

logical,private :: dos ! acessa prompt dos

character*220,private :: comando ! linha de comando

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

contains

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

  subroutine EXATO_MACH_abre_Newton_M

    character*(150) :: trabalho ! nome do diretório para saídas
    character*(50)  :: caso     ! nome do caso

    if ( w_it == 0 .or. lista_itera ==0 ) return

    call ARQUIVO_passa_diretorio_trabalho ( trabalho )
    call ARQUIVO_passa_caso ( caso )
    comando = 'notepad '// trim(trabalho) // &
              trim(caso) // '.Mach1D.Newton_M.txt'
    dos = systemqq ( comando )

  end subroutine EXATO_MACH_abre_Newton_M

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

  subroutine EXATO_MACH_fluxo_massa ( gama, Rg, p_cam, T_cam, Sg, fm_exato )

    real*8,intent(in)  :: gama     ! razão de calores específicos (adimensional)
    real*8,intent(in)  :: Rg       ! 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)  :: T_cam    ! temperatura na câmara de combustão (K)
    real*8,intent(in)  :: Sg       ! área da garganta (m2)
    real*8,intent(out) :: fm_exato ! fluxo de massa exato (kg/s)

    real*8 :: a1, a2, a3 ! auxiliares

    a1 = gama / ( Rg * T_cam )

    a2 = 2.0d0 / ( gama + 1.0d0 )

    a3 = (gama + 1.0d0 ) / ( gama - 1.0d0 )

    fm_exato = p_cam * Sg * dsqrt ( a1 * (a2**a3) )

  end subroutine EXATO_MACH_fluxo_massa

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

  subroutine EXATO_MACH_le_dados ( unit )

    integer,intent(in) :: unit ! unidade de leitura do arquivo de dados

    open ( unit, file = trim(nome_exato_mach) )

    read ( unit, * ) gama
    read ( unit, * ) Rg
    read ( unit, * ) p_cam
    read ( unit, * ) T_cam
    read ( unit, * ) nitm
    read ( unit, * ) tol
    read ( unit, * ) w_it
    read ( unit, * ) M_inicio_sub
    read ( unit, * ) M_inicio_super
    read ( unit, * ) M_inicio_negativo
    read ( unit, * ) lista_itera
    read ( unit, * ) g

    close ( unit )

  end subroutine EXATO_MACH_le_dados 

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

  subroutine EXATO_MACH_le_inicio ( unit, nome )

    integer,      intent(in) :: unit ! unidade do arquivo
    character*(*),intent(in) :: nome ! nome do arquivo de dados de configuração do módulo

    open ( unit, file = trim(nome) )

    read ( unit, * ) nome_exato_mach 
    read ( unit, * ) lista
    read ( unit, * ) le
    read ( unit, * ) copia

    close ( unit )

    if ( lista == 1 ) call ARQUIVO_notepad_abre_dado ( nome_exato_mach )

    if ( le    == 1 ) call EXATO_MACH_le_dados ( unit )

    if ( copia == 1 ) call ARQUIVO_copia_arquivo ( nome_exato_mach )

  end subroutine EXATO_MACH_le_inicio

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

  subroutine EXATO_MACH_Mach ( n, ng, gama, a, M )

    integer,intent(in) :: n    ! número de pontos a calcular
    integer,intent(in) :: ng   ! índice da garganta
    real*8, intent(in) :: gama ! razão de calores específicos (adimensional)
    real*8, intent(in), dimension(:) :: a ! razão de áreas (adimensional)
    real*8, intent(out),dimension(:) :: M ! número de Mach (adimensional)

    integer :: i             ! número do ponto
    integer :: ind           ! 0=escoamento subsônico; 1=supersônico
    integer :: n_it          ! número de iterações usadas para obter M
    integer :: unit_rot = 71 ! identificador do arquivo

    real*8 :: erro_M ! estimativa do erro na solução de M

    character*(50) :: caso ! nome do caso

    if ( w_it /= 0 ) then
       call ARQUIVO_passa_caso ( caso )
       comando = trim(caso) // '.Mach1D.Newton_M.txt'
       call ARQUIVO_cria_arquivo ( unit_rot, comando )
    end if

    do i = 1, n

       if ( i <= ng ) then
          ind = 0
       else
          ind = 1
       end if
 
       call EXATO_MACH_newton ( gama, a(i), ind, M(i), n_it, erro_M )

       if ( w_it /= 0 ) then
          if ( i==1 .or. i==n .or. mod(i,w_it)==0 ) &
             write(unit_rot,1) i, M(i), n_it, erro_M
             1 format('i=',i6,5x,'M=',1pe10.3,5x,'iterações=',i4,5x,'erro=',1pe10.3)
       end if

    end do

    if ( w_it /= 0 ) close(unit_rot)

  end subroutine EXATO_MACH_Mach

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

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

    real*8, intent(in) :: gama ! razão de calores específicos (adimensional)
    real*8, intent(in) :: Rg   ! constante do gás (J/kg.K)
    real*8, intent(in), dimension(:) :: u ! velocidade (m/s)
    real*8, intent(in), dimension(:) :: T ! temperatura (K)
    real*8, intent(out),dimension(:) :: M ! número de Mach (adimensional)

    M = u / dsqrt ( gama * Rg * T )
  
  end subroutine EXATO_MACH_Mach_u

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

  subroutine EXATO_MACH_massa_especifica ( Rg, p, T, ro )

    real*8,intent(in) :: Rg ! constante do gás (J/kg.K)

    real*8,dimension(:),intent(in)  :: p  ! pressão (Pa)
    real*8,dimension(:),intent(in)  :: T  ! temperatura (K)

    real*8,dimension(:),intent(out) :: ro ! massa específica (kg/m3)

    ro = p / ( Rg * T )

  end subroutine EXATO_MACH_massa_especifica

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

  subroutine EXATO_MACH_newton ( gama, a, ind, Mp, j, deltam )

    real*8,intent(in)   :: a      ! razão de áreas (adimensional)
    real*8,intent(in)   :: gama   ! razão de calores específicos (adimensional)
    integer,intent(in)  :: ind    ! 0=escoamento subsônico; 1=supersônico
    real*8,intent(out)  :: Mp     ! solução do número de Mach analítico (adimensional)
    real*8,intent(out)  :: deltam ! estimativa do erro na solução de Mp
    integer,intent(out) :: j      ! número da iteração
    real*8 :: c, d, e, den, fat1, flinha ! variáveis auxiliares
 
    if ( ind == 1 ) Mp = M_inicio_super
    if ( ind == 0 ) Mp = M_inicio_sub

    e = (gama - 1.0d0) / 2.0d0
    c = (gama + 1.0d0) / ( 2.0d0 - 2.0D0 * gama)
    d = ((gama + 1.0d0) / 2.0d0) ** c

    do j = 1, nitm
      fat1 = Mp * ( c * (1.0d0+e*Mp*Mp) ** (c - 1.0d0)    &
             * 2.0d0 * e * Mp)
      den = Mp * ( 1.0d0 + e * Mp * Mp) ** c
      flinha = - d * ((1.0d0+e*Mp*Mp) ** c + fat1) / (den * den)
      deltam = - (d/den-a) / flinha
      Mp = Mp + deltam
      if ( Mp < 0.0d0 ) Mp = M_inicio_negativo
      if ( dabs(deltam) < tol ) exit
    end do

  end subroutine EXATO_MACH_newton

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

  subroutine EXATO_MACH_p_tot ( n, gama, M, p, p_tot )

    integer,intent(in) :: n    ! número de pontos a calcular
    real*8, intent(in) :: gama ! razão de calores específicos (adimensional)
    real*8, intent(in), dimension(:) :: M     ! número de Mach (adimensional)
    real*8, intent(in), dimension(:) :: p     ! pressão (Pa)
    real*8, intent(out),dimension(:) :: p_tot ! pressão de estagnação (Pa)

    integer :: i  ! número do ponto

    real*8 :: cte ! variável auxiliar
    real*8 :: exp ! variável auxiliar

    do i = 1, n

       exp = gama / (gama - 1.0d0)

       cte = 1.0d0 + (gama - 1.0d0) * (M(i)**2) / 2.0d0
  
       p_tot(i) = p(i) * (cte ** exp)

    end do

  end subroutine EXATO_MACH_p_tot

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

  subroutine EXATO_MACH_passa_dados ( gama_r, Rg_r, p_cam_r, T_cam_r, g_r )

    real*8, intent(out) :: gama_r  ! razão de calores específicos (adimensional)
    real*8, intent(out) :: Rg_r    ! constante do gás (J/kg.K)
    real*8, intent(out) :: p_cam_r ! pressão na câmara de combustão (Pa)
    real*8, intent(out) :: T_cam_r ! temperatura na câmara de combustão (K)
    real*8, intent(out) :: g_r     ! aceleração gravitacional (m/s2)

    gama_r = gama
    Rg_r   = Rg
    p_cam_r = p_cam
    T_cam_r = T_cam
    g_r = g

  end subroutine EXATO_MACH_passa_dados

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

  subroutine EXATO_MACH_pressao ( n, gama, p_cam, M, p )

    integer,intent(in) :: n     ! número de pontos a calcular
    real*8, intent(in) :: gama  ! razão de calores específicos (adimensional)
    real*8, intent(in) :: p_cam ! pressão na câmara de combustão (Pa)
    real*8, intent(in), dimension(:) :: M ! número de Mach (adimensional)
    real*8, intent(out),dimension(:) :: p ! pressão (Pa)

    integer :: i  ! número do ponto

    real*8 :: cte ! variável auxiliar
    real*8 :: exp ! variável auxiliar

    do i = 1, n

       exp = gama / (gama - 1.0d0)

       cte = 1.0d0 + (gama - 1.0d0) * (M(i)**2) / 2.0d0
  
       p(i) = p_cam / (cte ** exp)

    end do

  end subroutine EXATO_MACH_pressao

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

  subroutine EXATO_MACH_T_tot ( n, gama, M, T, T_tot )

    integer,intent(in) :: n    ! número de pontos a calcular
    real*8, intent(in) :: gama ! razão de calores específicos (adimensional)
    real*8, intent(in), dimension(:) :: M     ! número de Mach (adimensional)
    real*8, intent(in), dimension(:) :: T     ! temperatura (K)
    real*8, intent(out),dimension(:) :: T_tot ! temperatura de estagnação (K)

    integer :: i  ! número do ponto

    real*8 :: cte ! variável auxiliar

    do i = 1, n

       cte = 1.0d0 + (gama - 1.0d0) * (M(i)**2) / 2.0d0
  
       T_tot(i) = T(i) * cte
  
    end do

  end subroutine EXATO_MACH_T_tot

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

  subroutine EXATO_MACH_temperatura ( n, gama, T_cam, M, T )

    integer,intent(in) :: n     ! número de pontos a calcular
    real*8, intent(in) :: gama  ! razão de calores específicos (adimensional)
    real*8, intent(in) :: T_cam ! temperatura na câmara de combustão (K)
    real*8, intent(in), dimension(:) :: M ! número de Mach (adimensional)
    real*8, intent(out),dimension(:) :: T ! temperatura (K)

    integer :: i  ! número do ponto

    real*8 :: cte ! variável auxiliar

    do i = 1, n

       cte = 1.0d0 + (gama - 1.0d0) * (M(i)**2) / 2.0d0
  
       T(i) = T_cam / cte
  
    end do

  end subroutine EXATO_MACH_temperatura

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

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

    real*8, intent(in) :: gama ! razão de calores específicos (adimensional)
    real*8, intent(in) :: Rg   ! constante do gás (J/kg.K)
    real*8, intent(in), dimension(:) :: M ! número de Mach (adimensional)
    real*8, intent(in), dimension(:) :: T ! temperatura (K)
    real*8, intent(out),dimension(:) :: u ! velocidade (m/s)

    u = M * dsqrt ( gama * Rg * T )

  end subroutine EXATO_MACH_u_Mach

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

end module EXATO_MACH