!> \mainpage{Derivate_Field_IcGrid2D_Manager モジュールのテストプログラム}
!!
!! \brief Derivate_Field_IcGrid2D_Manager モジュールのテストプログラム
!!
!! - 正二十面体格子上の物理場に対する微分演算を提供するモジュールのテスト  
!!  - 微分演算のテストでは, 水平発散, 水平勾配, 回転の鉛直成分, ラプラシアンの微分演算提供手続きをテストを行います. 
!!  - 数値誤差を定量的に見積もるために, 数値解と解析解から誤差ノルムを計算します.  
!!  - このテストプログラムは, MPI 並列に対応しています(実行時オプションを参照). 
!!
!! - 実行時オプションについて
!!  - MPI 並列を利用する場合は, 実行時オプションとして --enable-mpi を指定する必要があります. 
!! 
!! Authors::   Yuta Kawai
!! Copyright:: Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved.
!!
program diffOptr_test_main

  ! モジュール引用 ; Use statements
  !

#ifdef LIB_MPI
  use mpi
#endif    
  
  !* gtool *********************
  !*
  
  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP, &
    &                 TOKEN, &
    &                 STRING
  
  !
  !
  !
  use dc_message, only: &
    & MessageNotify
        
  !* IGMBaseLib *****************
  !*     
  
  ! 正二十面体格子データの管理
  ! Manager of Icosahedral grid data
  !
  use IcGrid2D_FVM_Manager
  
  !
  !
  !
  use Field_IcGrid2D_Manager

  !
  !
  !
  use IcGrid_ncWriter_mod
  
  !
  !
  !
  use IcGrid2D_FVM_Builder, only: &
    & construct_icosahedral_grid, &
    & GTYPE_STDGRID, GTYPE_STDGCGRID, GTYPE_SPRGCGRID

  !
  !
  !
  use IcGrid_ncGridLoader, only: load_gridData

  !
  !
  !
  use Derivate_Field_IcGrid2D_manager, only: &
    & Derivate_Field_IcGrid2D, &
    & Derivate_Field_IcGrid2D_Init, Derivate_Field_IcGrid2D_Final, &
    & divergence_op, gradient_op, vertical_curl_op, laplacian_op

  !
  !
  !
  use Field_IcGrid2D_error_norm_check
  
  !
  !
  !
  use test_diff_field, only: &
    & init_test_diff_field, final_test_diff_field,      &
    &  field_scalar, dfield_scalar, ana_dfield_scalar,  &
    &  field_vector, dfield_vector, ana_dfield_vector,  &
    &  create_field, calc_ana_dfield, &
    &  DIV_OPERATOR, GRADIENT_OPERATOR, ROT_OPERATOR, LAPLACIAN_OPERATOR 
  
  ! 宣言文 ; Declaration statements
  !
  implicit none
  
  !> \~japanese 正二十面格子データを管理する構造体 IcGrid_Field2D の変数. 
  !! \english 
  type(IcGrid2D_FVM) :: ic_grid
  
  !> \~japanese 
  !! \~english
  integer, parameter :: grid_type_id = GTYPE_SPRGCGRID
  
  !> \~japanese 
  !! \~english
  character(STRING) :: grid_filename = './SPR-GC-grid_glevel5.nc'
  
  !> \~japanese 
  !! \~english
  logical :: load_flag = .true.
 
  !> \~japanese 正二十面格子の水平分割レベル.  
  !! \~english  Horizontal division level of icosahedral grid(glevel). 　 
  integer, parameter :: glevel = 5
  
  !> \~japanese 正二十面格子の半径. 
  !! \~english Radius of icosahedral grid. 
  real(DP), parameter :: ic_radius = 1.0d0
 
  !> \~japanese モード数. 
  !! \~english  Number of mode. 
  integer :: m = 3
  
  !> \~japanese モード数. 
  !! \~english  Number of mode. 
  integer :: n = 3
  
  !> \~japanese 
  !! \~english
  type(Derivate_Field_IcGrid2D) :: diff_eval

  logical :: mpi_flag
  integer :: ierr

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! メインプログラム
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  ! 準備
  ! Preparation
  !
  
  ! MPI を使用するかを実行時引数で確認する. 
  !
  mpi_flag = check_MPIFLAG()
  
  if ( mpi_flag ) then

#ifdef LIB_MPI    
  	
  	! MPI の初期化 
	  ! Initialization of MPI
    call MPI_Init(ierr)
    
    call MessageNotify( 'M', 'diffOptr_test', &
      & 'MPI will be used in diffOptr_test program !', rank_mpi=0 )
