module usuario

! *** ESTE MÓDULO DEVE SER REVISADO A CADA NOVA SIMULAÇÃO
!     PARA EVENTUAIS ADAPTAÇÕES QUE SEJAM NECESSÁRIAS ***

! Objetivos: definir a estimativa inicial das incógnitas,
!            aplicar as condições de contorno,
!            definir funções para propriedades,
!            fazer pós-processamento.

! Versão atual: 27 Jan 09 (alterado em 25 Mai 11)

! Relação das sub-rotinas em ordem de apresentação:
!    1) calor_especifico
!    2) estimativa_inicial
!    3) ccu
!    4) ccv
!    5) ccT
!    6) ccp
!    7) Tin_pin_entrada
!    8) plin_entrada_p_ficticio
!    9) u_v_nos_ficticios_com_pl
!   10) Uce_Vcn_de_dn_faces_contornos
!   11) Uce_Vcn_faces_contornos_com_pl
!   12) contornos_pos_processamento
!   13) EXATO_MACH_fluxo_massa
!   14) EXATO_MACH_newton

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

use coeficientes

implicit none

real*8 :: pr = 1.01325d5   ! pressão atmosférica ao nível do mar

real*8 :: go = 9.80665d+0  ! aceleração gravitacional ao nível do mar

contains

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

  subroutine calor_especifico

    ! definição de propriedades variáveis ou constantes

    cp = gama * Rg / ( gama - 1 )

  end subroutine calor_especifico

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

  subroutine estimativa_inicial

    ! solução analítica Q1D para tubeiras circulares

    integer :: ind, npe, npw, npww
    real*8  :: dyf, dyr, ys,  xw,   dxf,  dxr,  aux 
    real*8  :: Si,  exp, cte, ue1D, Me1D, Te1D, pe1D

    ! determina o fluxo de massa da solução analítica 1D
    call EXATO_MACH_fluxo_massa ( gama, Rg, po, T0, Sg, fm1D )

    ! zeragem dos campos
    u    = 0.0d0
    v    = 0.0d0
    ue   = 0.0d0
    ve   = 0.0d0
    un   = 0.0d0
    vn   = 0.0d0
!    p    = 0.0d0
    p    = po
!    T    = 0.0d0
    T    = T0
    ro   = 0.0d0
    pl   = 0.0d0
    Uce  = 0.0d0
    Vcn  = 0.0d0
    de   = 0.0d0
    dn   = 0.0d0
    due  = 0.0d0
    dun  = 0.0d0
    dve  = 0.0d0
    dvn  = 0.0d0
    roe  = 0.0d0
    ron  = 0.0d0
    cup  = 0.0d0
    cvp  = 0.0d0
    plin = 0.0d0
    cp   = 0.0d0
    M    = 0.0d0

    ! estimativa inicial de u, p e T nos nós
    ! e cálculo da solução analítica 1D nos nós reais

    exp = gama / (gama - 1)

    do i = 2, nx-1

       j = ny-1
       np = (j-1)*nx + i

       Si = pi * ( rn(np) ** 2 ) 

       ind = 0
       if ( i > ig ) ind = 1
       call EXATO_MACH_newton ( gama, Si/Sg, ind, M1D(i) )

       cte = 1 + (gama - 1) * (M1D(i)**2) / 2

       p1D(i) = po / (cte ** exp)

       T1D(i) = T0 / cte

       u1D(i) = M1D(i) * dsqrt ( gama * Rg * T1D(i) )

       do j = 1, ny

          np = (j-1)*nx + i

          u(np) = u1D(i)
          p(np) = p1D(i)
          T(np) = T1D(i)
          if ( j > 1 .and. j < ny-1 ) Vcn(np) = v(np)*xk(np) - u(np)*yk(np)
          if ( j > 1 .and. j < ny-1 ) un(np) = u1D(i)
  
       end do

    end do

    ! calcula u e v nos nós fictícios de acordo com as condições de contorno
    call u_v_nos_ficticios_com_pl

    ! inicializa T e p na entrada da tubeira com a estimativa inicial
    call Tin_pin_entrada

    ! calcula p e T nos nós fictícios de acordo com as condições de contorno na entrada
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       p(np)  = 2 * pin(j) - p(npe)
       T(np)  = 2 * Tin(j) - T(npe)
    end do

    ! calcula p e T nos nós fictícios de acordo com as condições de contorno na saída
    i = nx
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npw = np - 1
       npww  = np - 2
       xw    = x(npw)
       dxf   = 2 * (xw - xp(npw))
       dxr   = xp(npw) - xp(npww)
       p(np) = p(npw) + (dxf/dxr) * (p(npw) - p(npww))
       T(np) = T(npw) + (dxf/dxr) * (T(npw) - T(npww))
    end do

