module malha

! *** ESTE MÓDULO DEVE SER ADAPTADO PARA CADA NOVA GEOMETRIA ***

! Objetivos: gerar, escrever e visualizar o domínio de cálculo,
!            gerar, escrever e visualizar a malha.

! Versão atual: 15 Mar 09

! Relação das sub-rotinas em ordem de apresentação:
!   1) pre_processamento
!   2) define_dominio
!   3) gera_malha_UL
!   4) gera_malha_PG
!   5) PG_malha_NS
!   6) parametros_malha
!   7) parametros_tubeira
!   8) escreve_dominio
!   9) grafico_dominio
!  10) escreve_malha
!  11) grafico_malha

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

! Descrição do domínio:

!   Tipo:       tubeira câmara-cossenóide
!   Origem:     Marchi et al. (2000)
!   RA-conv.:   1.1025
!   RA-div.:    1.1025
!   Malha-base: 5 volumes na direção axial
!               2 volumes na direção radial

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

use DADOS

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

implicit none

real*8 :: rcg ! raio de curvatura na garganta da tubeira

contains

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

  subroutine pre_processamento

    call define_dominio

    select case ( grid )
       case ( 1 ) ! uniforme por linha (U-L)
          call gera_malha_UL
       case ( 2 ) ! PG
          call gera_malha_PG
    end select

    call parametros_malha

    call parametros_tubeira

    call escreve_dominio

    call grafico_dominio

    call escreve_malha

    call grafico_malha

  end subroutine pre_processamento

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

  subroutine define_dominio
 
    integer   NXc,NXt,NPS
    real*8    dx,dy
    real*8 :: Lc ! comprimento da câmara (m)
    real*8 :: Ln ! comprimento da tubeira (m)
    real*8 :: L  ! comprimento total: câmara + tubeira (m)
    real*8 :: Ri ! raio na entrada da tubeira (m)
    real*8 :: Rg ! raio na garganta da tubeira (m)
    real*8 :: d  ! amplitude da cossenóide (m)
 
    ! *** inicialização ***

    x = 0.0d0
    y = 0.0d0

    NXc = (nx-2) / 5
    NXt = nx - 2 - Nxc 

    ! *** dados ***

    Lc = 0.100d0
    Ln = 0.400d0
    
    select case ( razao_ig )
    case (1) ! Razão de refino: 9
        Ri = 0.120d0  ! para razão 9: 0.120d0;  para razão 1.1025: 0.105d0; para razão 4: 0.120d0; para razão 1.5625: 0.05
        Rg = 0.040d0  ! para razão 9: 0.040d0;  para razão 1.1025: 0.100d0; para razão 4: 0.060d0; para razão 1.5625: 0.04
    case (2) ! Razão de refino: 1.1025
        Ri = 0.105d0
        Rg = 0.100d0
    end select

    L = Lc + Ln
    d = ( ri - rg ) / 2
 
    ! *** CONTORNO SUL ***

    ! câmara
    J = 1
    I = 1
    np = (j-1)*nx + i
    DX = Lc / Nxc
    X(np) = 0.0d0
    Y(np) = 0.0d0
    DO I = 2, NXc
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    I = NXc + 1
    np = (j-1)*nx + i
    X(np) = Lc
    Y(np) = 0.0d0
 
    ! tubeira cossenoidal
    DX = Ln / NXt
    DO I = NXc+2, NX-2
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    I = NX-1
    np = (j-1)*nx + i
    X(np) = L
    Y(np) = 0.0d0
 
    ! *** CONTORNO NORTE ***
    
    ! câmara
	J = ny-1
    DO I = 2, NXc
       np  = (j-1)*nx + i
       nps = I
       X(np) = X(nps)
       Y(np) = Ri
    end do
    I = NXc + 1
    np = (j-1)*nx + i
    X(np) = Lc
    Y(np) = Ri
 
    ! índice i da garganta
    ig = Nxc * 3 + 1

    ! tubeira cossenoidal
    DO I = NXc+2, NX-2
       np  = (j-1)*nx + i
       nps = I
       X(np) = X(nps)
       Y(np) = Rg + d * ( 1 + dcos ( 2*pi*(x(np)-Lc) / Ln ) )
    end do
    I = NX-1
    np = (j-1)*nx + i
    X(np) = L
    Y(np) = Ri
 
    ! *** CONTORNO OESTE ***

    DY = Ri / (ny-2)
    I = 1
    DO J = 2, ny-2
       np = (j-1)*nx + i
       X(np) = 0.0d0
       Y(np) = Y(np-nx) + DY
    end do
    j = ny-1
    np = (j-1)*nx + i
    X(np) = 0.0d0
    Y(np) = Ri
 
    ! *** CONTORNO LESTE ***
    
	DY = Ri / (ny-2)
    I = nx-1
    DO J = 2, ny-2
       np = (j-1)*nx + i
       X(np) = L
       Y(np) = Y(np-nx) + DY
    end do
    j = ny-1
    np = (j-1)*nx + i
    X(np) = L
    Y(np) = Ri

    ! *** raio de curvatura na garganta ***

    rcg = ( Ln ** 2 ) / ( 2 * (pi**2) * ( Ri - Rg ) )
 
  end subroutine define_dominio

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

  subroutine gera_malha_UL

    integer :: npn, nps
    real*8  :: dy

    ! *** PONTOS INTERNOS ***
    
	DO I = 2, nx-2

       npn = (ny-2)*nx + i

       nps = i

       DY = (Y(npn) - Y(nps)) / (ny-2)

       DO J = 2, ny-2

          np = (j-1)*nx + i

          X(np) = X(nps)

          Y(np) = Y(np-nx) + DY

       end do  

    end do

  end subroutine gera_malha_UL

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

  subroutine gera_malha_PG

    integer :: npn, nps
    real*8  :: Ly, qpg, dy

    ! *** PONTOS INTERNOS ***
    
	DO I = 1, nx-1

       npn = (ny-2)*nx + i

       nps = i

       Ly = Y(npn) - Y(nps)

       call PG_malha_NS ( A1, Ly, Ny-2, qpg )

       DO J = ny-2, 2, -1

          np = (j-1)*nx + i

          X(np) = X(nps)

          dy = a1 * ( qpg ** (ny-2-j) )

          Y(np) = Y(np+nx) - DY

       end do  

    end do

  end subroutine gera_malha_PG

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

      SUBROUTINE PG_malha_NS ( A1, SN, N, Q )
     
      real*8  :: tol, qa, sa, a1, sn, q, ra, fq, flq, dq
      integer :: N, IM

      TOL = 1.0d-15
      IM = 500
      QA = 2.0d0
      SA = SN/A1
      RA = 1.0d0*N
      IF(SA.LT.RA) QA = 0.9d0

      DO k = 1,IM
         FQ = A1 * (QA**N - 1)/SN + 1 - QA
         FLQ = (A1 * N * QA**(N-1))/SN - 1
         Q = QA - FQ/FLQ
         DQ = dABS(Q-QA)
         IF(DQ.LT.TOL) exit
         QA = Q
      end do

      END subroutine PG_malha_NS

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

  subroutine parametros_malha

    real*8  :: dy, dx, dxmax, dxmin, dymax, dymin

    i  = nx-1
    j  = ny-1
    np = (j-1)*nx + i
    dxmin = x(np)
    dxmax = 0
    dymin = y(np)
    dymax = 0
    
	DO j = 1, ny-1

       DO i = 2, nx-1

          np = (j-1)*nx + i

          dx = dabs ( x(np) - x(np-1) )

          if ( dx < dxmin ) dxmin = dx

          if ( dx > dxmax ) dxmax = dx

       end do  

    end do

	DO i = 1, nx-1

       DO j = 2, ny-1

          np = (j-1)*nx + i

          dy = dabs ( y(np) - y(np-nx) )

          if ( dy < dymin ) dymin = dy

          if ( dy > dymax ) dymax = dy

       end do  

    end do

    write(10,1) dxmin, dxmax, dxmax/dxmin, dymin, dymax, dymax/dymin

    1 format(//,'*** Parâmetros da malha ***', //, &
             1pe22.9,' = dxmin: tamanho mínimo dos volumes de controle na direção x (m)',/, &
             1pe22.9,' = dxmax: tamanho máximo dos volumes de controle na direção x (m)',/, &
             1pe22.9,' = dxmax / dxmin',/, &
             1pe22.9,' = dymin: tamanho mínimo dos volumes de controle na direção y (m)',/, &
             1pe22.9,' = dymax: tamanho máximo dos volumes de controle na direção y (m)',/, &
             1pe22.9,' = dymax / dymin')

  end subroutine parametros_malha

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

  subroutine parametros_tubeira

    integer :: npg,  npin, npex, npw, npe
    real*8  :: rain, raex

    i    = ig
    j    = ny-1
    npg  = (j-1)*nx + i
    npw  = npg - 1
    npe  = npg + 1

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

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

    rcag = rcg / y(npg)

    Sg   = pi * ( y(npg) ** 2 )

    rain = ( y(npin) / y(npg) ) ** 2

    raex = ( y(npex) / y(npg) ) ** 2

    write(10,1) ig, x(npg), y(npg), rcg, rcag, Sg, rain, raex

    1 format(//,'*** Parâmetros da Tubeira ***', //, &
              i8,12x,'   = ig:   número do volume na direção X cuja face leste', &
                           ' coincide com a garganta da tubeira',/, &
             1pe22.9,' = Xg:   coordenada X da garganta da tubeira (m)',/, &
             1pe22.9,' = rg:   raio da garganta da tubeira (m)',/, &
             1pe22.9,' = rcg:  raio de curvatura na garganta da tubeira (m)',/, &
             1pe22.9,' = Rcag: raio de curvatura adimensional na garganta da tubeira',/, &
             1pe22.9,' = Sg:   área da garganta da tubeira (m2)',/, &
             1pe22.9,' = rain: razão de áreas do convergente',/, &
             1pe22.9,' = raex: razão de áreas do divergente')

  end subroutine parametros_tubeira

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

  subroutine escreve_dominio

    ! escreve contorno sul

    open(7,file=trim(adjustl(caso))//'.contorno_sul.dat')

    write(7,1) 

    j = 1
    do i = 1, nx-1
       np = (j-1)*nx + i
       if ( i==1 .or. i==nx-1 .or. mod(i,w_g)==0 ) &
          write(7,2) i, j, np, x(np), y(np)
    end do

    close(7)

    ! escreve contorno norte

    open(7,file=trim(adjustl(caso))//'.contorno_norte.dat')

    write(7,1) 

    j = ny-1
    do i = 1, nx-1
       np = (j-1)*nx + i
       if ( i==1 .or. i==nx-1 .or. mod(i,w_g)==0 ) &
          write(7,2) i, j, np, x(np), y(np)
    end do

    close(7)

    ! escreve contorno oeste

    open(7,file=trim(adjustl(caso))//'.contorno_oeste.dat')

    write(7,1) 

    i = 1
    do j = 1, ny-1
       np = (j-1)*nx + i
       if ( j==1 .or. j==ny-1 .or. mod(j,w_g)==0 ) &
          write(7,2) i, j, np, x(np), y(np)
    end do

    close(7)

    ! escreve contorno leste

    open(7,file=trim(adjustl(caso))//'.contorno_leste.dat')

    write(7,1) 

    i = nx-1
    do j = 1, ny-1
       np = (j-1)*nx + i
       if ( j==1 .or. j==ny-1 .or. mod(j,w_g)==0 ) &
          write(7,2) i, j, np, x(np), y(np)
    end do

    close(7)

    1 format('#', t5,'i', t10,'j', t19,'np', t24,'x (m)', t48,'y (m)')

    2 format(2(i5),i10,2(1pe24.15))

  end subroutine escreve_dominio

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

  subroutine grafico_dominio

    nome = 'Mach2D_6p1_dominio.gnu'

    open(7,file=trim(adjustl(nome)))

    do k = 1, 6
       read(7,*)
    end do

    write(7,1) trim(adjustl(caso)) // '.contorno_sul.dat'
    1 format("plot '", a, "' using 4:5 title 'sul'")

    write(7,2) trim(adjustl(caso)) // '.contorno_norte.dat'
    2 format("replot '", a, "' using 4:5 title 'norte'")

    write(7,3) trim(adjustl(caso)) // '.contorno_oeste.dat'
    3 format("replot '", a, "' using 4:5 title 'oeste'")

    write(7,4) trim(adjustl(caso)) // '.contorno_leste.dat'
    4 format("replot '", a, "' using 4:5 title 'leste'")

    write(7,5) trim(adjustl(title))
    5 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot ' // trim(adjustl(nome))

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_dominio

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

  subroutine escreve_malha

    ! cria o arquivo da malha

    open(7,file=trim(adjustl(caso))//'.malha.dat')

    ! escreve linhas entre os contornos oeste e leste

    do j = 1, ny-1

       write(7,1) 

       do i = 1, nx-1
          np = (j-1)*nx + i
          if ( i==1 .or. i==nx-1 .or. mod(i,w_g)==0 ) &
             write(7,2) x(np), y(np)
       end do
   
    end do

    ! escreve linhas entre os contornos sul e norte

    do i = 1, nx-1

       write(7,1) 

       do j = 1, ny-1
          np = (j-1)*nx + i
          if ( j==1 .or. j==ny-1 .or. mod(j,w_g)==0 ) &
             write(7,2) x(np), y(np)
       end do
   
    end do

    close(7)

    1 format(x)

    2 format(2(1pe24.15))

  end subroutine escreve_malha

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

  subroutine grafico_malha

    nome = 'Mach2D_6p1_malha.gnu'

    open(7,file=trim(adjustl(nome)))

    do k = 1, 5
       read(7,*)
    end do

    write(7,1) trim(adjustl(caso)) // '.malha.dat'
    1 format("plot '", a, "' notitle")

    write(7,2) trim(adjustl(title))
    2 format("set title '", a, "'", /, "replot")

    close(7)

    nome = 'wgnuplot ' // trim(adjustl(nome))

    if ( sem_g == 0 ) dos = system(trim(adjustl(nome)))

  end subroutine grafico_malha

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

end module malha
