module malha

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

! Versão atual:  10 Out 06

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


use DADOS

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

implicit none

contains

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

  subroutine define_dominio

	if (geom==1.0) call BACK_1965
	if (geom==2.0) call cossenoide

  end subroutine define_dominio
 

!-------------------------------------------------
  
  subroutine BACK_1965

! Descrição do domínio:

!   Tipo:       tubeira convergente-divergente
!   Origem:     artigo de L. H. Back, P. F. Massier e H. L. Gier,
!                  AIAA Journal, vol. 3, p. 1606-1614, 1965. 
!   RA-conv.:   9.76
!   RA-div.:    6.63
!   Malha-base: 90 volumes na direção axial

    integer fator ! razão de refino da malha-base; múltiplo de 90 volumes reais em X
    integer nvk1,nvk2,nvk3,nvk4,nvk5,nvk6,nps
    integer nk1,nk2,nk3,nk4,nk5,nk6,nk7,nk8,nk9,nk10,nk11,nk12,nk13,nk14,nk15,nk16,nk17
    real*8  L1,L2,L3,L4,L5,L6,dx,dy
    real*8  ri,rt,re,rl1,rl2,rl3,rl5,alfa,beta
 
    ! *** inicialização ***

    x = 0.0d0
    y = 0.0d0

    fator = (nx-2)/90

    NVK1 =  4 * fator
    NVK2 =  7 * fator
    NVK3 = 17 * fator
    NVK4 =  5 * fator
    NVK5 =  2 * fator
    NVK6 = 55 * fator

    L1 = 0.007874d0
    L2 = 0.022250d0
    L3 = 0.055880d0
    L4 = 0.064872d0
    L5 = 0.068148d0
    L6 = 0.185039d0

    RI = 0.063482d0
    RT = 0.020320d0
    RE = 0.052322d0
    RC = 0.012700d0
 
    RL1 = 0.020320d0
    RL2 = 0.057523d0
    RL3 = 0.024051d0
    RL5 = 0.020750d0
 
    ALFA = 44.865
    BETA = 15.115
    ALFA = ALFA * PI / 180.0d0
    BETA = BETA * PI / 180.0d0
 
    ! *** CONTORNO SUL ***

    J = 1
    I = 1
    np = (j-1)*nx + i
    DX = L1 / NVK1
    X(np) = 0.0d0
    Y(np) = 0.0d0
    NK1 = NVK1
    DO I = 2, NK1
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    NK2 = NK1 + 1
    i = nk2
    np = (j-1)*nx + i
    X(np) = L1
    Y(np) = 0.0d0
 
    DX = (L2 - L1) / NVK2
    NK3 = NK2 + 1
    NK4 = NK2 + NVK2 - 1
    DO I = NK3, NK4
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    NK5 = NK4 + 1
    i = nk5
    np = (j-1)*nx + i
    X(np) = L2
    Y(np) = 0.0d0
 
    DX = (L3 - L2) / NVK3
    NK6 = NK5 + 1
    NK7 = NK5 + NVK3 - 1
    DO I = NK6, NK7
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    NK8 = NK7 + 1
    i = nk8
    np = (j-1)*nx + i
    X(np) = L3
    Y(np) = 0.0d0
 
    DX = (L4 - L3)/ NVK4
    NK9 = NK8 + 1
    NK10 = NK8 + NVK4 - 1
    DO I = NK9, NK10
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    NK11 = NK10 + 1
    i = nk11
    np = (j-1)*nx + i
    X(np) = L4
    Y(np) = 0.0d0
 
    DX = (L5 - L4)/ NVK5
    NK12 = NK11 + 1
    NK13 = NK11 + NVK5 - 1
    DO I = NK12, NK13
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    NK14 = NK13 + 1
    i = nk14
    np = (j-1)*nx + i
    X(np) = L5
    Y(np) = 0.0d0
 
    DX = (L6 - L5) / NVK6
    NK15 = NK14 + 1
    NK16 = NK14 + NVK6 - 1
    DO I = NK15, NK16
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    NK17 = NK16 + 1
    i = nk17
    np = (j-1)*nx + i
    X(np) = L6
    Y(np) = 0.0d0
 
    ! *** 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 = RE / (ny-2)
    I = nx-1
    DO J = 2, ny-2
       np = (j-1)*nx + i
       X(np) = L6
       Y(np) = Y(np-nx) + DY
    end do
    j = ny-1
    np = (j-1)*nx + i
    X(np) = L6
    Y(np) = RE
 
    ! *** CONTORNO NORTE ***
    
	J = ny-1
    NK1 = NVK1
    DO I = 2, NK1
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       Y(np) = RI
    end do
    NK2 = NK1 + 1
    i = nk2
    np = (j-1)*nx + i
    X(np) = L1
    Y(np) = RI
 
    NK3 = NK2 + 1
    NK4 = NK2 + NVK2 - 1
    DO I = NK3, NK4
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       DY = dSQRT( RL1**2 - (X(np)-L1)**2 )
       Y(np) = RI - RL1 + DY
    end do
    NK5 = NK4 + 1
    i = nk5
    np = (j-1)*nx + i
    X(np) = L2
    Y(np) = RL2
 
    NK6 = NK5 + 1
    NK7 = NK5 + NVK3 - 1
    DO I = NK6, NK7
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       DY = (X(np) - L2) * dTAN(ALFA)
       Y(np) = RL2 - DY
    end do
    NK8 = NK7 + 1
    i = nk8
    np = (j-1)*nx + i
    X(np) = L3
    Y(np) = RL3
 
    NK9 = NK8 + 1
    NK10 = NK8 + NVK4 - 1
    DO I = NK9, NK10
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       DY = dSQRT( RC**2 - (X(np)-L4)**2 )
       Y(np) = RT + RC - DY
    end do
    NK11 = NK10 + 1
    i = nk11
    np = (j-1)*nx + i
    X(np) = L4
    Y(np) = RT

    ! índice i da garganta

    ig = i
 
    NK12 = NK11 + 1
    NK13 = NK11 + NVK5 - 1
    DO I = NK12, NK13
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       DY = dSQRT( RC**2 - (X(np)-L4)**2 )
       Y(np) = RT + RC - DY
    end do
    NK14 = NK13 + 1
    i = nk14
    np = (j-1)*nx + i
    X(np) = L5
    Y(np) = RL5
 
    NK15 = NK14 + 1
    NK16 = NK14 + NVK6 - 1
    DO I = NK15, NK16
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       DY = (L6 - X(np)) * dTAN(BETA)
       Y(np) = RE - DY
    end do
    NK17 = NK16 + 1
    i = nk17
    np = (j-1)*nx + i
    X(np) = L6
    Y(np) = RE

  end subroutine BACK_1965

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

  subroutine cossenoide

