module coeficientes

! objetivo: calcular os coeficientes e termos fontes das equações discretizadas

use variaveis

implicit none

contains

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

  subroutine calcula_coeficientes_contornos

    ! volume fictício esquerdo
    aw(0) =  0
    ae(0) = -1
    ap(0) =  1
    bp(0) =  0

    ! volume fictício direito
    aw(N+1) = -1
    ae(N+1) =  0
    ap(N+1) =  1
    bp(N+1) =  2
  
  end subroutine calcula_coeficientes_contornos

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

  subroutine calcula_fonte

    select case ( fonte )

       case ( 0 )
          bp = 0

       case ( 1 )
          bp = Pe

       case ( 2 )
          bp = 2*Pe*X - 2

       case ( 3 )
          bp = 3*Pe*(X**2) - 6*X

       case ( 4 )
          bp = 4*Pe*(X**3) - 12*(X**2)

    end select
  
  end subroutine calcula_fonte

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

  subroutine calcula_coeficientes_uds

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    do P = 1, N

       aw(P) = 1 + Pe*h

       ae(P) = 1

       ap(P) = aw(P) + ae(P)

       bp(P) = bp(P)*(h**2)

    end do

  end subroutine calcula_coeficientes_uds

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

  subroutine calcula_coeficientes_cds

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    do P = 1, N

       aw(P) = 2 + Pe*h

       ae(P) = 2 - Pe*h

       ap(P) = aw(P) + ae(P)

       bp(P) = 2*bp(P)*(h**2)

    end do

  end subroutine calcula_coeficientes_cds

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

  subroutine calcula_coeficientes_wuds

    call calcula_fonte

    ! volume fictício esquerdo
    aw(0) =  0
    ae(0) =  alfa - 0.5q0
    ap(0) =  0.5q0 + alfa
    bp(0) =  0

    ! volumes reais
    do P = 1, N
       aw(P) = beta + (0.5q0+alfa)*Pe*h
       ae(P) = beta - (0.5q0-alfa)*Pe*h
       ap(P) = aw(P) + ae(P)
       bp(P) = bp(P)*(h**2)
    end do

    ! volume fictício direito
    aw(N+1) = -(0.5q0 + alfa)
    ae(N+1) =   0
    ap(N+1) =   0.5q0 - alfa
    bp(N+1) =   1

  end subroutine calcula_coeficientes_wuds

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

  subroutine calcula_coeficientes_plds

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    do P = 1, N

       aw(P) = pow + Pe*h

       ae(P) = pow

       ap(P) = aw(P) + ae(P)

       bp(P) = bp(P)*(h**2)

    end do

  end subroutine calcula_coeficientes_plds

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

  subroutine calcula_coeficientes_quick

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    P = 1

       aw(P) = 0

       ae(P) = 4

       ap(P) = 12 + 3*Pe*h

    do P = 2, N-1

       aw(P) = 1 + Pe*h

       ae(P) = 1

       ap(P) = aw(P) + ae(P)

    end do

    P = N

       aw(P) =  4 + 3*Pe*h

       ae(P) =  0

       ap(P) = 12

  end subroutine calcula_coeficientes_quick

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

  subroutine calcula_fontes_quick

    call calcula_fonte

    ! volumes de controle reais

    P = 1

       bp(P) = 3*bp(P)*(h**2) + 3*Pe*h*( T(P) - 3*T(P+1) ) / 8.0q0

    P = 2

       bp(P) = bp(P)*(h**2) + Pe*h*( 5*T(P) - 3*T(P+1) ) / 8.0q0

    do P = 3, N-1

       bp(P) = bp(P)*(h**2) + Pe*h*( 5*T(P) - T(P-1) - T(P-2) - 3*T(P+1) ) / 8.0q0

    end do

    P = N

       bp(P) = 3*bp(P)*(h**2) + 8 + 3*Pe*h*( ( (3*T(P) - 2*T(P-1) &
             - T(P-2) ) / 8.0q0 ) - 1.0q0 )

  end subroutine calcula_fontes_quick

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

  subroutine calcula_coeficientes_uds2

    ! volumes de controle reais

    P = 1

       aw(P) = 1 + Pe*h/2

       ae(P) = 1

       ap(P) = 2 + 3*Pe*h/2.0q0

    do P = 2, N-1

       aw(P) = 1 + Pe*h

       ae(P) = 1

       ap(P) = aw(P) + ae(P)

    end do

    P = N

       aw(P) = 1 + Pe*h

       ae(P) = 1

       ap(P) = 2

  end subroutine calcula_coeficientes_uds2

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

  subroutine calcula_fontes_uds2

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    P = 1

       bp(P) = bp(P)*(h**2)

    do P = 2, N-1

       bp(P) = bp(P)*(h**2) + Pe*h*( 2*T(P-1) - T(P) - T(P-2) ) / 2.0q0

    end do

    P = N

       bp(P) = bp(P)*(h**2) + Pe*h*( ((T(P-1)-T(P-2))/2) - 1 )

  end subroutine calcula_fontes_uds2

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

  subroutine calcula_coeficientes_ads

    real*16 :: Tmp, Tmn, Tmc, Tme

    ! cálculo dos coeficientes alfm

    alfm(0) = 0

    do P = 1, N-1
       if ( T(P+1) == T(P) ) then
          alfm(P) = 0.5q0
       else
          Tmn = T(P)   + 0.5q0*(T(P)   - T(P-1))
          Tmc =          0.5q0*(T(P)   + T(P+1))
          Tmp = T(P+1) - 0.5q0*(T(P+2) - T(P+1))
          Tme = qmax1( qmin1(Tmn,Tmc), qmin1(Tmn,Tmp), qmin1(Tmc,Tmp) )
          alfm(P) = 0.5q0*(2*Tme-T(P)-T(P+1))/(T(P)-T(P+1))
          alfm(P) = qmin1 ( 0.5q0, qmax1(-0.5q0,alfm(P)) )
       end if
    end do

    alfm(N) = 0

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    do P = 1, N

       aw(P) = 1 + (0.5q0 + alfm(P-1))*Pe*h

       ae(P) = 1 - (0.5q0 - alfm(P)  )*Pe*h

       ap(P) = aw(P) + ae(P)

       bp(P) = bp(P)*(h**2)

    end do

  end subroutine calcula_coeficientes_ads

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

  subroutine calcula_coeficientes_tvd

    real*16 :: psi, r

    ! cálculo dos coeficientes alfm

    alfm(0) = 0

    do P = 1, N-1
       if ( T(P+1) == T(P) ) then
          alfm(P) = 0.5q0
       else
          r   = ( T(P) - T(P-1) ) / ( T(P+1) - T(P) )
          psi = qmax1(0.0q0,qmin1(2*r,1.0q0),qmin1(r,2.0q0))
          alfm(P) = 0.5q0*(1 - psi)
       end if
    end do

    alfm(N) = 0

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    do P = 1, N

       aw(P) = 1 + (0.5q0 + alfm(P-1))*Pe*h

       ae(P) = 1 - (0.5q0 - alfm(P)  )*Pe*h

       ap(P) = aw(P) + ae(P)

       bp(P) = bp(P)*(h**2)

    end do

  end subroutine calcula_coeficientes_tvd

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

  subroutine calcula_coeficientes_alfa

    call calcula_fonte

    ! volume fictício esquerdo
    aw(0) =  0
    ae(0) =  alf - 0.5q0
    ap(0) =  0.5q0 + alf
    bp(0) =  0

    ! volumes reais
    do P = 1, N
       aw(P) = 1 + (0.5q0+alf)*Pe*h
       ae(P) = 1 - (0.5q0-alf)*Pe*h
       ap(P) = aw(P) + ae(P)
       bp(P) = bp(P)*(h**2)
    end do

    ! volume fictício direito
    aw(N+1) = -(0.5q0 + alf)
    ae(N+1) =   0
    ap(N+1) =   0.5q0 - alf
    bp(N+1) =   1

  end subroutine calcula_coeficientes_alfa

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

  subroutine calcula_coeficientes_cds4

    ! volumes de controle reais

    P = 1

       aw(P) = 0

       ae(P) = 1.0q0 + 81/24.0q0

       ap(P) = 9.0q0 + 81/24.0q0 + 3*Pe*h

    do P = 2, N-1

       aw(P) = 27/24.0q0 + Pe*h

       ae(P) = 27/24.0q0

       ap(P) = aw(P) + ae(P)

    end do

    P = N

       aw(P) = 1.0q0 + 81/24.0q0 + 3*Pe*h

       ae(P) = 0

       ap(P) = 9.0q0 + 81/24.0q0

  end subroutine calcula_coeficientes_cds4

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

  subroutine calcula_fontes_cds4

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    P = 1

       bp(P) = 3*bp(P)*(h**2) + 3*Pe*h*( T(P) - 3*T(P+1) ) / 8.0q0 &
             + ( T(P-1) - T(P+2) ) / 8.0q0

    P = 2

       bp(P) = bp(P)*(h**2) + Pe*h*( 5*T(P) - 3*T(P+1) ) / 8.0q0 &
             + ( T(P-1) + T(P+1) - T(P-2) - T(P+2) ) / 24.0q0

    do P = 3, N-1

       bp(P) = bp(P)*(h**2) + Pe*h*( 5*T(P) - T(P-1) - T(P-2) - 3*T(P+1) ) / 8.0q0 &
             + ( T(P-1) + T(P+1) - T(P-2) - T(P+2) ) / 24.0q0

    end do

    P = N

       bp(P) = 3*bp(P)*(h**2) - 3*Pe*h + 8 + ( T(P+1) - T(P-2) ) / 8.0q0 &
             + 3*Pe*h*( 3*T(P) - 2*T(P-1) - T(P-2) ) / 8.0q0  

  end subroutine calcula_fontes_cds4

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

  subroutine calcula_coeficientes_quick3

    ! volumes de controle reais

    P = 1

       aw(P) = 0

       ae(P) = 111/60.0q0

       ap(P) = 295/60.0q0 + Pe*h

    P = 2

       aw(P) = 70/60.0q0 + Pe*h

       ae(P) = 27/24.0q0

       ap(P) = 275/120.0q0 + Pe*h

    do P = 3, N-2

       aw(P) = 27/24.0q0 + Pe*h

       ae(P) = 27/24.0q0

       ap(P) = aw(P) + ae(P)

    end do

    P = N-1

       aw(P) = 27/24.0q0 + Pe*h

       ae(P) = 70/60.0q0

       ap(P) = 275/120.0q0 + Pe*h

    P = N

       aw(P) = 111/60.0q0 + Pe*h

       ae(P) = 0

       ap(P) = 295/60.0q0

  end subroutine calcula_coeficientes_quick3

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

  subroutine calcula_fontes_quick3

    call calcula_fonte

    call calcula_coeficientes_contornos

    ! volumes de controle reais

    P = 1

       bp(P) = bp(P)*(h**2) + Pe*h*( T(P) - 3*T(P+1) ) / 8.0q0 &
             + ( 9*T(P+1) - 5*T(P) - 12*T(P+2) ) / 60.0q0

    P = 2

       bp(P) = bp(P)*(h**2) + Pe*h*( 5*T(P) - 3*T(P+1) ) / 8.0q0 &
             + ( 15*T(P-1) + 6*T(P+1) - 5*T(P+2) ) / 120.0q0

    do P = 3, N-2

       bp(P) = bp(P)*(h**2) + Pe*h*( 5*T(P) - T(P-1) - T(P-2) - 3*T(P+1) ) / 8.0q0 &
             + ( T(P-1) + T(P+1) - T(P-2) - T(P+2) ) / 24.0q0

    end do

    P = N-1

       bp(P) = bp(P)*(h**2) + Pe*h*( 5*T(P) - T(P-1) - T(P-2) - 3*T(P+1) ) / 8.0q0 &
             + ( 15*T(P+1) + 6*T(P-1) - 5*T(P-2) - 16.0q0 ) / 120.0q0

    P = N

       bp(P) = bp(P)*(h**2) - Pe*h + Pe*h*( 3*T(P) - 2*T(P-1) - T(P-2) ) / 8.0q0 &
             + 184/60.0q0 + ( 9*T(P-1) + 8.0q0 - 12*T(P-2) - 5*T(P) ) / 60.0q0

  end subroutine calcula_fontes_quick3

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

end module coeficientes
