module calculos

! objetivo: calcular as variáveis.

use dados

implicit none

contains

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

  subroutine calcula

    call extrapolacoes_e_ordens_aparente

    call estimador_Richardson

    call multiextrapolacao
    
    call estimadores  !aplica cálculo de estimadores sobre solução com MER
    
    call estimadores1 !aplica cálculo de estimadores sobre solução sem MER

  end subroutine calcula

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

  subroutine extrapolacoes_e_ordens_aparente
  
    real*16 :: psi, psi1       ! taxa de convergência e equivalente
    real*16 :: del, del1       ! auxiliares

    ! pU_h = ordem aparente com base em Th
    pU_h = great
 	it_pU_h = 0
    if ( nm >= 3 ) then
       do i = 3, nm
          del = Th(i) - Th(i-1)
          if ( del /= 0.0q0 ) then
             psi = (Th(i-1)-Th(i-2))/(Th(i)-Th(i-1))
          else
             psi = 0.0q0
          end if
          if ( tipo_refino == 1 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             pU_h(i) = qlog10(psi)/qlog10(r21)
          end if
          if ( tipo_refino == 2 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             r32 = h(i-2) / h(i-1)
             call bissecao ( pU_min, pU_max, r21, r32, psi, it_pU_h(i), pU_h(i) )
          end if
       end do
    end if
    
    !Cálculo da ordem aparente equivalente: pUh1 
     if ( nm >= 3 ) then
         do i = 3, nm
            del = Th(i) - Th(i-1)
            if ( del /= 0.0q0 ) then
                psi = (Th(i-1)-Th(i-2))/del
                psi_v(i)=psi
            end if
            if ( tipo_refino == 1 ) then
                if ( (psi_v(i)/=0.0q0) .and. (psi_v(i)/=111) ) then
                    r21 = h(i-1) / h(i)
                    pU_h1(i) = qlog10(qabs(psi))/qlog10(r21)
                end if
            end if
         end do
      end if
    

    ! Ti_12 = Th extrapolado com as malhas 1 e 2 e pL
    do i = 2, nm
       r21 = h(i-1) / h(i)
       Ti_12(i) = Th(i) + (Th(i) - Th(i-1)) / ((r21**pL) - 1.0q0)
    end do

    ! Ti_13 = Th extrapolado com as malhas 1 e 3 e pL
    ! Ti_23 = Th extrapolado com as malhas 2 e 3 e pL
    ! Ti_pU = Th extrapolado com as malhas 1 a 3 e pU_h
    ! Ti_pU1 = Th extrapolado com as malhas 1 a 3 e pU_h equivalente
    ! Tc_12 = Média de Ti_12 com Ti_pU
    ! Tc_13 = Média de Ti_13 com Ti_pU
    if ( nm >= 3 ) then
       do i = 3, nm
          r21 = h(i-1) / h(i)
          r32 = h(i-2) / h(i-1)
          Ti_23(i) = Th(i-1) + (Th(i-1) - Th(i-2)) / ((r32**pL) - 1.0q0)
          Ti_13(i) = Th(i)   + (Th(i)   - Th(i-2)) / (((r21*r32)**pL) - 1.0q0)
          if ( pU_h(i) <= 0.0q0 .or. pU_h(i) > limite ) then
             Ti_pU(i) = great
             Ti_pU1(i) = great
             Tc_12(i) = great
             Tc_13(i) = great
          else
             Ti_pU(i) = Th(i) + (Th(i) - Th(i-1)) / ((r21**pU_h(i)) - 1.0q0)
             Ti_pU1(i)= Th(i) + (Th(i) - Th(i-1)) / ((r21**pU_h1(i)) - 1.0q0) 
             Tc_12(i) = ( Ti_12(i) + Ti_pU(i) ) / 2
             Tc_13(i) = ( Ti_13(i) + Ti_pU(i) ) / 2
          end if
          if ( pU_h(i) <= 0.0q0 ) then
             Ti_pU1(i) = Th(i) + (Th(i) - Th(i-1)) / ((r21**pU_h1(i)) - 1.0q0)
          end if
       end do
    end if

    ! pU_i = ordem aparente com base em Ti_pU
    ! pU_i1 = ordem aparente equivalente com base em Ti_pU*
    pU_i = great
    pU_i1 = great !Equivalente
    it_pU_i = 0
    if ( nm >= 5 ) then
       do i = 5, nm
          del  = Ti_pU(i) - Ti_pU(i-1)
          del1 = Ti_pU1(i) - Ti_pU1(i-1)
          if ( del /= 0.0q0 ) then
             psi = (Ti_pU(i-1)-Ti_pU(i-2))/(Ti_pU(i)-Ti_pU(i-1))
          else
             psi = great
          end if
          if ( del1 /= 0.0q0 ) then
             psi1 = (Ti_pU1(i-1)-Ti_pU1(i-2))/(Ti_pU1(i)-Ti_pU1(i-1))
          else
             psi1 = great
          end if
          if ( tipo_refino == 1 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             pU_i(i) = qlog10(psi)/qlog10(r21)
          end if
          if ( psi1 /= 0 .and. psi1 /= great) then
             r21 = h(i-1) / h(i)
             psi_i1(i) = psi1
             pU_i1(i) = qlog10(qabs(psi1))/qlog10(r21) !(*) Equivalente
          end if
          
          
          if ( tipo_refino == 2 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             r32 = h(i-2) / h(i-1)
             call bissecao ( pU_min, pU_max, r21, r32, psi, it_pU_i(i), pU_i(i) )
          end if
       end do
    end if

    ! pU_c12 = ordem aparente com base em Tc_12
    pU_c12 = great
    it_pU_c12 = 0
    if ( nm >= 5 ) then
       do i = 5, nm
          del = Tc_12(i) - Tc_12(i-1)
          if ( del /= 0.0q0 ) then
             psi = (Tc_12(i-1)-Tc_12(i-2))/(Tc_12(i)-Tc_12(i-1))
          else
             psi = great
          end if
          if ( tipo_refino == 1 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             pU_c12(i) = qlog10(psi)/qlog10(r21)
          end if
          if ( tipo_refino == 2 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             r32 = h(i-2) / h(i-1)
             call bissecao ( pU_min, pU_max, r21, r32, psi, it_pU_c12(i), pU_c12(i) )
          end if
       end do
    end if

    ! pU_c13 = ordem aparente com base em Tc_13
    pU_c13 = great
    it_pU_c13 = 0
    if ( nm >= 5 ) then
       do i = 5, nm
          del = Tc_13(i) - Tc_13(i-1)
          if ( del /= 0.0q0 ) then
             psi = (Tc_13(i-1)-Tc_13(i-2))/(Tc_13(i)-Tc_13(i-1))
          else
             psi = 0.0q0
          end if
          if ( tipo_refino == 1 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             pU_c13(i) = qlog10(psi)/qlog10(r21)
          end if
          if ( tipo_refino == 2 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             r32 = h(i-2) / h(i-1)
             call bissecao ( pU_min, pU_max, r21, r32, psi, it_pU_c13(i), pU_c13(i) )
          end if
       end do
    end if

    ! Tbi_pU  = Th biextrapolado com malhas 1 a 5 e pU_i
    ! Tbi_pU1 = Th biextrapolado com malhas 1 a 5 e pU_i equivalente
    if ( nm >= 5 ) then
       do i = 5, nm
          r21 = h(i-1) / h(i)
          if ( pU_i(i) <= 0.0q0 .or. pU_i(i) > limite ) then
             Tbi_pU(i) = great
             Tbi_pU1(i) = Tbi_pU(i)
          else
             Tbi_pU(i) = Ti_pU(i) + (Ti_pU(i) - Ti_pU(i-1)) / ((r21**pU_i(i)) - 1.0q0)
             Tbi_pU1(i) = Ti_pU1(i) + (Ti_pU1(i) - Ti_pU1(i-1)) / ((r21**pU_i1(i)) - 1.0q0)
          end if
          if(pU_i(i) < 0.0q0) then
            Tbi_pU1(i) = Ti_pU1(i) + (Ti_pU1(i) - Ti_pU1(i-1)) / ((r21**pU_i1(i)) - 1.0q0)
          end if
       end do
    end if

    ! pU_bi = ordem aparente com base em Tbi_pU
    pU_bi    = great
    pU_bi1   = great !Equivalente
    it_pU_bi = 0
    if ( nm >= 7 ) then
       do i = 7, nm
          del  = Tbi_pU(i) - Tbi_pU(i-1)
          del1 = Tbi_pU1(i) - Tbi_pU1(i-1)
          if ( del /= 0.0q0 ) then
             psi = (Tbi_pU(i-1)-Tbi_pU(i-2))/(Tbi_pU(i)-Tbi_pU(i-1))
          else
             psi = great
          end if
          if ( del1 /= 0.0q0 ) then
             psi1 = (Tbi_pU1(i-1)-Tbi_pU1(i-2))/(Tbi_pU1(i)-Tbi_pU1(i-1))
             psi_bi1(i)=psi1
          else
             psi1 = great
             psi_bi1(i)=psi1
          end if
          if ( tipo_refino == 1 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             pU_bi(i) = qlog10(psi)/qlog10(r21)
          end if
          if ( psi1 /= 0  .and. psi1 /= great ) then
             r21 = h(i-1) / h(i)
             psi_bi1(i)=psi1
             pU_bi1(i) = qlog10(qabs(psi1))/qlog10(r21)
          end if
          
          if ( tipo_refino == 2 .and. psi > 0.0q0 ) then
             r21 = h(i-1) / h(i)
             r32 = h(i-2) / h(i-1)
             call bissecao ( pU_min, pU_max, r21, r32, psi, it_pU_bi(i), pU_bi(i) )
          end if
       end do
    end if

  end subroutine extrapolacoes_e_ordens_aparente

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

  subroutine bissecao ( pUi, pUs, r21, r32, psi, i_pU, pUm )

    integer :: it   ! número de iterações corrente na obtenção de pU

    real*16 :: f_pUi ! resto da função de pU para pUi
    real*16 :: f_pUm ! resto da função de pU para pUm
    real*16 :: dp    ! intervalo de teste de pU
    real*16 :: pi    ! limite inferior corrente de pU
    real*16 :: ps    ! limite superior corrente de pU

    real*16,intent(in) :: pUi ! limite inferior de pU (dado inicial)
    real*16,intent(in) :: pUs ! limite superior de pU (dado inicial)
    real*16,intent(in) :: r21 ! razão de h entre as malhas 2 e 1
    real*16,intent(in) :: r32 ! razão de h entre as malhas 3 e 2
    real*16,intent(in) :: psi ! razão de convergência

    integer,intent(out) :: i_pU ! número de iterações usado na obtenção de pU
    real*16, intent(out) :: pUm  ! solução de pU

    pi = pUi
    ps = pUs

    do it = 1, numero_iteracoes
      pUm = 0.5q0 * (pi + ps)
      f_pUi = pi  * qlog10(r21) - qlog10(psi*((r21**pi )-1.0q0)/((r32**pi )-1.0q0))
      f_pUm = pUm * qlog10(r21) - qlog10(psi*((r21**pUm)-1.0q0)/((r32**pUm)-1.0q0))
      if ( f_pUi * f_pUm > 0.0q0 ) then
         pi = pUm
      else
         dp = ps - pi
         i_pU = it
         if ( dp <= tol ) exit
         ps = pUm
      end if
    end do

    if ( it == numero_iteracoes .or. it == (numero_iteracoes+1) ) pUm = great

  end subroutine bissecao


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

  subroutine estimador_Richardson 

    ! Uri_12 = erro estimado de Th com as malhas 1 e 2 e pL
    do i = 2, nm
       Uri_12(i) = Ti_12(i) - Th(i)
    end do

    ! Uri_13 = erro estimado de Th com as malhas 1 e 3 e pL
    ! Uri_23 = erro estimado de Th com as malhas 2 e 3 e pL
    ! Uri_pU = erro estimado de Th com as malhas 1 a 3 e pU_h
    ! Uri_p  = erro estimado de Th com as malhas 1 a 3, pU_h e pL
    ! Uc_12  = erro estimado de Tc_12
    ! Uc_13  = erro estimado de Tc_13
    if ( nm >= 3 ) then
       do i = 3, nm
          Uri_13(i) = Ti_13(i) - Th(i)
          Uri_23(i) = Ti_23(i) - Th(i)
          if ( Ti_pU(i) > limite ) then
             Uri_pU(i) = great
             Uri_p(i)  = great
             Uc_12(i)  = great
             Uc_13(i)  = great
          else
             Uri_pU(i) = Ti_pU(i) - Th(i)
             Uri_p(i)  = qsign ( qmax1 ( qabs(Uri_12(i)), qabs(Uri_pU(i)) ), &
                                       ( Th(i) - Th(i-1) ) )
             Uc_12(i)  = qabs ( Ti_12(i) - Ti_pU(i) ) / 2
             Uc_13(i)  = qabs ( Ti_13(i) - Ti_pU(i) ) / 2
          end if
       end do
    end if

  end subroutine estimador_Richardson

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

  subroutine multiextrapolacao

    ! Tm_12 = Th multiextrapolado com as malhas 1 e 2 e pV

    m = 1
    do i = 2, nm
       Tm_12(i,m) = Ti_12(i)
    end do

    do m = 2, nm-1

       do i = m+1, nm

          r21 = h(i-1) / h(i)

          Tm_12(i,m) = Tm_12(i,m-1) + ( Tm_12(i,m-1) - Tm_12(i-1,m-1) ) &
                     / (  r21 ** pV(m) - 1.0q0 )

       end do

    end do

    ! Um_12 = erro estimado de Tm_12 com as malhas 1 e 2 e pV

    do m = 1, nm-2

       do i = m+2, nm

          Um_12(i,m) = Tm_12(i,m+1) - Tm_12(i,m) 

       end do

    end do

  end subroutine multiextrapolacao

!--------------------------------------------------------------------------------------------
!aplica cálculo de estimador sobre solução com MER
!--------------------------------------------------------------------------------------------
  subroutine estimadores            

    real*16 :: psi1       ! taxa de convergência internível 
    real*16 :: del1,del2  ! variáveis auxiliares
           
    do i = 2, nm  
       r21 = h(i-1) / h(i)
       !Atribui valores ao vetor T_MER
       T_MER(i) = Tm_12(i,i-1) 
        if (( nm >= 3 ).and.(i>=3)) then
          !variáveis auxiliares
          del1 = T_MER(i-1) - T_MER(i-2) 
          del2 = T_MER(i) - T_MER(i-1)       
          !Cálculo da taxa de convergência interníveis (máximo de extrapolações em cada malha, U(Tm1))
          if ( del2 /= 0.0q0 ) then
              psi1 = del1/del2       
              psi_inter(i)=psi1 
              !Cálculo de PU MER equivalente interníveis
              if (psi1 /= 0.0q0) then           
                 pU_MER(i) = qlog10(qabs(psi1))/qlog10(r21)
              end if 
          else
              psi1 = great
              psi_inter(i)=psi1 
          end if
        end if                  
                                 
          if ( (qabs(psi1) > 1) .and. (psi1 /= great)  ) then
            ! Calcula Estimador Richardson com a razão de convergência interníveis média, U_Ri_mpsi_MER
            ! Utiliza m(PU MER), ordem aparente interníveis média através do cálculo de psi1 atribuido à malha intermediária 
              U_Ri_mpsi_MER(i-1) = ((T_MER(i-1)-T_MER(i-2)) / (psi1-1.0q0 )) !estimativa na malha 2 e PU_MER atribuida à malha 2       
            !Calcula U_Ri_mpsi_MER na malha mais fina
            if(i==nm) then
                psi1=(psi_inter(i)**2.0q0)/(psi_inter(i-1))
                U_Ri_mpsi_MER(i) = ((T_MER(i)-T_MER(i-1)) / (psi1-1.0q0 ))
            end if
          end if
        
    end do
    
  end subroutine estimadores
  
  !-------------------------------------------------------------------------------
  !aplica cálculo de estimadores sobre solução sem MER
  !-------------------------------------------------------------------------------
  
  subroutine estimadores1 
 
    real*16 :: psi        ! razão de convergência (auxiliar)
    real*16 :: del        ! auxiliar   
      
    !Calculo dos estimadores de Eh 'sem MER'
   
     do i = 2, nm
     
          r21 = h(i-1) / h(i)    !razão de refino
          del = Th(i) - Th(i-1)  !auxiliar                          
        ! Calcula estimador Richardson com p0, U_Ri_p0
          U_Ri_p0(i) = del/((r21**pV(1))-1.0q0)             
        ! Calcula solução extrapolada com p0
          T_inf_p0(i) = Th(i) + U_Ri_p0(i)         
        ! Calcula estimador GCI com p0
          if (i==2) then
              U_GCI(i)= 3.0q0 * qabs ( U_Ri_p0(i) ) !Para 2 malhas, GCI com p0 e com Fs = 3!
          end if
          if(i>2) then                      
             if ( (qabs(psi_v(i))>1.0q0) .and. (psi_v(i)/=great) ) then
                   ! Calcula estimador Richardson com pU, U_Ri_pU
                      U_Ri_pU(i) = del/((r21**pU_h1(i))-1.0q0) 
                   ! Calcula solução extrapolada com pU
                      T_inf_pU(i) = Th(i) + U_Ri_pU(i)
                   ! Calcula solução numérica convergente
                      T_c(i) = (T_inf_p0(i)+T_inf_pU(i))/2.0q0 
                   ! Calcula ordem aparente e razão de convergência para estimativa da solução convergente    
                      if(i>4)then
                         del = T_c(i) - T_c(i-1)
                         if ( del /= 0.0q0 ) then
                            psi = (T_c(i-1)-T_c(i-2))/del
                            psi_c(i) = psi
                         else
                            psi = great
                         end if
                         if ( (psi/=0.0q0) .and. (psi/=great) ) then
                                r21 = h(i-1) / h(i)
                                pU_c1(i) = qlog10(qabs(psi))/qlog10(r21) 
                         else
                                pU_c1(i) = great 
                         end if
                       end if
                   ! Calcula estimador convergente
                       U_c(i) = qabs(U_Ri_p0(i)-U_Ri_pU(i))/2.0q0                  
                   ! Calcula estimador GCI para 3 ou mais malhas, p = min{p0,pU} e Fs=1,25 ou com p0 se pU não existir
                      if (pU_h(i) > 0 .and. pU_h(i) < great) then
                          p = qmin1 ( pL, pU_h(i) )
                          U_GCI(i)= 1.25q0 *  qabs (( Th(i) - Th(i-1) ) / ( (r21**p) - 1.0q0 ))
                      else
                          p = pL
                          U_GCI(i)= 1.25q0 *  qabs (( Th(i) - Th(i-1) ) / ( (r21**p) - 1.0q0 ))
                     end if            
             end if
           end if
           ! Calcula estimador multicoeficiente, com base em MER
             U_mc(i) = Tm_12(i,i-1) - Th(i)
    end do

  end subroutine estimadores1


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

end module calculos