! Descrição do domínio:

!   Tipo:       tubeira convergente-divergente
!   Origem:     artigo de C. H. Marchi,  F. Laroca, A. F. C. Silva e J. N. Hinckel,
!                  Numerical Heat Transfer, part A, vol. 45, p. 699-717, 2004. 
!   RA-conv.:   9.00
!   RA-div.:    9.00
!   Malha-base: 5 volumes na direção axial

    integer fator ! razão de refino da malha-base; múltiplo de 5 volumes reais em X
	integer :: nvk1, nvk2, nps
	integer :: nk1, nk2, nk3, nk4, nk5
	real*8 :: L1, L2, dx, dy
	real*8 :: ri, rg, rc
	real*8 :: pi
  
    ! *** inicialização ***

	x = 0.0d0
    y = 0.0d0

	fator = (nx-2)/5

    NVK1 =  1 * fator
    NVK2 =  4 * fator 

	L1 = 0.1d0
	L2 = 0.5d0

	ri = 0.30d0
	rg = 0.10d0

	pi = dacos(-1.0d0)

    ! *** CONTORNO SUL ***

	! cilíndrico
    J = 1
    I = 1
    np = (j-1)*nx + i
    DX = L1 / NVK1
    X(np) = 0.0d0
    Y(np) = 0.0d0
    NK1 = NVK1
    DO I = 2, NK1
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
    end do
    NK2 = NK1 + 1
    i = nk2
    np = (j-1)*nx + i
    X(np) = L1
    Y(np) = 0.0d0
 
	! cossenóide
    DX = (L2 - L1) / NVK2
    NK3 = NK2 + 1
    NK4 = NK2 + NVK2 - 1
    DO I = NK3, NK4
       np = (j-1)*nx + i
       X(np) = X(np-1) + DX
       Y(np) = 0.0d0
	end do

    NK5 = NK4 + 1
    i = nk5
    np = (j-1)*nx + i
    X(np) = L2
    Y(np) = 0.0d0

    ! *** 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) = L2
       Y(np) = Y(np-nx) + DY
    end do
    j = ny-1
    np = (j-1)*nx + i
    X(np) = L2
    Y(np) = RI
 
    ! *** CONTORNO NORTE ***
    
	!cilíndrico
	J = ny-1
    NK1 = NVK1
    DO I = 2, NK1
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       Y(np) = RI
    end do
    NK2 = NK1 + 1
    i = nk2
    np = (j-1)*nx + i
    X(np) = L1
    Y(np) = RI

	! cossenóide
    NK3 = NK2 + 1
    NK4 = NK2 + NVK2 - 1
    DO I = NK3, NK4
       np  = (j-1)*nx + i
       nps = i
       X(np) = X(nps)
       DY = ( RI - RG ) * ( 1.0d0 + dcos(2.0d0*pi*(x(np)-L1) / (L2-L1)) ) / 2.0d0
	   Y(np) = RG + DY
    end do

	! índice i da garganta
	ig = (NK4+NK3)/2 

	rc = ( (L2-L1) ** 2 ) / ( 2.0d0 * (pi**2) * ( ri - rg ) )

  end subroutine cossenoide

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

  subroutine gera_malha

    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

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

  subroutine parametros_tubeira

    integer :: npg,  npin, npex, npw, npe
    real*8  :: rain, raex, rcg,  yl,  y2l, hw, he

    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

    yl   = ( y(npe) - y(npw) ) / ( x(npe) - x(npw) )

    hw   = x(npg) - x(npw)

    he   = x(npe) - x(npg)

    y2l  = 2 * (he*y(npw) + hw*y(npe) - (hw+he)*y(npg)) / ((hw+he)*hw*he)

    rcg = ( ( 1 + (yl**2) ) ** (1.5d0) ) / y2l

    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 da 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

    integer :: k ! auxiliar

    nome = 'Mach2D_6p0_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 pre_processamento

    call parametros_tubeira

    call escreve_dominio

    call grafico_dominio

  end subroutine pre_processamento

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

end module malha