#else
    call MessageNotify( 'E', 'diffOptr_test', &
      & 'diffOptr_test program is executed with MPI option,' & 
      & // 'but IGMBaseLib was not built with MPI Library')
#endif
  end if
   
  ! 正二十面体格子の生成を行う. 
  !
  
  ! 正二十面格子データ管理用の構造体を初期化する. 
  call IcGrid2D_FVM_Init(ic_grid, glevel, ic_radius, mpi_flag )
  
  ! load_flag が真の時は, 格子データをファイルを NetCDF ファイルから読み取る.  
  ! そうでなければ, 指定された正二十面格子タイプの格子を生成し, NetCDF ファイルに格子データを出力する. 
  if ( load_flag ) then
    call load_gridData( grid_filename, ic_grid )
  else
    call construct_grid( grid_filename, grid_type_id )
  end if

  ! 微分演算用変数の初期化を行う. 
  call Derivate_Field_IcGrid2D_Init(diff_eval, ic_grid)

  ! テスト用の物理場を生成するためのモジュールを初期化する. 
  call init_test_diff_field(ic_grid)
  
  ! 各種微分演算子の数値誤差の検証を行う. 
  !
  !
  
  ! 発散演算子のテスト. 
  call test_differential_operator(DIV_OPERATOR, m, n)
  ! 回転演算子のテスト. 
  call test_differential_operator(ROT_OPERATOR, m, n)
  ! 勾配演算子のテスト. 
  call test_differential_operator(GRADIENT_OPERATOR, m, n)
  ! ラプラス演算子のテスト. 
  call test_differential_operator(LAPLACIAN_OPERATOR, 1, 1)

  ! 後片付け
  !
  !
  call final_test_diff_field()
  
  call Derivate_Field_IcGrid2D_Final(diff_eval)

  if ( mpi_flag ) then
#ifdef LIB_MPI
    call MPI_Finalize(ierr)
#endif
  end if
  
contains

function check_MPIFLAG() result(mpi_flag)
  ! モジュール引用 ; Use statements
  !

  !
  !
  !
  use dc_string, only: StoA
  
  ! コマンドライン引数の解析
  ! Analization of command-line arguments
  !
  use dc_args, only: &
    & ARGS, &
    & DCArgsOpen, DCArgsOption, DCArgsClose
    
  !
  !  
  logical :: mpi_flag
  
  !
  !
  !
  type(ARGS) :: arg

  ! 実行文 ; Executable statements
  ! 
  
  call DCArgsOpen(arg)
  
  call DCArgsOption( &
    & arg, StoA('--enable-mpi'), mpi_flag )
    
  call DCArgsClose(arg)

end function check_MPIFLAG

subroutine construct_grid(file_name, grid_type_id)

  ! モジュール引用 ; Use statements
  !

  !
  !
  use IcGrid_ncWriter_mod, only: IcGrid_ncWriter, &
    & IcGrid_ncWriter_Init, open_ncFile, ncdef_GridData, end_ncdef, &
    & write_GridData, close_ncFile

  ! 宣言文; Declaration statements
  !
  character(*), intent(in) :: file_name
  integer, intent(in) :: grid_type_id

  !
  !
  type(IcGrid_ncWriter) :: writer

  ! 実行文 ; Executable statements
  ! 
 
  !
  call construct_icosahedral_grid(grid_type_id, ic_grid)

  !
  call IcGrid_ncWriter_Init(writer, ic_grid)
  call open_ncFile(writer, file_name)
  call ncdef_GridData(writer)
  call end_ncdef(writer)
  call write_GridData(writer)
  call close_ncFile(writer)

end subroutine construct_grid

subroutine test_differential_operator(diffOptrType, m, n)

  ! 宣言文; Declaration statements
  !
  integer, intent(in) :: diffOptrType 
  integer, intent(in) :: m, n
   
  !
  !
  character(len=10) :: str_diffOptrType
  integer :: dfield_dim
  type(Field_IcGrid2D), pointer :: ana_dfield, dfield
  
  ! 実行文 ; Executable statements
  ! 

  ! 場の微分を数値的に評価する
  !
  
  ! 
  select case(diffOptrType)
    case(DIV_OPERATOR)
      ! ベクトル場の作成
      call create_field( field_vector, diffOptrType, m, n )
      call paste_field2D_margin(field_vector)

      dfield_dim = 1
      call divergence_op(diff_eval, field_vector, dfield_scalar)
      dfield     => dfield_scalar
      ana_dfield => ana_dfield_scalar
       
    case(ROT_OPERATOR)

      ! ベクトル場の作成
      call create_field( field_vector, diffOptrType, m, n )
      call paste_field2D_margin(field_vector)

      dfield_dim = 1
      call vertical_curl_op(diff_eval, field_vector, dfield_scalar)
      dfield     => dfield_scalar
      ana_dfield => ana_dfield_scalar

    case(GRADIENT_OPERATOR)

      ! スカラー場の作成
      call create_field( field_scalar, diffOptrType, m, n )
      call paste_field2D_margin(field_scalar)

      dfield_dim = 3
      call gradient_op(diff_eval, field_scalar, dfield_vector)
      dfield     => dfield_vector
      ana_dfield => ana_dfield_vector

    case(LAPLACIAN_OPERATOR)

      ! スカラー場の作成(勾配演算子のテストに使った場(beta)を用いる). 

      call create_field( field_scalar, GRADIENT_OPERATOR, m, n )
      !field_scalar%val(:,:,:,:) = 10d0
      call paste_field2D_margin(field_scalar)

      dfield_dim = 1
      call laplacian_op(diff_eval, field_scalar, dfield_scalar)
      dfield     => dfield_scalar
      ana_dfield => ana_dfield_scalar
      ! dfield_scalar%val(:,:,:,:) = abs(dfield_scalar%val(:,:,:,:))
  end select

  ! 場の微分を解析的に評価する
  call calc_ana_dfield( ana_dfield, diffOptrType, dfield_dim, m, n )
  
  if (diffOptrType ==LAPLACIAN_OPERATOR)  then
    call write_out_fieldData(dfield, ana_dfield)
  end if

  ! 数値的に求めた場と解析に求めた場のデータから, 誤差ノルムを計算する. 
  write(str_diffOptrType, '(i1)') diffOptrType
  call check_norm(str_diffOptrType , ana_dfield, dfield)

end subroutine test_differential_operator

subroutine check_norm( &
	& str_optr_type, ana_dfield, dfield &
	& )

  ! モジュール引用 ; Use statements
  !

  !
  !
  use rcregions_norm_check

  ! 宣言文; Declaration statements
  !
  character(*), intent(in) :: str_optr_type
  type(Field_IcGrid2D), intent(in) :: ana_dfield
  type(Field_IcGrid2D), intent(in) :: dfield  

  ! 作業変数
  ! Work variables
  !
  real(DP) :: norm_1, norm_2, norm_3
  integer :: my_rank, ierr
  
  ! 実行文 ; Executable statements
  ! 

  norm_1 = numerical_error_norm1( ana_dfield, dfield) 
  norm_2 = numerical_error_norm2( ana_dfield, dfield)
  norm_3 = numerical_error_norminfinity( ana_dfield, dfield, .true.)

  if ( mpi_flag ) then
#ifdef LIB_MPI
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)
#endif
  else 
    my_rank = 0
  end if
  
  if ( my_rank == 0 ) then
    write(*,*) '===== Numerical error test: ', str_optr_type
    write(*,*) 'norm1:', norm_1
    write(*,*) 'norm2:', norm_2
    write(*,*) 'normInfinity:', norm_3
  end if
  
