!> \~japanese
!! \brief 数値計算に伴う誤差を評価するためのサブルーチンを提供するモジュール.
!! \~english
!! \brief 数値計算に伴う誤差を評価するためのサブルーチンを提供するモジュール.
!!
!! \~
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module rcregions_norm_check
  
  ! モジュール引用 ; Use statement
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision. 
  
  ! 宣言文 ; Declaration statement
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public calc_norm1, calc_norm2, calc_norm3

  ! 公開変数
  ! Public variables
  !
  integer, parameter :: RC_REGIONS_NUM = 10

contains

!
!
!
function calc_norm1(true_field, num_field) result(normVal)
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: true_field(:,:,:,:)
  real(DP), intent(in) :: num_field(:,:,:,:)
  real(DP) normVal

  ! 作業変数
  ! Work varaibles
  !
  real(DP) fieldT_sum
  real(DP) fieldError_sum
  integer :: rcID, i,j
  integer :: EMax
  
  ! 実行文 ; Executable statements
  !

  EMax = size(true_field, 2)
  fieldT_sum = 0.0d0
  fieldError_sum = 0.0d0

  do rcID=1, RC_REGIONS_NUM
    do j=1,EMax-1
      do i=2,EMax
        fieldError_sum = fieldError_sum + abs_Val(true_field(:, i,j,rcID) - num_field(:, i,j,rcID))
        fieldT_sum = fieldT_sum + abs_Val(true_field(:, i,j,rcID))
        
      end do
    end do
  end do
  
  fieldError_sum = fieldError_sum + abs_Val(true_field(:,1,1,1) - num_field(:,1,1,1))
  fieldT_sum = fieldT_sum + abs_Val(true_field(:,1,1,1))
    
  fieldError_sum = fieldError_sum + abs_Val(true_field(:,EMax,EMax,6) - num_field(:, EMax,EMax,6))
  fieldT_sum = fieldT_sum + abs_Val(true_field(:, EMax,EMax,6))
  
  normVal = fieldError_sum / fieldT_sum

end function calc_norm1

!
!>
!!
!!
function calc_norm2(true_field, num_field) result(normVal)
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: true_field(:,:,:,:)
  real(DP), intent(in) :: num_field(:,:,:,:)
  real(DP) normVal

  ! 作業変数
  ! Work varaibles
  !
  real(DP) fieldT_sum
  real(DP) fieldError_sum
  integer :: rcID, i,j
  integer :: EMax


  ! 実行文 ; Executable statements
  !

  EMax = size(true_field, 2)
  fieldT_sum = 0.0d0
  fieldError_sum = 0.0d0

  do rcID=1, RC_REGIONS_NUM
    do j=1,EMax-1
      do i=2,EMax
        fieldError_sum = fieldError_sum + abs_Val(true_field(:, i,j,rcID) - num_field(:, i,j,rcID))**2
        fieldT_sum = fieldT_sum + abs_Val(true_field(:, i,j,rcID))**2
      end do
    end do
  end do

  fieldError_sum = fieldError_sum + abs_Val(true_field(:, 1,1,1) - num_field(:, 1,1,1))**2
  fieldT_sum = fieldT_sum + abs_Val(true_field(:, 1,1,1))**2
  
  fieldError_sum = fieldError_sum + abs_Val(true_field(:,EMax,EMax,6) - num_field(:,EMax,EMax,6))**2
  fieldT_sum = fieldT_sum + abs_Val(true_field(:,EMax,EMax,6))**2
  
  normVal = dsqrt(fieldError_sum) / dsqrt(fieldT_sum)
    
end function calc_norm2

!
!>
!!
function calc_norm3(true_field, num_field, output_info_flag ) result(normVal)
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: true_field(:,:,:,:)
  real(DP), intent(in) :: num_field(:,:,:,:)
  logical, intent(in), optional :: output_info_flag
  real(DP) normVal

  ! 作業変数
  ! Work varaibles
  !
  real(DP) :: fieldTMax
  real(DP) :: fieldErrorMax
  real(DP) absfieldError, absfieldT 
  integer :: rcID, i,j

  integer :: max_rcID, max_i, max_j
  integer :: EMax

  ! 実行文 ; Executable statements
  !

  EMax = size(true_field, 2)
  fieldTMax = 0.0d0
  fieldErrorMax = 0.0d0
  max_rcID = 1; max_i = 1; max_j = 1
 
  do rcID=1, RC_REGIONS_NUM
    do j=1,EMax
      do i=1,EMax
        
        absfieldError = abs_Val(true_field(:, i,j,rcID) - num_field(:, i,j,rcID))
        absfieldT = abs_Val(true_field(:, i,j,rcID))
        
        if( absfieldError >= fieldErrorMax ) then
          fieldErrorMax = absfieldError; 
          max_rcID=rcID; max_i=i; max_j=j;
        end if
        
        if( absfieldT > fieldTMax ) then
          fieldTMax = absfieldT
         ! write(*,*) rcID, i,j, absfieldError, absfieldT
        end if
        
      end do
    end do
  end do

  normVal = fieldErrorMax / fieldTMax

  if ( present(output_info_flag) .and. output_info_flag ) then
     write(*,*) '= norm3 maximum error info:'
     write(*,'(a,3(i6))') '* location(rcID, i, j):', max_rcID, max_i, max_j
     write(*,*) '* Max Val:', fieldTMax
     write(*,*) '* True value:', true_field(:,max_i,max_j, max_rcID)
     write(*,*) '* Calculated value:', num_field(:,max_i,max_j, max_rcID)

  end if
end function  calc_norm3

function abs_Val(errorVal) result(absVal)
  !
  !
  !
  real(DP), intent(in)  :: errorVal(:) 
  real(DP)absVal

  integer :: valDim
  valDim = size(errorVal, 1)

  if( valDim == 1 )then
    absVal = dabs(errorVal(1))
  else if( valDim == 3 ) then
    absVal = dsqrt(dot_product(errorVal, errorVal))
  else
    absVal = 0.0d0
    write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!'
  end if

end function abs_Val

end module rcregions_norm_check