!write(10,*) 'passou do ponto 1'

    ! estimativa inicial de u nas faces leste
    ! e cálculo da solução analítica 1D nos contornos

    do i = 1, nx-1

       j = ny-1
       np = (j-1)*nx + i

       Si = pi * ( raio(np) ** 2 ) 

       ind = 0
       if ( i > ig ) ind = 1
       call EXATO_MACH_newton ( gama, Si/Sg, ind, Me1D )

       cte = 1 + (gama - 1) * (Me1D**2) / 2

       pe1D = po / (cte ** exp)

       Te1D = T0 / cte

       ue1D = Me1D * dsqrt ( gama * Rg * Te1D )

       do j = 2, ny-1

          np = (j-1)*nx + i

          ue(np)  = ue1D
          Uce(np) = ue1D*ye(np) - v(np)*xe(np)
  
       end do

       ! contorno de entrada
       if ( i == 1 ) then
          u1D(i) = ue1D
          p1D(i) = pe1D
          T1D(i) = Te1D
          M1D(i) = Me1D
       end if

       ! contorno de saída
       if ( i == nx-1 ) then
          u1D(nx) = ue1D
          p1D(nx) = pe1D
          T1D(nx) = Te1D
          M1D(nx) = Me1D
          Fd1D  = fm1D * ue1D
          Fpv1D = pe1D * Si
       end if

    end do

!write(10,*) 'passou do ponto 2'

    ! inicializa as contravariantes nos contornos
    call Uce_Vcn_de_dn_faces_contornos

!write(10,*) 'passou do ponto 5'

    ! inicializa a massa específica
    call massa_especifica_nos

!write(10,*) 'passou do ponto 6'

    call massa_especifica_faces

!write(10,*) 'passou do ponto 3'

    ! inicializa os coeficientes e fontes
    a = 0.0d0
    b = 0.0d0
    do np = 1, nxy
       a(np,3) = 1.0d0
    end do