end subroutine check_norm

subroutine write_out_fieldData(dfield, anadfield)
  use IcGrid_ncWriter_mod
  
  type(Field_IcGrid2D), intent(in) :: dfield, anadfield
  
  integer :: dfieldID, anadfieldID, diffID
  type(IcGrid_ncWriter) :: writer
  type(Field_ICGrid2D) :: diff
 
  call Field_Icgrid2D_Init(diff, get_icgrid(dfield), 'diff', 1)
  diff%val = dfield%val - anadfield%val
  
  call IcGrid_ncWriter_Init(writer, get_icgrid(dfield), mpi_flag)
  call open_ncFile(writer, 'fieldData.nc')
  
  call ncdef_GridData(writer)  
  dfieldID = ncdef_FieldData(writer, dfield)
  anadfieldID = ncdef_FieldData(writer, anadfield)
  diffID = ncdef_FieldData(writer, diff)

  call ncdef_Simulation_Parameter(writer, 1d0, 0d0, 1d0)
  call end_ncdef(writer)
  
  call write_GridData(writer)
  call write_FieldData(writer, dfieldID, dfield)
  call write_FieldData(writer, anadfieldID, anadfield)
  call write_FieldData(writer, diffID, diff)

  call close_ncFile(writer)
  
  call IcGrid_ncWriter_Final(writer)
  
  call Field_IcGrid2D_Final(diff)
  
end subroutine write_out_fieldData

end program diffOptr_test_main
 
