module coeficientes

! Objetivos: calcular métricas da transformação de coordenadas,
!            calcular os coeficientes e fontes dos sistemas de
!               equações dos volumes reais,
!            calcular velocidades contravariantes nas faces internas, 
!            calcular coeficientes do SIMPLEC,
!            calcular correções de p, ro, u, v, Uce e Vcn com pl.

! Versão atual: 15 Mar 09

! Relação das sub-rotinas em ordem de apresentação:
!    1) x_y_nodal
!    2) metricas
!    3) raios
!    4) massa_especifica_nos
!    5) massa_especifica_faces
!    6) coeficientes_velocidades
!    7) fonte_velocidade_u
!    8) fonte_velocidade_v
!    9) coeficientes_fonte_temperatura
!   10) coeficientes_fonte_pressao
!   11) velocidades_contravariantes_faces
!   12) coeficientes_SIMPLEC
!   13) pressao_e_massa_especifica_com_pl
!   14) u_v_nos_reais_com_pl
!   15) Uce_Vcn_faces_internas_com_pl

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

use malha

implicit none

contains

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

  subroutine x_y_nodal

    integer :: npw, npe, nps, npsw

    xp = 0.0d0
    yp = 0.0d0

    do i = 2, nx-1

       do j = 2, ny-1

          np   = (j-1)*nx + i

          npw  = np - 1
          npe  = np + 1
          nps  = np - nx
          npsw = nps - 1

          xp(np) = ( x(np) + x(npw) + x(npsw) + x(nps) ) / 4

          yp(np) = ( y(np) + y(npw) + y(npsw) + y(nps) ) / 4

       end do

    end do

  end subroutine x_y_nodal

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

  subroutine metricas

    integer :: npw, nps, npe, npn, npnw
    real*8  :: xkp, ykp, xep, yep

    xe  = 0.0d0
    ye  = 0.0d0
    xen = 0.0d0
    yen = 0.0d0
    xk  = 0.0d0
    yk  = 0.0d0
    xke = 0.0d0
    yke = 0.0d0
    Jp  = 0.0d0
    Je  = 0.0d0
    Jn  = 0.0d0
    Alfae = 0.0d0
    Gaman = 0.0d0
    Betae = 0.0d0
    Betan = 0.0d0

    ! métricas na face leste

    do i = 1, nx-1

       do j = 2, ny-1

          np   = (j-1)*nx + i

          nps  = np - nx

          xe(np) = x(np) - x(nps)

          ye(np) = y(np) - y(nps)

          Alfae(np) = xe(np)**2 + ye(np)**2

       end do

    end do

    ! métricas na face norte

    do i = 2, nx-1

       do j = 1, ny-1

          np   = (j-1)*nx + i

          npw  = np - 1

          xk(np) = x(np) - x(npw)

          yk(np) = y(np) - y(npw)

          Gaman(np) = xk(np)**2 + yk(np)**2

       end do

    end do

    ! Jacobiano no centro do volume de controle

    do i = 2, nx-1

       do j = 2, ny-1

          np  = (j-1)*nx + i

          npw = np - 1
          nps = np - nx

          xkp = ( xk(np) + xk(nps) ) / 2

          ykp = ( yk(np) + yk(nps) ) / 2

          xep = ( xe(np) + xe(npw) ) / 2

          yep = ( ye(np) + ye(npw) ) / 2

          Jp(np) = 1 / dabs ( xkp * yep - ykp * xep )

       end do

    end do

    ! Jacobiano e outras métricas no centro da face leste

    do i = 2, nx-2

       do j = 2, ny-1

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

          Je(np)  = ( Jp(np) + Jp(npe) ) / 2

          xke(np) = xp(npe) - xp(np)

          yke(np) = yp(npe) - yp(np)

       end do

    end do

    i = 1

    do j = 2, ny-1

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

       Je(np)  = Jp(npe)

       xke(np) = ( xk(npe) + xk(npe-nx) ) / 2

       yke(np) = ( yk(npe) + yk(npe-nx) ) / 2

    end do

    i = nx-1

    do j = 2, ny-1

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

       Je(np)  = Jp(np)

       xke(np) = ( xk(np) + xk(nps) ) / 2

       yke(np) = ( yk(np) + yk(nps) ) / 2

    end do

    ! Jacobiano e outras métricas no centro da face norte

    do i = 2, nx-1

       do j = 2, ny-2

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

          Jn(np) = ( Jp(np) + Jp(npn) ) / 2

          xen(np) = xp(npn) - xp(np)

          yen(np) = yp(npn) - yp(np)

       end do

    end do

    j = 1

    do i = 2, nx-1

       np   = (j-1)*nx + i
       npn  = np  + nx
       npnw = npn - 1

       Jn(np) = Jp(npn)

       xen(np) = ( xe(npn) + xe(npnw) ) / 2

       yen(np) = ( ye(npn) + ye(npnw) ) / 2

    end do

    j = ny-1

    do i = 2, nx-1

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

       Jn(np) = Jp(np)

       xen(np) = ( xe(np) + xe(npw) ) / 2

       yen(np) = ( ye(np) + ye(npw) ) / 2

    end do

    ! Beta do tensor métrico no centro da face leste

    do i = 2, nx-2

       do j = 2, ny-1

          np  = (j-1)*nx + i

          Betae(np) = xke(np)*xe(np) + yke(np)*ye(np)

       end do

    end do

    i = 1

    do j = 2, ny-1

       np  = (j-1)*nx + i

       Betae(np) = Betae(np+1)

    end do

    i = nx-1

    do j = 2, ny-1

       np  = (j-1)*nx + i

       Betae(np) = Betae(np-1)

    end do

    ! Beta do tensor métrico no centro da face norte

    do i = 2, nx-1

       do j = 2, ny-2

          np  = (j-1)*nx + i

          Betan(np) = xk(np)*xen(np) + yk(np)*yen(np)

       end do

    end do

    j = 1

    do i = 2, nx-1

       np  = (j-1)*nx + i

       Betan(np) = Betan(np+nx)

    end do

    j = ny-1

    do i = 2, nx-1

       np  = (j-1)*nx + i

       Betan(np) = Betan(np-nx)

    end do

  end subroutine metricas

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

  subroutine raios

    integer :: npw, npe, nps, npsw, npww, npee, npss, npnn, npn
    real*8  :: dr

    ! caso de geometria plana e inicialização

    raio = 1.0d0
	re   = 1.0d0
    rn   = 1.0d0
    rp   = 1.0d0

    ! caso de geometria axissimétrica

    if ( coord == 1 ) then

       ! canto nordeste

	   raio = y

       ! centro do volume de controle

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i

             npw  = np - 1
             npe  = np + 1
             nps  = np - nx
             npsw = nps - 1

             rp(np) = ( raio(np) + raio(npw) + raio(npsw) + raio(nps) ) / 4

          end do

       end do

       ! rp nos fictícios para os fontes Su5 e Su6

       ! fictícios do contorno esquerdo (entrada)
       i = 1
       do j = 2, ny-1
          np   = (j-1)*nx + i
          npe  = np + 1
          npee = np + 2
          dr = rp(npee) - rp(npe)
          rp(np) = rp(npe) - dr
       end do

       ! contorno direito (saída)
       i = nx
       do j = 2, ny-1
          np   = (j-1)*nx + i
          npw  = np - 1
          npww = np - 2
          dr = rp(npw) - rp(npww)
          rp(np) = rp(npw) + dr
       end do

       ! contorno inferior (linha de simetria)
       j = 1
       do i = 1, nx
          np   = (j-1)*nx + i
          npn  = np + nx
          npnn = npn + nx
          dr = rp(npnn) - rp(npn)
          rp(np) = rp(npn) - dr
       end do

       ! contorno superior (parede)
       j = ny
       do i = 1, nx
          np   = (j-1)*nx + i
          nps  = np - nx
          npss = nps - nx
          dr = rp(nps) - rp(npss)
          rp(np) = rp(nps) + dr
       end do

       ! centro da face leste

       do i = 1, nx-1

          do j = 2, ny-1

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

             re(np) = ( raio(np) + raio(nps) ) / 2

          end do

       end do

       ! centro da face norte

       do i = 2, nx-1

          do j = 1, ny-1

             np   = (j-1)*nx + i

             npw  = np - 1

             rn(np) = ( raio(np) + raio(npw) ) / 2

          end do

       end do

    end if

  end subroutine raios

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

  subroutine massa_especifica_nos

    ro = p / ( Rg * T )

  end subroutine massa_especifica_nos

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

  subroutine massa_especifica_faces

    integer :: npe,npn
    real*8  :: ae, an ! coeficientes alfa do UDS

    if ( beta == 1.0d0 ) then 

       do i = 1, nx-1
          do j = 1, ny-1
             np = (j-1)*nx + i
             npe  = np + 1
             npn  = np + nx
             roe(np) = ( ro(np) + ro(npe) ) / 2
             ron(np) = ( ro(np) + ro(npn) ) / 2
          end do
       end do

    else

       do i = 1, nx-1
          do j = 1, ny-1
             np = (j-1)*nx + i
             npe  = np + 1
             npn  = np + nx
             ae = dsign(0.5d0,Uce(np))
             roe(np) = (0.5d0+ae)*ro(np) + (0.5d0-ae)*ro(npe) &
                     + ae*(ro(npe)-ro(np))*beta
             an = dsign(0.5d0,Vcn(np))
             ron(np) = (0.5d0+an)*ro(np) + (0.5d0-an)*ro(npn) &
                     + an*(ro(npn)-ro(np))*beta
          end do
       end do

    end if

  end subroutine massa_especifica_faces

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

  subroutine coeficientes_velocidades

    integer :: npw, nps ! VC vizinhos a np
    real*8  :: fme, fmw, fmn, fms ! fluxos de massa
    real*8  :: ae,  aw,  an,  as  ! coeficientes alfa do UDS
    real*8  :: mpa ! massa antiga do VC np

    ! volumes reais

    if ( modvis == 0 ) then ! Euler

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i

             npw  = np - 1
             nps  = np - nx

             mpa = rp(np) * roa(np) / Jp(np)

             fmw = re(npw) * roe(npw) * Uce(npw)
             fme = re(np)  * roe(np)  * Uce(np)
             fms = rn(nps) * ron(nps) * Vcn(nps)
             fmn = rn(np)  * ron(np)  * Vcn(np)

             aw = dsign(0.5d0,Uce(npw))
             ae = dsign(0.5d0,Uce(np))
             as = dsign(0.5d0,Vcn(nps))
             an = dsign(0.5d0,Vcn(np))

             a(np,1) = 0
             a(np,3) = 0
             a(np,7) = 0
             a(np,9) = 0

             a(np,2) = - fms * ( 0.5d0 + as )
             a(np,4) = - fmw * ( 0.5d0 + aw )
             a(np,6) =   fme * ( 0.5d0 - ae )
             a(np,8) =   fmn * ( 0.5d0 - an )

          end do

       end do

    else ! Navier-Stokes

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i

             npw  = np - 1
             nps  = np - nx

             mpa = rp(np) * roa(np) / Jp(np)

             fmw = re(npw) * roe(npw) * Uce(npw)
             fme = re(np)  * roe(np)  * Uce(np)
             fms = rn(nps) * ron(nps) * Vcn(nps)
             fmn = rn(np)  * ron(np)  * Vcn(np)

             aw = dsign(0.5d0,Uce(npw))
             ae = dsign(0.5d0,Uce(np))
             as = dsign(0.5d0,Vcn(nps))
             an = dsign(0.5d0,Vcn(np))

             a(np,1) =   ( re(npw) * Je(npw) * Betae(npw) * vle(npw) &
                       +   rn(nps) * Jn(nps) * Betan(nps) * vln(nps) ) / 4

             a(np,3) = - ( re(np ) * Je(np ) * Betae(np ) * vle(np ) &
                       +   rn(nps) * Jn(nps) * Betan(nps) * vln(nps) ) / 4

             a(np,7) = - ( re(npw) * Je(npw) * Betae(npw) * vle(npw) &
                       +   rn(np ) * Jn(np ) * Betan(np ) * vln(np ) ) / 4

             a(np,9) =   ( re(np ) * Je(np ) * Betae(np ) * vle(np ) &
                       +   rn(np ) * Jn(np ) * Betan(np ) * vln(np ) ) / 4

             a(np,2) = - fms * ( 0.5d0 + as ) &
                       -   rn(nps) * Jn(nps) * Gaman(nps) * vln(nps) &
                       + ( re(npw) * Je(npw) * Betae(npw) * vle(npw) &
                       -   re(np ) * Je(np ) * Betae(np ) * vle(np ) ) / 4

             a(np,4) = - fmw * ( 0.5d0 + aw ) &
                       -   re(npw) * Je(npw) * Alfae(npw) * vle(npw) &
                       + ( rn(nps) * Jn(nps) * Betan(nps) * vln(nps) &
                       -   rn(np ) * Jn(np ) * Betan(np ) * vln(np ) ) / 4

             a(np,6) =   fme * ( 0.5d0 - ae ) &
                       -   re(np ) * Je(np ) * Alfae(np ) * vle(np ) &
                       + ( rn(np ) * Jn(np ) * Betan(np ) * vln(np ) &
                       -   rn(nps) * Jn(nps) * Betan(nps) * vln(nps) ) / 4

             a(np,8) =   fmn * ( 0.5d0 - an ) &
                       -   rn(np ) * Jn(np ) * Gaman(np ) * vln(np ) &
                       + ( re(np ) * Je(np ) * Betae(np ) * vle(np ) &
                       -   re(npw) * Je(npw) * Betae(npw) * vle(npw) ) / 4

          end do

       end do

    end if

  end subroutine coeficientes_velocidades

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

  subroutine fonte_velocidade_u

    integer :: npw, npe, nps, npn, npne, npse, npnw, npsw ! VC vizinhos a NP
    real*8  :: fme, fmw, fmn, fms ! fluxos de massa
    real*8  :: ae,  aw,  an,  as  ! coeficientes alfa do UDS
    real*8  :: mpa ! massa antiga do VC NP
    real*8  :: Pup ! termo fonte de pressão de u do VC NP
    real*8  :: S1, S2, S3, S4, S5, S6, S1p, S2p

    ! volumes reais, sem correção adiada, só invíscido

    do i = 2, nx-1

       do j = 2, ny-1

          np   = (j-1)*nx + i
          npw  = np - 1
          npe  = np + 1
          nps  = np - nx
          npn  = np + nx

          mpa = rp(np) * roa(np) / Jp(np)

          Pup = rp(np) * ( (p(np)+p(npn))*yk(np) - (p(np)+p(nps))*yk(nps) &
              + (p(np)+p(npw))*ye(npw) - (p(np)+p(npe))*ye(np) ) / 2

          b(np) = mpa * ua(np) / dt + Pup

          a(np,5) = mpa/dt - (a(np,2) + a(np,4) + a(np,6) + a(np,8))

       end do

    end do

    ! acréscimo devido à correção adiada

    if ( beta > 0.0d0 ) then

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i
             npw  = np - 1
             npe  = np + 1
             nps  = np - nx
             npn  = np + nx

             fmw = re(npw) * roe(npw) * Uce(npw)
             fme = re(np)  * roe(np)  * Uce(np)
             fms = rn(nps) * ron(nps) * Vcn(nps)
             fmn = rn(np)  * ron(np)  * Vcn(np)

             aw = dsign(0.5d0,Uce(npw))
             ae = dsign(0.5d0,Uce(np))
             as = dsign(0.5d0,Vcn(nps))
             an = dsign(0.5d0,Vcn(np))

             cup(np) = (fme*ae*(u(np)-u(npe)) + fmw*aw*(u(np)-u(npw)) &
                     +  fmn*an*(u(np)-u(npn)) + fms*as*(u(np)-u(nps)) ) * beta

             b(np) = b(np) + cup(np)

          end do

       end do

    end if

    ! acréscimo devido à viscosidade laminar

    if ( modvis == 1 ) then

       do j = 2, ny-1

          do i = 2, nx-1

             np   = (j-1)*nx + i
             npw  = np - 1
             npe  = np + 1
             nps  = np - nx
             npn  = np + nx
             npne = npn + 1
             npnw = npn - 1
             npsw = nps - 1
             npse = nps + 1

             S1 = (1.0d0/3) * (vle(np)*Je(np)) * ( (ye(np)**2)*u(npe)  &
                - (yke(np)*ye(np))*(u(npn)+u(npne)-u(nps)-u(npse))/4 ) &
                + (1.0d0/3) * (vle(npw)*Je(npw)) * ( (ye(npw)**2)*u(npw) &
                + (yke(npw)*ye(npw))*(u(npn)+u(npnw)-u(nps)-u(npsw))/4 )

             S1p = vle(np ) * Je(np ) * (ye(np )**2) &
                 + vle(npw) * Je(npw) * (ye(npw)**2)

             S2 = (1.0d0/3) * (vln(np)*Jn(np)) * ( (yk(np)**2)*u(npn)  &
                - (yk(np)*yen(np))*(u(npe)+u(npne)-u(npw)-u(npnw))/4 ) &
                + (1.0d0/3) * (vln(nps)*Jn(nps)) * ( (yk(nps)**2)*u(nps) &
                + (yk(nps)*yen(nps))*(u(npe)+u(npse)-u(npw)-u(npsw))/4 )

             S2p = vln(np ) * Jn(np ) * (yk(np )**2) &
                 + vln(nps) * Jn(nps) * (yk(nps)**2)

             S3 = (1/rp(np))*(re(np)*vle(np)*xe(np)*Je(np)) &
                * ( yke(np)*(v(npn)+v(npne)-v(nps)-v(npse))/4 - ye(np)*(v(npe)-v(np)) ) &
                - (1/rp(np))*(re(npw)*vle(npw)*xe(npw)*Je(npw)) &
                * ( yke(npw)*(v(npn)+v(npnw)-v(nps)-v(npsw))/4 - ye(npw)*(v(np)-v(npw)) )

             S4 = (1/rp(np))*(rn(np)*vln(np)*xk(np)*Jn(np)) &
                * ( yen(np)*(v(npe)+v(npne)-v(npw)-v(npnw))/4 - yk(np)*(v(npn)-v(np)) ) &
                - (1/rp(np))*(rn(nps)*vln(nps)*xk(nps)*Jn(nps)) &
                * ( yen(nps)*(v(npe)+v(npse)-v(npw)-v(npsw))/4 - yk(nps)*(v(np)-v(nps)) )

             S5 = - (2.0d0/3)*(1/re(np))*(vle(np)*ye(np)*Je(np)) &
                * ( xke(np)*(rp(npn)*v(npn)+rp(npne)*v(npne)-rp(nps)*v(nps)-rp(npse)*v(npse))/4 &
                    - xe(np)*(rp(npe)*v(npe)-rp(np)*v(np)) ) &
                  + (2.0d0/3)*(1/re(npw))*(vle(npw)*ye(npw)*Je(npw)) &
                * ( xke(npw)*(rp(npn)*v(npn)+rp(npnw)*v(npnw)-rp(nps)*v(nps)-rp(npsw)*v(npsw))/4 &
                    - xe(npw)*(rp(np)*v(np)-rp(npw)*v(npw)) )

             if ( j == 2 ) then

                S6 = - (2.0d0/3)*(1/rn(np))*(vln(np)*yk(np)*Jn(np)) &
                   * ( xen(np)*(rp(npe)*v(npe)+rp(npne)*v(npne)-rp(npw)*v(npw)-rp(npnw)*v(npnw))/4 &
                     - xk(np)*(rp(npn)*v(npn)-rp(np)*v(np)) )

             else

                S6 = - (2.0d0/3)*(1/rn(np))*(vln(np)*yk(np)*Jn(np)) &
                   * ( xen(np)*(rp(npe)*v(npe)+rp(npne)*v(npne)-rp(npw)*v(npw)-rp(npnw)*v(npnw))/4 &
                     - xk(np)*(rp(npn)*v(npn)-rp(np)*v(np)) ) &
                     + (2.0d0/3)*(1/rn(nps))*(vln(nps)*yk(nps)*Jn(nps)) &
                   * ( xen(nps)*(rp(npe)*v(npe)+rp(npse)*v(npse)-rp(npw)*v(npw)-rp(npsw)*v(npsw))/4 &
                     - xk(nps)*(rp(np)*v(np)-rp(nps)*v(nps)) )

             end if

             Sup(np) = rp(np) * (S1 + S2 + S3 + S4 + S5 + S6)

             b(np) = b(np) + Sup(np)

             a(np,5) = a(np,5) + rp(np)*(S1p+S2p)/3

          end do

       end do

    end if

    apu = a(:,5)

  end subroutine fonte_velocidade_u

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

  subroutine fonte_velocidade_v

    integer :: npw, npe, nps, npn, npne, npse, npnw, npsw ! VC vizinhos a NP
    real*8  :: fme, fmw, fmn, fms ! fluxos de massa
    real*8  :: ae,  aw,  an,  as  ! coeficientes alfa do UDS
    real*8  :: mpa ! massa antiga do VC NP
    real*8  :: Pvp ! termo fonte de pressão de v do VC NP
    real*8  :: S1, S2, S3, S4, S5, S6, S1p, S2p, S7p, S8

    ! volumes reais, sem correção adiada, só invíscido

    do i = 2, nx-1

       do j = 2, ny-1

          np   = (j-1)*nx + i
          npw  = np - 1
          npe  = np + 1
          nps  = np - nx
          npn  = np + nx

          mpa = rp(np) * roa(np) / Jp(np)

          Pvp = rp(np) * ( (p(np)+p(npe))*xe(np) - (p(np)+p(npw))*xe(npw) &
              + (p(np)+p(nps))*xk(nps) - (p(np)+p(npn))*xk(np) ) / 2

          b(np) = mpa * va(np) / dt + Pvp

          a(np,5) = mpa/dt - (a(np,2) + a(np,4) + a(np,6) + a(np,8))

       end do

    end do

    ! acréscimo devido à correção adiada

    if ( beta > 0.0d0 ) then

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i
             npw  = np - 1
             npe  = np + 1
             nps  = np - nx
             npn  = np + nx

             fmw = re(npw) * roe(npw) * Uce(npw)
             fme = re(np)  * roe(np)  * Uce(np)
             fms = rn(nps) * ron(nps) * Vcn(nps)
             fmn = rn(np)  * ron(np)  * Vcn(np)

             aw = dsign(0.5d0,Uce(npw))
             ae = dsign(0.5d0,Uce(np))
             as = dsign(0.5d0,Vcn(nps))
             an = dsign(0.5d0,Vcn(np))

             cvp(np) = (fme*ae*(v(np)-v(npe)) + fmw*aw*(v(np)-v(npw)) &
                     +  fmn*an*(v(np)-v(npn)) + fms*as*(v(np)-v(nps)) ) * beta

             b(np) = b(np) + cvp(np)

          end do

       end do

    end if

    ! acréscimo devido à viscosidade laminar

    if ( modvis == 1 ) then

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i
             npw  = np - 1
             npe  = np + 1
             nps  = np - nx
             npn  = np + nx
             npne = npn + 1
             npnw = npn - 1
             npsw = nps - 1
             npse = nps + 1

             S1 = (1/(3*rp(np))) * re(np ) * vle(np ) * Je(np ) * ( (xe(np )**2)*v(npe) &
                - (xke(np)*xe(np))*(v(npn)+v(npne)-v(nps)-v(npse))/4 ) &
                + (1/(3*rp(np))) * re(npw) * vle(npw) * Je(npw) * ( (xe(npw)**2)*v(npw) &
                + (xke(npw)*xe(npw))*(v(npn)+v(npnw)-v(nps)-v(npsw))/4 )

             S1p = re(np ) * vle(np ) * Je(np ) * (xe(np )**2) &
                 + re(npw) * vle(npw) * Je(npw) * (xe(npw)**2)

             S2 = (1/(3*rp(np))) * rn(np ) * vln(np ) * Jn(np ) * ( (xk(np )**2)*v(npn)  &
                - (xk(np)*xen(np))*(v(npe)+v(npne)-v(npw)-v(npnw))/4 ) &
                + (1/(3*rp(np))) * rn(nps) * vln(nps) * Jn(nps) * ( (xk(nps)**2)*v(nps)  &
                + (xk(nps)*xen(nps))*(v(npe)+v(npse)-v(npw)-v(npsw))/4 )

             S2p = rn(np ) * vln(np ) * Jn(np ) * (xk(np )**2) &
                 + rn(nps) * vln(nps) * Jn(nps) * (xk(nps)**2)

             S3 = vle(np)*ye(np)*Je(np) &
                * ( xke(np)*(u(npn)+u(npne)-u(nps)-u(npse))/4 - xe(np)*(u(npe)-u(np)) ) &
                - vle(npw)*ye(npw)*Je(npw) &
                * ( xke(npw)*(u(npn)+u(npnw)-u(nps)-u(npsw))/4 - xe(npw)*(u(np)-u(npw)) )

             S4 = vln(np)*yk(np)*Jn(np) &
                * ( xen(np)*(u(npe)+u(npne)-u(npw)-u(npnw))/4 - xk(np)*(u(npn)-u(np)) ) &
                - vln(nps)*yk(nps)*Jn(nps) &
                * ( xen(nps)*(u(npe)+u(npse)-u(npw)-u(npsw))/4 - xk(nps)*(u(np)-u(nps)) )

             S5 = - (2.0d0/3)*vle(np)*xe(np)*Je(np) &
                * ( yke(np)*(u(npn)+u(npne)-u(nps)-u(npse))/4 - ye(np)*(u(npe)-u(np)) ) &
                  + (2.0d0/3)*vle(npw)*xe(npw)*Je(npw) &
                * ( yke(npw)*(u(npn)+u(npnw)-u(nps)-u(npsw))/4 - ye(npw)*(u(np)-u(npw)) )

             S6 = - (2.0d0/3)*vln(np)*xk(np)*Jn(np) &
                * ( yen(np)*(u(npe)+u(npne)-u(npw)-u(npnw))/4 - yk(np)*(u(npn)-u(np)) ) &
                  + (2.0d0/3)*vln(nps)*xk(nps)*Jn(nps) &
                * ( yen(nps)*(u(npe)+u(npse)-u(npw)-u(npsw))/4 - yk(nps)*(u(np)-u(nps)) )

             S7p = 4*coord*vlp(np)/(rp(np)*Jp(np))

             S8 = 2*coord*v(np)/(3*rp(np)) * ( xe(np)*(vlp(np)+vlp(npe))/2 &
                - xe(npw)*(vlp(np)+vlp(npw))/2 ) - 2*coord*v(np)/(3*rp(np)) &
                * ( xk(np)*(vlp(np)+vlp(npn))/2 - xk(nps)*(vlp(np)+vlp(nps))/2 ) 

             Svp(np) = rp(np) * (S1 + S2 + S3 + S4 + S5 + S6 + S8)

             b(np) = b(np) + Svp(np)

             a(np,5) = a(np,5) + (S1p + S2p + S7p)/3

          end do

       end do

    end if

  end subroutine fonte_velocidade_v

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

  subroutine coeficientes_fonte_temperatura

    integer :: npw, npe, nps, npn, npne, npse, npnw, npsw ! VC vizinhos a NP
    integer :: ii ! contador para espécies químicas
    real*8  :: fme, fmw, fmn, fms ! fluxos de massa
    real*8  :: ae,  aw,  an,  as  ! coeficientes alfa do UDS
    real*8  :: mpa ! massa antiga do VC NP
    real*8  :: Pup ! termo fonte de pressão de u do VC NP
    real*8  :: Pvp ! termo fonte de pressão de v do VC NP
    real*8  :: STp ! termo fonte de pressão de T do VC NP
    real*8  :: SYip ! termo fonte de espécies químicas de T no VC NP
    real*8  :: cTp ! acréscimo devido à correção adiada de b de T
    real*8  :: S1, S2, S3, S4, S5, SYiaux, fat_eq

    ! COEFICIENTES dos volumes reais

    if ( modvis == 0 ) then ! Euler

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i

             npw  = np - 1
             nps  = np - nx

             mpa = rp(np) * roa(np) / Jp(np)

             fmw = re(npw) * roe(npw) * Uce(npw)
             fme = re(np)  * roe(np)  * Uce(np)
             fms = rn(nps) * ron(nps) * Vcn(nps)
             fmn = rn(np)  * ron(np)  * Vcn(np)

             aw = dsign(0.5d0,Uce(npw))
             ae = dsign(0.5d0,Uce(np))
             as = dsign(0.5d0,Vcn(nps))
             an = dsign(0.5d0,Vcn(np))

             a(np,1) = 0
             a(np,3) = 0
             a(np,7) = 0
             a(np,9) = 0

             a(np,2) = - fms * ( 0.5d0 + as ) * cp(np)
             a(np,4) = - fmw * ( 0.5d0 + aw ) * cp(np)
             a(np,6) =   fme * ( 0.5d0 - ae ) * cp(np)
             a(np,8) =   fmn * ( 0.5d0 - an ) * cp(np)

          end do

       end do

    else ! Navier-Stokes

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i

             npw  = np - 1
             nps  = np - nx

             mpa = rp(np) * roa(np) / Jp(np)

             fmw = re(npw) * roe(npw) * Uce(npw)
             fme = re(np)  * roe(np)  * Uce(np)
             fms = rn(nps) * ron(nps) * Vcn(nps)
             fmn = rn(np)  * ron(np)  * Vcn(np)

             aw = dsign(0.5d0,Uce(npw))
             ae = dsign(0.5d0,Uce(np))
             as = dsign(0.5d0,Vcn(nps))
             an = dsign(0.5d0,Vcn(np))

             a(np,1) =   ( re(npw) * Je(npw) * Betae(npw) * ke(npw) &
                       +   rn(nps) * Jn(nps) * Betan(nps) * kn(nps) ) / 4

             a(np,3) = - ( re(np ) * Je(np ) * Betae(np ) * ke(np ) &
                       +   rn(nps) * Jn(nps) * Betan(nps) * kn(nps) ) / 4

             a(np,7) = - ( re(npw) * Je(npw) * Betae(npw) * ke(npw) &
                       +   rn(np ) * Jn(np ) * Betan(np ) * kn(np ) ) / 4

             a(np,9) =   ( re(np ) * Je(np ) * Betae(np ) * ke(np ) &
                       +   rn(np ) * Jn(np ) * Betan(np ) * kn(np ) ) / 4

             a(np,2) = - fms * ( 0.5d0 + as ) * cp(np) &
                       -   rn(nps) * Jn(nps) * Gaman(nps) * kn(nps) &
                       + ( re(npw) * Je(npw) * Betae(npw) * ke(npw) &
                       -   re(np ) * Je(np ) * Betae(np ) * ke(np ) ) / 4

             a(np,4) = - fmw * ( 0.5d0 + aw ) * cp(np) &
                       -   re(npw) * Je(npw) * Alfae(npw) * ke(npw) &
                       + ( rn(nps) * Jn(nps) * Betan(nps) * kn(nps) &
                       -   rn(np ) * Jn(np ) * Betan(np ) * kn(np ) ) / 4

             a(np,6) =   fme * ( 0.5d0 - ae ) * cp(np) &
                       -   re(np ) * Je(np ) * Alfae(np ) * ke(np ) &
                       + ( rn(np ) * Jn(np ) * Betan(np ) * kn(np ) &
                       -   rn(nps) * Jn(nps) * Betan(nps) * kn(nps) ) / 4

             a(np,8) =   fmn * ( 0.5d0 - an ) * cp(np) &
                       -   rn(np ) * Jn(np ) * Gaman(np ) * kn(np ) &
                       + ( re(np ) * Je(np ) * Betae(np ) * ke(np ) &
                       -   re(npw) * Je(npw) * Betae(npw) * ke(npw) ) / 4

          end do

       end do

    end if

    ! FONTES dos volumes reais, sem correção adiada, só invíscido

    do i = 2, nx-1

       do j = 2, ny-1

          np   = (j-1)*nx + i
          npw  = np - 1
          npe  = np + 1
          nps  = np - nx
          npn  = np + nx

          mpa = rp(np) * roa(np) / Jp(np)

          Pup = rp(np) * ( (p(np)+p(npn))*yk(np) - (p(np)+p(nps))*yk(nps) &
              + (p(np)+p(npw))*ye(npw) - (p(np)+p(npe))*ye(np) ) / 2

          Pvp = rp(np) * ( (p(np)+p(npe))*xe(np) - (p(np)+p(npw))*xe(npw) &
              + (p(np)+p(nps))*xk(nps) - (p(np)+p(npn))*xk(np) ) / 2

          STp = rp(np) * ( p(np) - pa(np) ) / (Jp(np) * dt ) &
              - u(np) * Pup - v(np) * Pvp

          b(np) = mpa * cp(np) * Ta(np) / dt + STp

          a(np,5) = mpa*cp(np)/dt - (a(np,2) + a(np,4) + a(np,6) + a(np,8))

       end do

    end do
 
    ! acréscimo devido ao equilíbrio químico local
    
    if ( modf == 3 ) then
        do i = 2, nx-1

           do j = 2, ny-1

         	  SYip = 0.0d0
              
              np   = (j-1)*nx + i
              npw  = np - 1
              npe  = np + 1
              nps  = np - nx
              npn  = np + nx

              fmw = re(npw) * roe(npw) * Uce(npw)
              fme = re(np)  * roe(np)  * Uce(np)
              fms = rn(nps) * ron(nps) * Vcn(nps)
              fmn = rn(np)  * ron(np)  * Vcn(np)
  
              aw = dsign(beta3,Uce(npw))
              ae = dsign(beta3,Uce(np))
              as = dsign(beta3,Vcn(nps))
              an = dsign(beta3,Vcn(np))

              do ii = 1, Ne
		         SYiaux =   fme*((0.5d0 - ae)*Yi(npe,ii)+(0.5d0 + ae)*Yi(np,ii))  &
		                  - fmw*((0.5d0 - aw)*Yi(np,ii) +(0.5d0 + aw)*Yi(npw,ii)) &
			              + fmn*((0.5d0 - an)*Yi(npn,ii)+(0.5d0 + an)*Yi(np,ii))  &
			              - fms*((0.5d0 - as)*Yi(np,ii) +(0.5d0 + as)*Yi(nps,ii))
               
		         SYip = SYip - hi(np,ii)*SYiaux
	    	  
	    	  end do
	    	  
	          if ( it <= Ite1 ) fat_eq = 0.0d0 
              if ( it >  Ite1 .and. it < Ite2 ) &
                 fat_eq = 1.0d0*(it-Ite1)/(Ite2-Ite1)
              if ( it >= Ite2 ) fat_eq = 1.0d0 
              
	    	  b(np) = b(np) + fat_eq * SYip
	    	  
	       end do
	       
	    end do
	    
    end if

    ! acréscimo devido à correção adiada

    if ( beta > 0.0d0 ) then

       do i = 2, nx-1
 
          do j = 2, ny-1

             np   = (j-1)*nx + i
             npw  = np - 1
             npe  = np + 1
             nps  = np - nx
             npn  = np + nx

             fmw = re(npw) * roe(npw) * Uce(npw)
             fme = re(np)  * roe(np)  * Uce(np)
             fms = rn(nps) * ron(nps) * Vcn(nps)
             fmn = rn(np)  * ron(np)  * Vcn(np)

             aw = dsign(0.5d0,Uce(npw))
             ae = dsign(0.5d0,Uce(np))
             as = dsign(0.5d0,Vcn(nps))
             an = dsign(0.5d0,Vcn(np))

             cTp = (fme*ae*(T(np)-T(npe)) + fmw*aw*(T(np)-T(npw)) &
                 +  fmn*an*(T(np)-T(npn)) + fms*as*(T(np)-T(nps)) ) * beta * cp(np)

             b(np) = b(np) + cTp

          end do

       end do

    end if

    ! acréscimo devido à viscosidade laminar

    if ( modvis == 1 ) then

       do i = 2, nx-1

          do j = 2, ny-1

             np   = (j-1)*nx + i
             npw  = np - 1
             npe  = np + 1
             nps  = np - nx
             npn  = np + nx
             npne = npn + 1
             npnw = npn - 1
             npsw = nps - 1
             npse = nps + 1

             S1 = ( ye(np)*(u(np)+u(npe)) - ye(npw)*(u(np)+u(npw)) &
                -   yk(np)*(u(np)+u(npn)) + yk(nps)*(u(np)+u(nps)) ) / 2
             S1 = 2 * vlp(np) * Jp(np) * ( S1**2 )

             S2 = ( xk(np)*(v(np)+v(npn)) - xk(nps)*(v(np)+v(nps)) &
                -   xe(np)*(v(np)+v(npe)) + xe(npw)*(v(np)+v(npw)) ) / 2
             S2 = 2 * vlp(np) * Jp(np) * ( S2**2 )

             S3 = 2*coord*vlp(np)*((v(np)/rp(np))**2)/Jp(np)

             S4 = ( ye(np )*(v(np)+v(npe)) - xe(np )*(u(np)+u(npe)) &
                -   ye(npw)*(v(np)+v(npw)) + xe(npw)*(u(np)+u(npw)) &
                +   xk(np )*(u(np)+u(npn)) - yk(np )*(v(np)+v(npn)) &
                -   xk(nps)*(u(np)+u(nps)) + yk(nps)*(v(np)+v(nps)) ) / 2
             S4 = vlp(np) * Jp(np) * ( S4**2 )

             S5 = Uce(np) - Uce(npw) + Vcn(np) - Vcn(nps) &
                + coord*v(np)/(rp(np)*Jp(np))
             S5 = - 2 * vlp(np) * Jp(np) * ( S5**2 ) / 3

             STp = S1 + S2 + S3 + S4 + S5

             b(np) = b(np) + rp(np)*STp

          end do

       end do

    end if

  end subroutine coeficientes_fonte_temperatura

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

  subroutine coeficientes_fonte_pressao

    integer :: npw,npe,nps,npn ! VC vizinhos a NP
    real*8  :: ae,  aw,  an,  as  ! coeficientes alfa do UDS

    ! volumes reais

    do i = 2, nx-1

       do j = 2, ny-1

          np   = (j-1)*nx + i
          npw  = np - 1
          npe  = np + 1
          nps  = np - nx
          npn  = np + nx

          aw = dsign(0.5d0,Uce(npw))
          ae = dsign(0.5d0,Uce(np))
          as = dsign(0.5d0,Vcn(nps))
          an = dsign(0.5d0,Vcn(np))

          ap(np,1) = - rn(nps) * ( (0.5d0+as)*ro(nps) + (0.5d0-as)*ro(np)  ) * dn(nps) &
                     - rn(nps) * ( (0.5d0+as)*Vcn(nps)/(Rg*T(nps)) )
          ap(np,2) = - re(npw) * ( (0.5d0+aw)*ro(npw) + (0.5d0-aw)*ro(np)  ) * de(npw) &
                     - re(npw) * ( (0.5d0+aw)*Uce(npw)/(Rg*T(npw)) )
          ap(np,4) = - re(np)  * ( (0.5d0+ae)*ro(np)  + (0.5d0-ae)*ro(npe) ) * de(np)  &
                     + re(np)  * ( (0.5d0-ae)*Uce(np) /(Rg*T(npe)) )
          ap(np,5) = - rn(np)  * ( (0.5d0+an)*ro(np)  + (0.5d0-an)*ro(npn) ) * dn(np)  &
                     + rn(np)  * ( (0.5d0-an)*Vcn(np) /(Rg*T(npn)) )

          ap(np,3) = rp(np) / ( Jp(np) * Rg * T(np) * dt ) &
                   + rn(nps) * ( (0.5d0+as)*ro(nps) + (0.5d0-as)*ro(np)  ) * dn(nps) &
                   - rn(nps) * ( (0.5d0-as)*Vcn(nps)/(Rg*T(np)) ) &
                   + re(npw) * ( (0.5d0+aw)*ro(npw) + (0.5d0-aw)*ro(np)  ) * de(npw) &
                   - re(npw) * ( (0.5d0-aw)*Uce(npw)/(Rg*T(np)) ) &
                   + re(np)  * ( (0.5d0+ae)*ro(np)  + (0.5d0-ae)*ro(npe) ) * de(np)  &
                   + re(np)  * ( (0.5d0+ae)*Uce(np) /(Rg*T(np)) ) &
                   + rn(np)  * ( (0.5d0+an)*ro(np)  + (0.5d0-an)*ro(npn) ) * dn(np)  &
                   + rn(np)  * ( (0.5d0+an)*Vcn(np) /(Rg*T(np)) )

          b(np) = rp(np) * ( roa(np) - ro(np) ) / ( Jp(np) * dt ) &
                + rn(nps) * ( (0.5d0+as)*ro(nps) + (0.5d0-as)*ro(np)  ) * Vcn(nps) &
                + re(npw) * ( (0.5d0+aw)*ro(npw) + (0.5d0-aw)*ro(np)  ) * Uce(npw) &
                - re(np)  * ( (0.5d0+ae)*ro(np)  + (0.5d0-ae)*ro(npe) ) * Uce(np)  &
                - rn(np)  * ( (0.5d0+an)*ro(np)  + (0.5d0-an)*ro(npn) ) * Vcn(np)  &
                + beta    * (                                                     &
                  re(npw)*aw*(ro(np)-ro(npw))*(de(npw)*(pl(npw)-pl(np))+Uce(npw)) &
                - re(np )*ae*(ro(npe)-ro(np))*(de(np )*(pl(np)-pl(npe))+Uce(np))  &
                + rn(nps)*as*(ro(np)-ro(nps))*(dn(nps)*(pl(nps)-pl(np))+Vcn(nps)) &
                - rn(np )*an*(ro(npn)-ro(np))*(dn(np )*(pl(np)-pl(npn))+Vcn(np))  )

       end do

    end do

  end subroutine coeficientes_fonte_pressao

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

  subroutine velocidades_contravariantes_faces

    integer :: npw,npe,nps,npn,npne,npse,npee,npnw,npnn,npsw,npsee,npnee,npnnw,npnne
    real*8  :: somaup, somaue, sapu, sapv, aux1, aux2, somavp, somave, somaun, somavn
    real*8  :: mpa, mea, mna, Pue, Pve, Pun, Pvn

    ! Uce nas faces internas da malha

    do i = 2, nx-2

       do j = 2, ny-1

          np = (j-1)*nx + i

          npw = np - 1
          npe = np + 1
          nps = np - nx
          npn = np + nx
          npne = npn + 1
          npse = nps + 1
          npee = npe + 1
          npsw = nps - 1
          npnw = npn - 1
          npsee = npse + 1
          npnee = npne + 1

          mpa = rp(np)  * roa(np)  / Jp(np)
          mea = rp(npe) * roa(npe) / Jp(npe)

          somaup = a(np,1)*u(npsw) + a(np,2)*u(nps)  + a(np,3)*u(npse) + a(np,4)*u(npw) &
                 + a(np,6)*u(npe)  + a(np,7)*u(npnw) + a(np,8)*u(npn)  + a(np,9)*u(npne)

          somaue = a(npe,1)*u(nps)  + a(npe,2)*u(npse) + a(npe,3)*u(npsee) + a(npe,4)*u(np) &
                 + a(npe,6)*u(npee) + a(npe,7)*u(npn)  + a(npe,8)*u(npne)  + a(npe,9)*u(npnee)

          somavp = a(np,1)*v(npsw) + a(np,2)*v(nps)  + a(np,3)*v(npse) + a(np,4)*v(npw) &
                 + a(np,6)*v(npe)  + a(np,7)*v(npnw) + a(np,8)*v(npn)  + a(np,9)*v(npne)

          somave = a(npe,1)*v(nps)  + a(npe,2)*v(npse) + a(npe,3)*v(npsee) + a(npe,4)*v(np) &
                 + a(npe,6)*v(npee) + a(npe,7)*v(npn)  + a(npe,8)*v(npne)  + a(npe,9)*v(npnee)

          sapu = apu(np) + apu(npe)

          sapv = a(np,5) + a(npe,5)

          Pue = re(np) * ( (p(np)-p(npe))*ye(np) &
              + (p(npn)+p(npne)-p(nps)-p(npse))*yke(np)/4 )

          Pve = re(np) * ( (p(npe)-p(np))*xe(np) &
              + (p(nps)+p(npse)-p(npn)-p(npne))*xke(np)/4 )

          ue(np) = ( (mpa+mea)*uea(np)/dt + cup(np) + cup(npe) &
                 + Sup(np) + Sup(npe) - somaup - somaue + 2*Pue ) / sapu

          ve(np) = ( (mpa+mea)*vea(np)/dt + cvp(np) + cvp(npe) &
                 + Svp(np) + Svp(npe) - somavp - somave + 2*Pve ) / sapv

          Uce(np) = ue(np) * ye(np) - ve(np) * xe(np)

          aux1 = a(np,1) + a(np,2) + a(np,3) + a(np,4) &
               + a(np,6) + a(np,7) + a(np,8) + a(np,9)

          aux2 = a(npe,1) + a(npe,2) + a(npe,3) + a(npe,4) &
               + a(npe,6) + a(npe,7) + a(npe,8) + a(npe,9)

          due(np) = 2 * re(np) * ye(np) / ( apu(np) + aux1 + apu(npe) + aux2 )

          dve(np) = 2 * re(np) * xe(np) / ( a(np,5) + aux1 + a(npe,5) + aux2 )

       end do

    end do

    ! Vcn nas faces internas da malha

    do i = 2, nx-1

       do j = 2, ny-2

          np = (j-1)*nx + i

          npw = np - 1
          npe = np + 1
          nps = np - nx
          npn = np + nx
          npne = npn + 1
          npnw = npn - 1
          npnn = npn + nx
          npnnw = npnn - 1
          npnne = npnn + 1

          mpa = rp(np)  * roa(np)  / Jp(np)
          mna = rp(npn) * roa(npn) / Jp(npn)


          somaup = a(np,1)*u(npsw) + a(np,2)*u(nps)  + a(np,3)*u(npse) + a(np,4)*u(npw) &
                 + a(np,6)*u(npe)  + a(np,7)*u(npnw) + a(np,8)*u(npn)  + a(np,9)*u(npne)

          somaun = a(npn,1)*u(npw)  + a(npn,2)*u(np) + a(npn,3)*u(npe) + a(npn,4)*u(npnw) &
                 + a(npn,6)*u(npne) + a(npn,7)*u(npnnw) + a(npn,8)*u(npnn) + a(npn,9)*u(npnne)

          somavp = a(np,1)*v(npsw) + a(np,2)*v(nps)  + a(np,3)*v(npse) + a(np,4)*v(npw) &
                 + a(np,6)*v(npe)  + a(np,7)*v(npnw) + a(np,8)*v(npn)  + a(np,9)*v(npne)

          somavn = a(npn,1)*v(npw)  + a(npn,2)*v(np) + a(npn,3)*v(npe) + a(npn,4)*v(npnw) &
                 + a(npn,6)*v(npne) + a(npn,7)*v(npnnw) + a(npn,8)*v(npnn) + a(npn,9)*v(npnne)

          sapu = apu(np) + apu(npn)

          sapv = a(np,5) + a(npn,5)

          Pun = rn(np) * ( (p(npn)-p(np))*yk(np) &
              + (p(npw)+p(npnw)-p(npe)-p(npne))*yen(np)/4 )

          Pvn = rn(np) * ( (p(np)-p(npn))*xk(np) &
              + (p(npe)+p(npne)-p(npw)-p(npnw))*xen(np)/4 )

          un(np) = ( (mpa+mna)*una(np)/dt + cup(np) + cup(npn) &
                 + Sup(np) + Sup(npn) - somaup - somaun + 2*Pun ) / sapu

          vn(np) = ( (mpa+mna)*vna(np)/dt + cvp(np) + cvp(npn) &
                 + Svp(np) + Svp(npn) - somavp - somavn + 2*Pvn ) / sapv

          Vcn(np) = vn(np) * xk(np) - un(np) * yk(np)

          aux1 = a(np,1) + a(np,2) + a(np,3) + a(np,4) &
               + a(np,6) + a(np,7) + a(np,8) + a(np,9)

          aux2 = a(npn,1) + a(npn,2) + a(npn,3) + a(npn,4) &
               + a(npn,6) + a(npn,7) + a(npn,8) + a(npn,9)

          dun(np) = 2 * rn(np) * yk(np) / ( apu(np) + aux1 + apu(npn) + aux2 )

          dvn(np) = 2 * rn(np) * xk(np) / ( a(np,5) + aux1 + a(npn,5) + aux2 )

       end do

    end do

  end subroutine velocidades_contravariantes_faces

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

  subroutine coeficientes_SIMPLEC

    ! de nas faces internas da malha

    do i = 2, nx-2

       do j = 2, ny-1

          np  = (j-1)*nx + i

          de(np) = due(np) * ye(np) + dve(np) * xe(np)

       end do

    end do

    ! dn nas faces internas da malha

    do i = 2, nx-1

       do j = 2, ny-2

          np  = (j-1)*nx + i

          dn(np) = dun(np) * yk(np) + dvn(np) * xk(np)

       end do

    end do

  end subroutine coeficientes_SIMPLEC

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

  subroutine pressao_e_massa_especifica_com_pl

    ! correção da pressão com pl

    p = p + pl

    ! correção da massa especifica nodal com pl

    ro = p / ( Rg * T )
  
  end subroutine pressao_e_massa_especifica_com_pl

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

  subroutine u_v_nos_reais_com_pl

    integer :: npw, npe, nps, npn ! VC vizinhos a NP
    real*8  :: Supl ! termo fonte de pressão de u do VC NP
    real*8  :: Svpl ! termo fonte de pressão de v do VC NP
    real*8  :: aux

    ! volumes internos

    do i = 2, nx-1

       do j = 2, ny-1

          np = (j-1)*nx + i

          npw = np - 1
          npe = np + 1
          nps = np - nx
          npn = np + nx

          Supl = rp(np) * ( (pl(np)+pl(npn))*yk(np) - (pl(np)+pl(nps))*yk(nps) &
               + (pl(np)+pl(npw))*ye(npw) - (pl(np)+pl(npe))*ye(np) ) / 2

          Svpl = rp(np) * ( (pl(np)+pl(npe))*xe(np) - (pl(np)+pl(npw))*xe(npw) &
               + (pl(np)+pl(nps))*xk(nps) - (pl(np)+pl(npn))*xk(np) ) / 2

          aux = a(np,1) + a(np,2) + a(np,3) + a(np,4) &
              + a(np,6) + a(np,7) + a(np,8) + a(np,9)

          u(np) = u(np) + Supl / ( apu(np) + aux )

          v(np) = v(np) + Svpl / ( a(np,5) + aux )

       end do

    end do
  
  end subroutine u_v_nos_reais_com_pl

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

  subroutine Uce_Vcn_faces_internas_com_pl

    integer :: npe, npn ! VC vizinhos a NP

    ! velocidades Uce

    do i = 2, nx-2

       do j = 2, ny-1

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

          ue(np)  = ue(np)  + due(np) * (pl(np)  - pl(npe))
          ve(np)  = ve(np)  + dve(np) * (pl(npe) - pl(np))
          Uce(np) = ue(np) * ye(np) - ve(np) * xe(np)

       end do

    end do

    ! velocidades Vcn

    do i = 2, nx-1

       do j = 2, ny-2

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

          un(np)  = un(np)  + dun(np) * (pl(npn) - pl(np))
          vn(np)  = vn(np)  + dvn(np) * (pl(np)  - pl(npn))
          Vcn(np) = vn(np) * xk(np) - un(np) * yk(np)

       end do

    end do

  end subroutine Uce_Vcn_faces_internas_com_pl

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

end module coeficientes