!write(10,*) 'passou do ponto 4'

    if ( ci == 2 ) then
       open(99,file=trim(adjustl(caso))//'.bin',form='unformatted')
       read(99) ito,norma1,tcpuo,pin,u,v,p,pl,ue,vn,T,ro,ve,un,Uce,Vcn,roe,ron
       itmax = ito + itmax
       ito = ito + 1
       close(99)
    end if

!write(10,*) 'passou estimativa_inicial'

  end subroutine estimativa_inicial

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

  subroutine ccu

    integer :: npw, npww, npe, npee
    real*8  :: xw,  dxf,  dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       npe  = np + 1
       npee = np + 2
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) =-1.0d0
       a(np,5) = 0.0d0
       dxf = 2 * xp(npe)
       dxr = xp(npee) - xp(npe)
       b(np)   = (dxf/dxr) * (u(npe) - u(npee))
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (u(npw) - u(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np = (j-1)*nx + i
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) =-1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do
  
  end subroutine ccu

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

  subroutine ccv

    integer :: npw, npww
    real*8  :: xw,  dxf,  dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 1.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (v(npw) - v(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np = (j-1)*nx + i
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do
  
  end subroutine ccv
!-------------------------------------------------

  subroutine ccT

    integer :: npw, npww
    real*8  :: xw,  dxf,  dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 1.0d0
       a(np,5) = 0.0d0
       b(np) = 2 * Tin(j)
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (T(npw) - T(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 1, nx
       np = (j-1)*nx + i
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) =-1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 1, nx
       np = (j-1)*nx + i
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do
  
  end subroutine ccT

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

  subroutine ccp

    integer :: npn, npnn, nps, npss, npw,  npww
    real*8  :: xw, dxf, dxr

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 1.0d0
       a(np,5) = 0.0d0
       b(np)   = 2 * plin(j)
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       a(np,1) = 0.0d0
       a(np,2) =-1.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       b(np)   = (dxf/dxr) * (pl(npw) - pl(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np   = (j-1)*nx + i
       a(np,1) = 0.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) =-1.0d0
       b(np)   = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np = (j-1)*nx + i
       a(np,1) =-1.0d0
       a(np,2) = 0.0d0
       a(np,3) = 1.0d0
       a(np,4) = 0.0d0
       a(np,5) = 0.0d0
       b(np)   = 0.0d0
    end do
  
  end subroutine ccp

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

  subroutine Tin_pin_entrada

    integer :: npe
    real*8  :: aux

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       uin(j) = ( u(np) + u(npe) ) / 2
       vin(j) = 0.0d0
       Tin(j) = T0 - ( gama - 1 ) * ( uin(j)**2 + vin(j)**2 ) / ( 2 * gama * Rg )
       Mw(j)  = dsqrt ( (uin(j)**2 + vin(j)**2) / (gama * Rg * Tin(j)) )
       aux    = 1 + ( gama - 1 ) * ( Mw(j) ** 2 ) / 2
       pin(j) = po / ( aux ** ( gama / ( gama - 1 ) ) )
    end do

  end subroutine Tin_pin_entrada

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

  subroutine plin_entrada_p_ficticio

    integer :: npe

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       p(np)   = 2 * pin(j) - p(npe)
       plin(j) = pin(j) - pina(j)
    end do

    end subroutine plin_entrada_p_ficticio

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

  subroutine u_v_nos_ficticios_com_pl

    integer :: npn, nps, npw, npww, npe, npee
    real*8  :: xw,  dxf, dxr 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       npe  = np + 1
       npee = np + 2
       dxf = 2 * xp(npe)
       dxr = xp(npee) - xp(npe)
       u(np) = u(npe) + (dxf/dxr) * (u(npe) - u(npee))
       v(np) = - v(npe)
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       npww = np - 2
       xw  = x(npw)
       dxf = 2 * (xw - xp(npw))
       dxr = xp(npw) - xp(npww)
       u(np) = u(npw) + (dxf/dxr) * (u(npw) - u(npww))
       v(np) = v(npw) + (dxf/dxr) * (v(npw) - v(npww))
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np  = (j-1)*nx + i
       npn = np + nx
       u(np) = u(npn)
       v(np) = -v(npn)
    end do

    ! contorno superior (parede)
    j = ny
    do i = 2, nx-1
       np  = (j-1)*nx + i
       nps = np - nx
       u(np) = u(nps)
       v(np) = v(nps)
    end do
  
  end subroutine u_v_nos_ficticios_com_pl

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

  subroutine Uce_Vcn_de_dn_faces_contornos

    integer :: npw, npe
    real*8  :: ue,  uin

    ! velocidades Uce e de

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npe = np + 1
       uin = ( u(np) + u(npe) ) / 2
       Uce(np) = uin * ye(np)
       de(np)  = de(npe)
    end do

    ! contorno direito (saída)
    i = nx-1
    do j = 2, ny-1
       np  = (j-1)*nx + i
       npw = np - 1
       npe = np + 1
       ue  = ( u(np) + u(npe) ) / 2
       Uce(np) = ue * ye(np)
       de(np)  = de(npw)
    end do

    ! velocidades Vcn e dn

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np = (j-1)*nx + i
       Vcn(np) = 0.0d0
       dn(np)  = 0.0d0
    end do

    ! contorno superior (parede)
    j = ny-1
    do i = 2, nx-1
       np  = (j-1)*nx + i
       Vcn(np) = 0.0d0
       dn(np)  = 0.0d0
    end do

  end subroutine Uce_Vcn_de_dn_faces_contornos

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

  subroutine Uce_Vcn_faces_contornos_com_pl

    ! velocidades Uce

    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       Uce(np) = Uce(np) + de(np) * (pl(np) - pl(np+1))
    end do

    i = nx-1
    do j = 2, ny-1
       np = (j-1)*nx + i
       Uce(np) = Uce(np) + de(np) * (pl(np) - pl(np+1))
    end do

  end subroutine Uce_Vcn_faces_contornos_com_pl

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

  subroutine contornos_pos_processamento

    integer :: npn, npw, nps, npnn
    real*8  :: raz 

    ! contorno esquerdo (entrada)
    i = 1
    do j = 2, ny-1
       np = (j-1)*nx + i
       nps = np - nx
       u(np) = uin(j)
       v(np) = vin(j)
       p(np) = pin(j)
       T(np) = Tin(j)
       xp(np) = x(np)
       yp(np) = ( y(nps) + y(np) ) / 2
    end do

    ! contorno direito (saída)
    i = nx
    do j = 2, ny-1
       np = (j-1)*nx + i
       npw  = np - 1
       u(np) = ( u(npw) + u(np) ) / 2
       v(np) = ( v(npw) + v(np) ) / 2
       p(np) = ( p(npw) + p(np) ) / 2
       T(np) = ( T(npw) + T(np) ) / 2
       xp(np) = x(npw)
       yp(np) = ( y(npw-nx) + y(npw) ) / 2
    end do

    ! contorno inferior (linha de simetria)
    j = 1
    do i = 2, nx-1
       np   = (j-1)*nx + i
       npw  = np - 1
       npn  = np + nx
       u(np) = u(npn)
       v(np) = 0.0d0
       p(np) = p(npn)
       T(np) = T(npn)
       xp(np) = ( x(npw) + x(np) ) / 2
       yp(np) = ( y(npw) + y(np) ) / 2
    end do

    ! canto esquerdo inferior
    j = 1
    i = 1
       np   = (j-1)*nx + i
       npn  = np + nx
       u(np) = u(npn)
       v(np) = 0.0d0
       p(np) = p(npn)
       T(np) = T(npn)
       xp(np) = x(np)
       yp(np) = y(np)

    ! canto direito inferior
    j = 1
    i = nx
       np   = (j-1)*nx + i
       npw  = np - 1
       npn  = np + nx
       u(np) = u(npn)
       v(np) = 0.0d0
       p(np) = p(npn)
       T(np) = T(npn)
       xp(np) = x(npw)
       yp(np) = y(npw)

    ! contorno superior (parede)
    j = ny
    do i = 1, nx
       np  = (j-1)*nx + i
       nps = np - nx
       u(np) = u(nps)
       v(np) = v(nps)
       p(np) = p(nps)
       T(np) = T(nps)
       xp(np) = ( x(nps) + x(nps-1) ) / 2
       yp(np) = ( y(nps) + y(nps-1) ) / 2
    end do

    ! cantos superiores
    j = ny
    do i = 1, nx, nx-1
       np  = (j-1)*nx + i
       nps = np - nx
       p(np) = p(nps)
       if ( i == 1 ) then
          xp(np) = x(nps)
          yp(np) = y(nps)
       end if
       if ( i == nx ) then
          xp(np) = x(nps-1)
          yp(np) = y(nps-1)
       end if
    end do
  
  end subroutine contornos_pos_processamento

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

  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_newton ( gama, a, ind, Mp )

    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 :: c, d, e, den, fat1, flinha ! variáveis auxiliares
    real*8 :: deltam ! estimativa do erro na solução de Mp

    integer :: nitm = 500        ! número máximo de iterações no método de Newton-Raphson
    real*8  :: tol  =   1.0d-15  ! tolerância no método de Newton-Raphson
    real*8  :: M_inicio_sub      = 1.0d-01 ! Mach inicial subsônico
    real*8  :: M_inicio_super    = 1.0d+01 ! Mach inicial supersônico
    real*8  :: M_inicio_negativo = 1.0d-02 ! Mach inicial se negativo
 
    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

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

end module usuario
