module solvers
  ! Aims:
  !   1) Given the coefficients, sources and a initial guess, solve
  !      with the tdma method the system of algebraic equations like AX = B,
  !      where A has 5 non-null diagonals for pl and
  !                  9 non-null diagonals for u, v and T
  !   2) Calculate the L1 norm of the residual of the algebraic equations,
  !      for 5 and 9 non-null diagonals, considering only real volumes.
  !
  ! Latest version: 13 Jul 2011
  !
  ! Subroutines given by this module:
  !   1) norma_L1_5d
  !   2) norma_L1_9d
  !
  ! Subroutines given by module mtdma2d5
  !   3) set_tdma2d5_num_threads
  !   4) get_tdma2d_num_threads()
  !   5) tdma2d5
  !
  ! Subroutines given by module mtdma2d9
  !   6) set_tdma2d9_num_threads
  !   7) get_tdma2d9_num_threads()
  !   8) tdma2d9
  !

  use mtdma2d5
  use mtdma2d9

  implicit none

contains

  subroutine norm_l1_5d( nx, ny, var, b, a, norm)
    implicit none
    integer, intent(in) :: nx   ! Number of volumes in csi direction (real+fictitious)
    integer, intent(in) :: ny   ! Number of volumes in eta direction (real+fictitious)
    real(8), dimension (nx*ny),   intent(in) :: var ! Unknown variable
    real(8), dimension (nx*ny),   intent(in) :: b   ! Source vector of the linear system
    real(8), dimension (nx*ny,5), intent(in) :: a   ! Coefficients of the linear system
    real(8), intent(inout) :: norm

    ! Auxiliary variables
    integer :: i, j, np, nps, npn, npw, npe

    ! Norm is calculated taking into account only real volumes

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

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

          norm = norm + abs(        &
               + a(np,1) * var(nps) &
               + a(np,2) * var(npw) &
               + a(np,3) * var(np ) &
               + a(np,4) * var(npe) &
               + a(np,5) * var(npn) &
               - b(np) )

       end do
    end do

  end subroutine norm_l1_5d


  subroutine norm_l1_9d( nx, ny, var, b, a, norm)
    implicit none
    integer, intent(in) :: nx   ! Number of volumes in csi direction (real+fictitious)
    integer, intent(in) :: ny   ! Number of volumes in eta direction (real+fictitious)
    real(8), dimension (nx*ny),   intent(in) :: var ! Unknown variable
    real(8), dimension (nx*ny),   intent(in) :: b   ! Source vector of the linear system
    real(8), dimension (nx*ny,9), intent(in) :: a   ! Coefficients of the linear system
    real(8), intent(inout) :: norm

    ! Auxiliary variables
    integer :: i, j, np, nps, npn, npw, npe, npsw, npse, npnw, npne

    ! Norm is calculated taking into account only real volumes

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

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

          norm = norm + abs(        &
               + a(np,1) * var(npsw) &
               + a(np,2) * var(nps ) &
               + a(np,3) * var(npse) &
               + a(np,4) * var(npw ) &
               + a(np,5) * var(np  ) &
               + a(np,6) * var(npe ) &
               + a(np,7) * var(npnw) &
               + a(np,8) * var(npn ) &
               + a(np,9) * var(npne) &
               - b(np) )

       end do
    end do

  end subroutine norm_l1_9d


    !> \brief Calculates the coefficient for convergence criteria of a 9-diagonal matrix.
    !! This coefficient is calculated only for real volumes.
    real(8) function get_maxc_9d(nx, ny, a)
        implicit none
        integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
        real(8), dimension (nx*ny,9), intent(in) :: a   !< Coefficients of the linear system

        ! Inner variables

        integer :: i, j, np

        real(8) :: c, maxc

        ! Looking for the maxval of c = ( sum_nb | A_nb | ) / | A_P |

        maxc = 0.d0

        do i = 2, nx-1

            do j = 2, ny-1

                np   = nx * (j-1) + i

                c = ( sum(abs(a(np,1:4))) + sum(abs(a(np,6:9))) ) / abs( a(np,5) )

                if ( maxc < c ) maxc = c

            end do

        end do

        get_maxc_9d = maxc

    end function


    !> \brief Calculates the coefficient for convergence criteria of a 5-diagonal matrix.
    !! This coefficient is calculated only for real volumes.
    real(8) function get_maxc_5d(nx, ny, a)
        implicit none
        integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
        integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
        real(8), dimension (nx*ny,5), intent(in) :: a   !< Coefficients of the linear system

        ! Inner variables

        integer :: i, j, np

        real(8) :: c, maxc

        ! Looking for the maxval of c = ( sum_nb | A_nb | ) / | A_P |

        maxc = 0.d0

        do i = 2, nx-1

            do j = 2, ny-1

                np   = nx * (j-1) + i

                c = ( sum(abs(a(np,1:2))) + sum(abs(a(np,4:5))) ) / abs( a(np,3) )

                if ( maxc < c ) maxc = c

            end do

        end do

        get_maxc_5d = maxc

    end function
end module solvers
