!> \~japanese
!! \brief 正二十面体格子上で水平 2 次元的に分布する物理場に作用する空間微分演算子を評価するための手続きと
!! その微分演算に必要な情報を管理するための構造型を提供するモジュール.
!!
!! \par Derivate_Field_IcGrid2D_Manager モジュールについて
!!   - 本モジュールは, 正二十面格子上で水平 2 次元的に分布する物理場に対して, 空間微分演算を行うための手続きを提供する. 
!!     - 空間微分演算として, 球面上の水平勾配, 水平発散, 回転の鉛直成分, 水平ラプラシアンを利用できる.  
!!
!! \~english
!! \brief This module provides some subroutine to evaluate the spatial differential operators acting on the horizontally distoributed physical field,
!! and a derived type to manage the information necessary for the differential operations.
!!
!! \par About Derivate_Field_IcGrid2D_Manager
!!   - 本モジュールは, 正二十面格子上で水平 2 次元的に分布する物理場に対して空間微分演算を行うための手続きを提供する. 
!!     - 空間微分演算として, 球面上の水平勾配, 水平発散, 回転の鉛直成分, 水平ラプラシアンを利用できる.  
!!
!! \~
!! \since IGMBaseLib 1.0
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module Derivate_Field_IcGrid2D_Manager

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

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision

  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: &
    & vec_normalize, &
    & cross, dot

  ! 球面三角法
  ! Spherical trigonometry
  !
  use igmcore_spherical_trigonometry, only: &
    & spherical_tri_area, geodesic_arc_length

  ! 一つの格子点に対する微分演算子の提供
  ! Providing with the differential operation for a single grid-point
  !
  use FVM_HDiff_operator, only: &
    & eval_div_operator, eval_rot_operator, eval_grad_operator, &
    & FVM_HDerivate_CVInfo

  ! 正二十面体格子データの管理
  ! Manager of icosahedral grid data
  !
   use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, IcGrid2D_FVM_CV, &
    & paste_margin_width, generate_CV5_GPindex, generate_CV6_GPindex, &
    & get_IdMin, get_IdMax, &
    & get_EffSize_Min, get_EffSize_Max, get_IcRadius, &
    & get_LocalProcRectInfo, &
    & RC_REGIONS_NUM

  ! 正二十面格子における広域矩形領域のサブドメインである局所矩形領域(各 MPI プロセスと対応)の情報管理
  ! Manager of information about a local rectangle region (corresponding to each MPI process)
  !
  use LocalProcRect_Manager, only: &
    & LocalProcRectInfo, &
    & get_InteriorRegion_IndexInfo
    
  ! 正二十面体格子上の物理場データの管理
  ! Manager of the physical field data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & get_field2D_rank, get_icgrid, &
    & paste_field2d_margin

  ! OpenMP
  !
  !
  use omp_lib

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開変数
  ! Public variables
  !

  !
  !> \~japanese
  !! \brief 正二十面体格子上で水平 2 次元的に分布する物理場の水平微分(勾配, 発散, 回転)を評価する際に必要な情報を管理する構造型.
  !! \~english
  !! \brief The derived type to manage the information necessary for evaluating the horizontal differential operators(divegence, gradient, curl)
  !!acting on the horizontally distirbuted physical field.
  !!
  type, public :: Derivate_Field_IcGrid2D

    !> \~japanese 構造体 IcGrid2D_FVM の変数へのポインタ.
    !! \~english The pointer to derived type IcGrid2D_FVM.
    type(IcGrid2D_FVM), pointer :: icgrid
    ! f2003: type(IcGrid2D_FVM), pointer, private :: icgrid


    !> \~japanese
    !! コントロールボリューム上の物理場の値を求める際に重みとして必要となる,
    !! コントロールボリュームの頂点と近傍の格子点で作られる 3 つの三角形の面積を保持する配列.
    !! \~english
    !! コントロールボリューム上の物理場の値を求める際に重みとして必要となる,
    !! コントロールボリュームの頂点と近傍の格子点で作られる 3 つの三角形の面積を保持する配列.
    real(DP), pointer :: rcs_areaGPP(:,:,:,:,:)
    ! f2003: real(DP), allocatable, private :: rcs_areaGPP(:,:,:,:,:)

    !> \~japanese 北半球側の矩形領域において, コントロールボリュームの頂点を取り囲む近傍の 3 つの格子点のインデックスを保持する配列.
    !! \~english 北半球側の矩形領域において, コントロールボリュームの頂点を取り囲む近傍の 3 つの格子点のインデックスを保持する配列.
    integer, pointer :: rcU_CVGPindex(:,:,:,:,:)
    ! f2003: integer, allocatable, private :: rcU_CVGPindex(:,:,:,:,:)

    !> \~japanese 南半球側の矩形領域において, コントロールボリュームの頂点を取り囲む近傍の 3 つの格子点のインデックスを保持する配列.
    !! \~english 南半球側の矩形領域において, コントロールボリュームの頂点を取り囲む近傍の 3 つの格子点のインデックスを保持する配列.
    integer, pointer :: rcL_CVGPindex(:,:,:,:,:)
    ! integer, allocatable, private :: rcL_CVGPindex(:,:,:,:,:)

    type(FVM_HDerivate_CVInfo), pointer :: CVInfo(:,:,:)

  end type Derivate_Field_IcGrid2D

  ! 公開メンバ関数
  ! Public memeber function
  !
  public :: Derivate_Field_IcGrid2D_Init, Derivate_Field_IcGrid2D_Final
  public :: divergence_op, vertical_curl_op, gradient_op, laplacian_op

contains

!
!> \~japanese
!! \brief 構造型 Derivate_Field_IcGrid2D の変数を初期化する.
!!
!! @param[in,out] self 構造型 Derivate_Field_IcGrid2D の変数.
!! @param[in]     icgrid 正二十面格子を管理する構造型 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Initializes the variable of derived type Derivate_Field_IcGrid2D.
!!
!! @param[in,out] self The variable of derived type Derivate_Field_IcGrid2D.
!! @param[in]     icgrid The variable of derived type IcGrid2D_FVM which manages icosahedral grid.
!!
subroutine Derivate_Field_IcGrid2D_Init( &
  & self,              &  ! (inout)
  & icgrid             &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self
  type(IcGrid2D_FVM), intent(in), target :: icgrid

  ! 作業変数
  ! Work variables
  !
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax

  ! 実行文 ; Declaration statements
  !

  self%icgrid => icgrid
  lcRCInfo => get_LocalProcRectInfo(self%icgrid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  ! Derivate_Field_IcGrid2D クラスが作業変数として保持する配列を確保する.
  allocate(self%rcs_areaGPP(3, 6, iEMin:iEMax, jEMin:jEMax, rcIDMin:rcIDMax))
  allocate(self%rcU_CVGPindex(6,3,2, iEMin:iEMax, jEMin:jEMax ))
  allocate(self%rcL_CVGPindex(6,3,2, iEMin:iEMax, jEMin:jEMax ))
  allocate(self%CVInfo(iEMin:iEMax, jEMin:jEMax, rcIDMin:rcIDMax))

  ! コントールボリュームの頂点座標を計算するために必要な格子点のインデックスを, あらかじめ計算しておく.
  call calc_CVGPindex(self)
  ! コントールボリュームの頂点上での物理場の値を計算する際に必要となる, 面積データをあらかじめ計算しておく.
  call calc_GPP_area_weight(self)

  ! 各種水平微分演算のために用いられるコントロールボリュームの情報を計算する. 
  call set_HDerivate_CVInfo(self)

end subroutine Derivate_Field_IcGrid2D_Init

!
!> \~japanese
!! \brief 構造型 Derivate_Field_IcGrid2D の変数を最終化する.
!!
!! @param[in,out] self 構造型 Derivate_Field_IcGrid2D の変数.
!!
!! \~english
!! \brief Finalizes a variable of derived type Derivate_Field_IcGrid2D.
!!
!! @param[in,out] self A variable of derived type Derivate_Field_IcGrid2D.
!!
subroutine Derivate_Field_IcGrid2D_Final( &
  & self &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self

  ! 実行文 ; Declaration statements
  !

  deallocate(self%rcs_areaGPP)
  deallocate(self%rcU_CVGPindex)
  deallocate(self%rcL_CVGPindex)
  deallocate(self%CVInfo)

end subroutine Derivate_Field_IcGrid2D_Final

!
!> \~japanese
!! \brief 正二十面体格子上の物理場の水平発散を計算する.
!!
!! @param[in,out] self              構造型 Derivate_Field_IcGrid2D の変数.
!! @param[in]     vector_field      発散演算子を作用させる物理場データを管理する構造型 Field_IcGrid2D の変数.
!! @param[in,out] ret_scalar_field  計算結果の水平発散場を保存する構造型 Field_IcGrid2D の変数.
!!
!! \~english
!! \brief Calculates the horizontal divergence of the physical field on icosahedral grid.
!!
!! @param[in,out] self              The variable of derived type Derivate_Field_IcGrid2D.
!! @param[in]     vector_field      The variable of derived type Field_IcGrid2D to manage the physical field which is acted by the divergence operator.
!! @param[in,out] ret_scalar_field  The variable of derived type Field_IcGrid2D to save the resulting horizotal divergence field.
!!
subroutine divergence_op( &
  & self,                &  ! (inout)
  & vector_field,        &  ! (in)
  & ret_scalar_field     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: vector_field
  type(Field_IcGrid2D), intent(inout) :: ret_scalar_field
  
  ! 作業変数
  ! Work variables
  !

  ! 実行文 ; Executable statements
  !

  call eval_diffOptr( &
    & self=self, &
    & diffOptr_fnc=eval_div_operator, & 
    & optr_For_CV5Val=calc_CV5_val, optr_For_CV6Val=calc_CV6_val, &
    & in_field2D=vector_field, val_dim=3,                         &  ! (in)
    & out_dfield2D=ret_scalar_field, ret_dim=1                    &  
    & )                                     

end subroutine divergence_op

!
!> \~japanese
!! \brief 正二十面体格子上に分布する物理場の回転の鉛直成分を計算する.
!!
!! @param[in,out] self              構造型 Derivate_Field_IcGrid2D の変数.
!! @param[in]     vector_field      回転演算子を作用させる物理場データを保持する構造型 Field_IcGrid2D の変数.
!! @param[in,out] ret_scalar_field  計算結果の回転の鉛直成分を保存する構造型 Field_IcGrid2D の変数.
!!
!! \~english
!! \brief Calculates the vertical component of curl of the physical field on icosahedral grid.
!!
!! @param[in,out] self              The variable of derived type Derivate_Field_IcGrid2D.
!! @param[in]     vector_field      The variable of derived type Field_IcGrid2D to manage the physical field which is acted by the (vertical) curl operator.
!! @param[in,out] ret_scalar_field  The variable of derived type Field_IcGrid2D to save the resulting component of curl of the field.
!!
subroutine vertical_curl_op( &
  & self,                &  ! (inout)
  & vector_field,        &  ! (in)
  & ret_scalar_field     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: vector_field
  type(Field_IcGrid2D), intent(inout) :: ret_scalar_field

  ! 作業変数
  ! Work variables
  !

  ! 実行文 ; Executable statements
  !
  
  call eval_diffOptr( &
    & self=self, &
    & diffOptr_fnc=eval_rot_operator, & 
    & optr_For_CV5Val=calc_CV5_val, optr_For_CV6Val=calc_CV6_val, &
    & in_field2D=vector_field, val_dim=3,    &  ! (in)
    & out_dfield2D=ret_scalar_field,  ret_dim=1                              &
     & )                                     

end subroutine vertical_curl_op

!
!> \~japanese
!! \brief 正二十面体格子上に分布する物理場の水平勾配を計算する.
!!
!! @param[in,out] self              構造型 Derivate_Field_IcGrid2D の変数.
!! @param[in]     vector_field      水平勾配演算子を作用させる物理場データを保持する構造型 Field_IcGrid2D の変数.
!! @param[in,out] ret_scalar_field  水平勾配を保存する構造型 Field_IcGrid2D の変数.
!!
!! \~english
!! \brief Calculates the horizontal gradient of the physical field on icosahedral grid.
!!
!! @param[in,out] self              The variable of derived type Derivate_Field_IcGrid2D.
!! @param[in]     vector_field      The variable of derived type Field_IcGrid2D to manage the physical field which is acted by the gradient operator.
!! @param[in,out] ret_scalar_field  The variable of derived type Field_IcGrid2D to save the resulting horizotal gradient field.
!!
subroutine gradient_op( &
  & self,                &  ! (inout)
  & scalar_field,        &  ! (in)
  & ret_vector_field     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: scalar_field
  type(Field_IcGrid2D), intent(inout) :: ret_vector_field
  
  ! 作業変数
  ! Work variables
  !

  ! 実行文 ; Executable statements
  !

  call eval_diffOptr( &
    & self=self, &
    & diffOptr_fnc=eval_grad_operator, & 
    & optr_For_CV5Val=calc_CV5_val, optr_For_CV6Val=calc_CV6_val, &
    & in_field2D=scalar_field, val_dim=1,     &  ! (in)
    & out_dfield2D=ret_vector_field, ret_dim=3                                 &
    & )                                     

end subroutine gradient_op

subroutine laplacian_op( &
  & self,                &  ! (inout)
  & in_field,        &  ! (in)
  & ret_field     &  ! (inout)
  & )


  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: in_field
  type(Field_IcGrid2D), intent(inout) :: ret_field
  
  ! 作業変数
  ! Work variables
  !
  integer :: fieldRank, rankID
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid2D) :: field_comp, ret_field_comp
  
  ! 実行文 ; Executable statements
  !
  fieldRank = get_field2D_rank(in_field)
  icgrid2D => get_icgrid(in_field)
  
  call Field_IcGrid2D_Init(field_comp, icgrid2D, 'field_comp', 1)
  call Field_IcGrid2D_Init(ret_field_comp, icgrid2D, 'ret_field_comp', 1)
  

  do rankID =1, fieldRank
    field_comp%val(1,:,:,:) = in_field%val(rankID,:,:,:)
    call paste_field2d_margin(field_comp)
    
    call eval_diffOptr( &
      & self=self, &
      & diffOptr_fnc=eval_div_operator, & 
      & optr_For_CV5Val=calc_CV5_gradVec, optr_For_CV6Val=calc_CV6_gradVec, &
      & in_field2D=field_comp, cvVal_dim=3, val_dim=1,                      &  ! (in)
      & out_dfield2D=ret_field_comp, ret_dim=1                              &
       & )
    
    ret_field%val(rankID, :,:,:) = ret_field_comp%val(1, :,:,:)
                                            
  end do
  
  call Field_IcGrid2D_Final(field_comp)
  call Field_IcGrid2D_Final(ret_field_comp)
  
end subroutine laplacian_op

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
!
! モジュール外非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!> \~japanese
!! \brief 引数で与えられた物理場に対して微分演算を施す.
!!
!! 物理場の発散・回転・勾配のいずれかを計算する.
!!
!! @param[in,out] self        構造型 Derivate_Field_IcGrid2D の変数.
!! @param[in]     diffOptr    各微分演算を実装している呼び出し規則が共通なサブルーチン.
!! @param[in]     rcs_GP_val  微分演算子を作用させる正二十面体格子上の物理場データを保持する配列.
!! @param[in]     val_dim　    微分演算を作用させる物理場の次元.
!! @param[in,out] ret_rcs_val 微分演算の結果を保存する配列.
!! @param[in]     ret_dim     微分演算子の結果生じる物理場の次元.
!! @param[in]     idMin       物理場データを保持する配列 rcs_GP_val の第 2(or 3)インデックスの下限.
!!
!! \~english
!! \brief Calculates the differenatial of the specified physical field.
!!
!! One of the differential operation(divergence, curl and gradient) on the physical field will be performed.
!!
!! @param[in,out] self        The variable of derived type Derivate_Field_IcGrid2D.
!! @param[in]     diffOptr    各微分演算を実装している呼び出し規則が共通なサブルーチン.
!! @param[in]     rcs_GP_val  The array containing the physical field data on icosahedral grid which is acted by the differential operator.
!! @param[in]     val_dim　   The dimension of the physical field acted by differential operators.
!! @param[in,out] ret_rcs_val The array in which the result of differential operations is stored.
!! @param[in]     ret_dim     The dimension of the physcal field generated by the differential operation.
!! @param[in]     idMin       The lower bound index of array 'rcs_GP_val' in the second (or third) dimension. This array contains the physical field data.
!!
!!
subroutine eval_diffOptr( &
  & self,                                           &  ! (inout)
  & diffOptr_fnc, optr_For_CV5Val, optr_For_CV6Val, &  ! (in) 
  & in_field2D, val_dim, CVval_dim,                 &  ! (in)
  & out_dfield2D,                       &  ! (inout)
  & ret_dim                             &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: in_field2D
  integer, intent(in) :: val_dim
  integer, optional :: CVval_dim 
  type(Field_IcGrid2D), intent(inout) :: out_dfield2D    
  integer, intent(in) :: ret_dim

  interface
    !
    function diffOptr_fnc( &
      & val_GP0, val_CV, GP0, CVInfo,       &  ! (in)
      & CVArea, val_dim, ret_dim            &  ! (in)
      & ) result(val)
      
      use dc_types, only: DP
      use FVM_HDiff_operator, only: FVM_HDerivate_CVInfo
      implicit none

      ! 宣言文 ; Declaration statements
      !

      integer, intent(in) :: val_dim
      integer, intent(in) :: ret_dim
      type(FVM_HDerivate_CVInfo), intent(in) :: CVInfo
      real(DP), intent(in) :: val_CV(val_dim, CVInfo%CV_num)
      real(DP), intent(in) :: val_GP0(val_dim)
      real(DP), intent(in) :: GP0(3)
      real(DP), intent(in) :: CVArea
      real(DP) val(ret_dim)

    end function diffOptr_fnc
    
    !!
   subroutine optr_For_CV5Val( &
	   & self, &
	   & CV5_val, & 
	   & rcs_GP_val, dataNearNPole, dataNearSPole, &
	   & rcID, gp_i, gp_j, cvVal_dim, val_dim, &
	   & iIdMin, jIdMin, rcIDMin, glEMin, glEMax   & 
	    & )

      use dc_types, only: DP
      use Field_IcGrid2D_Manager
      import
      implicit none

       ! 宣言文 ; Declaration statement
       !
      type(Derivate_Field_IcGrid2D), intent(inout) :: self
      integer, intent(in) :: cvVal_dim, val_dim, iIdMin, jIdMin, rcIDMin, glEMin, glEMax
      real(DP), intent(out) :: CV5_val(val_dim, 5)
      real(DP), intent(in) :: rcs_GP_val(:,iIdMin:,jIdMin:,rcIDMin:)
      real(DP), intent(in) :: dataNearNPole(:,glEMin:,glEMin:,1:)
      real(DP), intent(in) :: dataNearSPole(:,glEMax-1:,glEMax-1:,6:) 
      integer, intent(in) :: rcID, gp_i, gp_j
    
    end subroutine optr_For_CV5Val
    
    subroutine optr_For_CV6Val( & 
	    & self, &
	    & CV6_val, &
	    & CV6GPVal9, CV6Area3, &
	    & rcID, gp_i, gp_j, cvVal_dim, val_dim &
	     & )
	    
      use dc_types, only: DP
      use Field_IcGrid2D_Manager
      import
      implicit none
      
       ! 宣言文 ; Declaration statement
       !
      
      type(Derivate_Field_IcGrid2D), intent(inout) :: self
      integer, intent(in) :: cvVal_dim, val_dim, gp_i, gp_j
      real(DP), intent(out) :: CV6_val(cvVal_dim, 6)
      real(DP), intent(in) :: CV6GPVal9(val_dim, gp_i-1:gp_i+1, gp_j-1:gp_j+1)
      real(DP), intent(in) :: CV6Area3(3,6)
      integer, intent(in) :: rcID
    end subroutine optr_For_CV6Val
  
  end interface

  
      
  ! 作業変数
  ! Work variables
  !
  real(DP) :: ic_radius
  real(DP), allocatable :: CV_val(:,:,:,:,:)
  real(DP) :: GPValCandiate(3,3,3)
  
  integer :: rcID, i,j
  integer :: CV_num
  integer :: CVValRank
  logical :: singularPtFlag = .false.
!$ integer :: chunk_size
  type(LocalProcRectInfo), pointer :: lcRCInfo
  
  integer :: glEMin, glEMax
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax, cvID, k
  real(DP) :: a(3,3,5)
  
  ! 実行文 ; Declaration statements
  !
  ic_radius = get_IcRadius(self%icgrid)
  glEMin = get_EffSize_Min(self%icgrid)
  glEMax = get_EffSize_Max(self%icgrid)
 
  lcRCInfo => get_LocalProcRectInfo(self%icgrid)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                          ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax )  ! (out) 

  ! コントロールボリュームの頂点上の物理場の値を保持するための配列のメモリを確保する.
  if ( .not. present(cvVal_dim) ) then
    CVValRank = val_dim
  else
    CVValRank = cvVal_dim
  end if
  
  allocate(CV_Val(1:CVValRank, 6, iEMin:iEMax, jEMin:jEMax, rcIDMin:rcIDMax))

  do rcID=rcIDMin, rcIDMax
   
   !$omp parallel
   !$ chunk_size = RC_REGIONS_NUM / omp_get_num_threads()
   !$omp do private(j, i, CV_num, singularPtFlag, GPValCandiate) schedule(static, chunk_size)
    do j=jEMin, jEMax
      do i=iEMin, iEMax

        if ( ( i==glEMin .and. j==glEMin ) .or. ( i==glEMax .and. j==glEMin ) .or. &
             ( i==glEMin .and. j==glEMax ) .or. ( i==glEMax .and. j==glEMax ) ) then
          ! 特異格子点のコントールボリュームの頂点数は 5 個
          CV_num = 5
          singularPtFlag = .true.

        else
          ! 通常ののコントールボリュームの頂点数は 6 個
          CV_num = 6
         singularPtFlag = .false.
        end if

        if ( singularPtFlag ) then
           ! 格子点(rcID,i,j) のコントールボリュームの頂点上における物理場の値を計算する.
          call optr_For_CV5Val( &
           & self, &
           & CV_val(1:CVValRank, 1:5, i,j,rcID),                   &  ! (inout)
           & in_field2D%val(1:val_dim,:,:,:),                    &  ! (in)
           & in_field2D%dataNearNPole, in_field2D%dataNearSPole, &  ! (in)
           & rcID, i, j, CVValRank, val_dim,                     &  ! (in)
           & lcRCInfo%iIdMin, lcRCInfo%jIdMin,                   &  ! (in) 
           & rcIDMin, glEMin, glEMax                             &  ! (in)
             & ) 

        else
  
           ! 格子点(rcID,i,j) のコントールボリュームの頂点上における物理場の値を計算する.
          GPValCandiate(1:val_dim, :,:) = in_field2D%val(1:val_dim, i-1:i+1, j-1:j+1, rcID)
          
          call optr_For_CV6Val( &
           & self, &
           & CV_val(1:CVValRank, 1:6, i,j,rcID),  &     ! (inout)
           & GPValCandiate(1:val_dim, 1:3, 1:3),  &     ! (in)
           & self%rcs_areaGPP(1:3,1:6, i,j,rcID), &     ! (in)
           & rcID, i,j, CVValRank, val_dim )            ! (in)
 
        end if

        ! 動的に接続された微分演算オブジェクトに対して微分演算を実行し, 格子点(rcID,i,j)上の微分値を計算する.
        ! また, その結果を保存する.
 
       out_dfield2D%val(1:ret_dim, i,j,rcID) = &       
         &    diffOptr_fnc( &
               & val_GP0=in_field2D%val(1:val_dim, i,j,rcID),     &  ! (in)
               & val_CV=CV_val(1:CVValRank, 1:CV_num, i,j,rcID),    &  ! (in)
               & GP0=self%icgrid%rcs_AGrid(1:3, i,j,rcID),        &  ! (in)
               & CVInfo=self%CVInfo(i,j,rcID),                    &  ! (in)
               & CVArea=self%icgrid%rcs_CVInfo(i,j,rcID)%cv_area, &  ! (in)
               & val_dim=val_dim, ret_dim=ret_dim )                  ! (in)

      end do
    end do  !$omp end do
  
    !$omp end parallel

  end do

end subroutine eval_diffOptr

!
!> \~japanese
!! \brief コントロールボリュームの頂点上の物理場の値を, それを取り囲む近傍の格子点の値を使って補間する.
!!
!!
!! @param[in,out] self   構造型 Derivate_Field_IcGrid2D の変数.
!! @param[in,out] CV_val コントロールボリュームの頂点上の物理場の値を格納する配列.
!! @param[in]     rcs_GP_val 格子点上の物理場の値を格納する配列.
!! @param[in]     rcID    矩形領域のインデックス.
!! @param[in]     gp_i    各矩形領域における水平第一成分のインデックス.
!! @param[in]     gp_j    各矩形領域における水平第二成分のインデックス.
!! @param[in]     val_dim 物理場の次元.
!! @param[in]     CV_num  コントロールボリュームの頂点数.
!! @param[in]     idMin 構造型 Derivate_Field_IcGrid2D_Manager の成分 rcs_AGrid の袖領域を含めた配列の最小インデックス.
!! @param[in]     EMin 構造型 Derivate_Field_IcGrid2D_Manager の成分 rcs_AGrid の袖領域を含めない配列の最小インデックス.
!!
!! \~english
!! \brief Interpolates the values at the vertices of a control volume from values at the grid-points surrounding it.
!!
!!
!! @param[in,out] self   The variable of derived type Derivate_Field_IcGrid2D.
!! @param[in,out] CV_val The array containg the values of the physical field at the vertices of control volumes.
!! @param[in]     rcs_GP_val The array containg the values of the physical field at the grid-points.
!! @param[in]     rcID    The index of rectagle region.
!! @param[in]     gp_i    各矩形領域における水平第一成分のインデックス.
!! @param[in]     gp_j    各矩形領域における水平第二成分のインデックス.
!! @param[in]     val_dim The dimension of physical field.
!! @param[in]     CV_num  The number of vertices that a control volume has.
!! @param[in]     idMin 構造型 Derivate_Field_IcGrid2D_Manager の成分 rcs_AGrid の袖領域を含めた配列の最小インデックス.
!! @param[in]     EMin　構造型 Derivate_Field_IcGrid2D_Manager の成分 rcs_AGrid の袖領域を含めない配列の最小インデックス.
!!
subroutine calc_CV5_val( &
	& self, &
	& CV5_val, & 
	& rcs_GP_val, dataNearNPole, dataNearSPole, &
	& rcID, gp_i, gp_j, cvVal_dim, val_dim,     & 
	& iIdMin, jIdMin, rcIDMin, glEMin, glEMax   & 
	& )

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self
  integer, intent(in) :: cvVal_dim, val_dim, iIdMin, jIdMin, rcIDMin, glEMin, glEMax
  real(DP), intent(out) :: CV5_val(cvVal_dim, 5)
  real(DP), intent(in) :: rcs_GP_val(:,iIdMin:,jIdMin:,rcIDMin:)
  real(DP), intent(in) :: dataNearNPole(:,glEMin:,glEMin:,1:)
  real(DP), intent(in) :: dataNearSPole(:,glEMax-1:,glEMax-1:,6:) 
  integer, intent(in) :: rcID, gp_i, gp_j

  ! 作業変数
  ! Work variable
  !
  integer :: cvID, k
  real(DP) :: GP_val3(3,val_dim)
  real(DP) :: area3(3)

  ! 実行文 ; Executable statement
  !

  do cvID=1, 5

    ! 矩形領域 rcID 内の格子点 (gp_i, gp_j) に付随するコントロールボリュームの一頂点 cvID と近傍の格子点が作る 3 つの三角形の
    ! 面積を取得する.
    area3(:) = self%rcs_areaGPP(:, cvID, gp_i,gp_j,rcID)

    !
    if( rcID <= 5 ) then
      if ( gp_i==glEMin .and. gp_j==glEMin )then
       GP_val3(1,:) = dataNearNPole(:, glEMin, glEMin, cvID); GP_val3(2,:) = dataNearNPole(:, glEMin+1, glEMin, cvID);
       GP_val3(3,:) = dataNearNPole(:, glEMin, glEMin+1, cvID);        
      else
        do k=1,3
          GP_val3(k,:) = rcs_GP_val(:, self%rcU_CVGPindex(cvID,k,1, gp_i,gp_j), self%rcU_CVGPindex(cvID,k,2, gp_i,gp_j), rcID)
        end do
        
      end if

    else
    
      if ( gp_i==glEMax .and. gp_j==glEMax )then
        GP_val3(1,:) = dataNearSPole(:, glEMax, glEMax, 11-cvID);   GP_val3(2,:) = dataNearSPole(:, glEMax-1, glEMax, 11-cvID);
        GP_val3(3,:) = dataNearSPole(:, glEMax, glEMax-1, 11-cvID);
      else
        do k=1,3
          GP_val3(k,:) = rcs_GP_val(:, self%rcL_CVGPindex(cvID,k,1, gp_i,gp_j), self%rcL_CVGPindex(cvID,k,2, gp_i,gp_j), rcID)
        end do
      end if

    end if

    ! コントロールボリュームの各頂点(点 G)上の物理場の値を, 近傍の格子点(点 P0,P1,P2)上の物理場の値から補間する.
    ! 各格子点の物理場の値は, コントロールボリュームの頂点と格子点が作る三角形(GP0P1, GP1P2, GP2P0)の面積によって重みを付けられる.
		 CV5_val(1:cvVal_dim, cvID) = &
      &   area3(1) * GP_val3(1,:) &
      & + area3(2) * GP_val3(2,:)  &
      & + area3(3) * GP_val3(3,:) 
      
  end do

end subroutine calc_CV5_val

subroutine calc_CV6_val( & 
	& self, &
	& CV6_val, &
	& CV6GPVal9, CV6Area3, &
	& rcID, gp_i, gp_j, cvVal_dim, val_dim &
	& )

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self
  integer, intent(in) :: cvVal_dim, val_dim, gp_i, gp_j
  real(DP), intent(out) :: CV6_val(cvVal_dim, 6)
  real(DP), intent(in) :: CV6GPVal9(val_dim, gp_i-1:gp_i+1, gp_j-1:gp_j+1)
  real(DP), intent(in) :: CV6Area3(3,6)
  integer, intent(in) :: rcID

  ! 作業変数
  ! Work variable
  !
  integer :: cvID
  integer, pointer :: CV6GPindex1(:,:)
  integer, pointer :: CV6GPindex2(:,:)

  ! 実行文 ; Executable statement
  !

  !CV6GPVal9 = CV6GPVal9_

  if( rcID <= 5 ) then
    CV6GPindex1 => self%rcU_CVGPindex(:,:,1, gp_i,gp_j)
    CV6GPindex2 => self%rcU_CVGPindex(:,:,2, gp_i,gp_j)
	else
	  CV6GPindex1 => self%rcL_CVGPindex(:,:,1, gp_i,gp_j)
	  CV6GPindex2 => self%rcL_CVGPindex(:,:,2, gp_i,gp_j)
	end if

   ! コントロールボリュームの各頂点(点 G)上の物理場の値を, 近傍の格子点(点 P0,P1,P2)上の物理場の値から補間する.
   ! 各格子点の物理場の値は, コントロールボリュームの頂点と格子点が作る三角形(GP0P1, GP1P2, GP2P0)の面積によって重みを付けられる.
  do cvID=1, 6

		CV6_val(1:cvVal_dim, cvID) = &
      &   CV6Area3(1, cvID) * CV6GPVal9(1:val_dim, CV6GPindex1(cvID,1), CV6GPindex2(cvID,1)) &
      & + CV6Area3(2, cvID) * CV6GPVal9(1:val_dim, CV6GPindex1(cvID,2), CV6GPindex2(cvID,2))  &
      & + CV6Area3(3, cvID) * CV6GPVal9(1:val_dim, CV6GPindex1(cvID,3), CV6GPindex2(cvID,3))

  end do

end subroutine calc_CV6_val


subroutine calc_CV5_gradVec( &
	& self, &
	& CV5_val, & 
	& rcs_GP_val, dataNearNPole, dataNearSPole, &
	& rcID, gp_i, gp_j, cvVal_dim, val_dim,     & 
	& iIdMin, jIdMin, rcIDMin, glEMin, glEMax   & 
	& )

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self
  integer, intent(in) :: cvVal_dim, val_dim, iIdMin, jIdMin, rcIDMin, glEMin, glEMax
  real(DP), intent(out) :: CV5_val(cvVal_dim, 5)
  real(DP), intent(in) :: rcs_GP_val(:,iIdMin:,jIdMin:,rcIDMin:)
  real(DP), intent(in) :: dataNearNPole(:,glEMin:,glEMin:,1:)
  real(DP), intent(in) :: dataNearSPole(:,glEMax-1:,glEMax-1:,6:) 
  integer, intent(in) :: rcID, gp_i, gp_j

  ! 作業変数
  ! Work variable
  !
  integer :: cvID, k, kup
  real(DP) :: GP_val3(3,val_dim)
  real(DP) :: area3(3), cv_val, CVPos(3,5), GP0(3), GP_DirVec(3)
  real(DP) :: grad_vec_tmp(3)
  type(FVM_HDerivate_CVInfo), pointer :: CVInfo => null()

  ! 実行文 ; Executable statement
  !

  CVInfo => self%CVInfo(gp_i, gp_j, rcID)
  GP0(:) = self%icgrid%rcs_AGrid(1:3, gp_i,gp_j,rcID)
  CVPos(:,:) = self%icgrid%rcs_CVInfo(gp_i,gp_j,rcID)%cv_pts(1:3, 1:5)
  
  do cvID=1, 5

    ! 矩形領域 rcID 内の格子点 (gp_i, gp_j) に付随するコントロールボリュームの一頂点 cvID と近傍の格子点が作る 3 つの三角形の
    ! 面積を取得する.
    area3(:) = self%rcs_areaGPP(:, cvID, gp_i,gp_j,rcID)

    !
    if( rcID <= 5 ) then
      if ( gp_i==glEMin .and. gp_j==glEMin )then
       GP_val3(1,:) = dataNearNPole(:, glEMin, glEMin, cvID); GP_val3(2,:) = dataNearNPole(:, glEMin+1, glEMin, cvID);
       GP_val3(3,:) = dataNearNPole(:, glEMin, glEMin+1, cvID);        
      else
        do k=1,3
          GP_val3(k,:) = rcs_GP_val(:, self%rcU_CVGPindex(cvID,k,1, gp_i,gp_j), self%rcU_CVGPindex(cvID,k,2, gp_i,gp_j), rcID)
        end do
      end if

    else
    
      if ( gp_i==glEMax .and. gp_j==glEMax )then
        GP_val3(1,:) = dataNearSPole(:, glEMax, glEMax, 11-cvID);   GP_val3(2,:) = dataNearSPole(:, glEMax-1, glEMax, 11-cvID);
        GP_val3(3,:) = dataNearSPole(:, glEMax, glEMax-1, 11-cvID);
      else
        do k=1,3
          GP_val3(k,:) = rcs_GP_val(:, self%rcL_CVGPindex(cvID,k,1, gp_i,gp_j), self%rcL_CVGPindex(cvID,k,2, gp_i,gp_j), rcID)
        end do
      end if

    end if

    ! コントロールボリュームの各頂点(点 G)上の物理場の値を, 近傍の格子点(点 P0,P1,P2)上の物理場の値から補間する.
    ! 各格子点の物理場の値は, コントロールボリュームの頂点と格子点が作る三角形(GP0P1, GP1P2, GP2P0)の面積によって重みを付けられる.
    grad_vec_tmp(:) = 0d0

		cv_val = &
      &   area3(1) * GP_val3(1,1)  &
      & + area3(2) * GP_val3(2,1)  &
      & + area3(3) * GP_val3(3,1) 
   
    do k=1, 3
       kup = 1 + mod(k, 3)
       grad_vec_tmp(1:3) = grad_vec_tmp(1:3) +  & 
         &  CVInfo%geoArcLength_Tri(k, cvID) * CVInfo%normalUniVec2_Tri(1:3, k,cvID) &
         &  * ( 0.5d0  * ( GP_val3(k,1) + GP_val3(kup,1) ) &
         & - cv_val &
         & )!

    end do

		 CV5_val(1:3, cvID) = grad_vec_tmp(:) / (1.0d0*CVInfo%CVTriArea(cvID))

  end do

end subroutine calc_CV5_gradVec

subroutine calc_CV6_gradVec( & 
	& self, &
	& CV6_val, &
	& CV6GPVal9, CV6Area3, &
	& rcID, gp_i, gp_j, cvVal_dim, val_dim &
	& )

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self
  integer, intent(in) :: cvVal_dim, val_dim, gp_i, gp_j
  real(DP), intent(out) :: CV6_val(cvVal_dim, 6)
  real(DP), intent(in) :: CV6GPVal9(val_dim, gp_i-1:gp_i+1, gp_j-1:gp_j+1)
  real(DP), intent(in) :: CV6Area3(3,6)
  integer, intent(in) :: rcID

  ! 作業変数
  ! Work variable
  !
  integer :: cvID, k, kup
  real(DP) :: cv_val
  integer, pointer :: CV6GPindex1(:,:)
  integer, pointer :: CV6GPindex2(:,:)
  
  real(DP) :: grad_vec_tmp(3)
  type(FVM_HDerivate_CVInfo), pointer :: CVInfo

  ! 実行文 ; Executable statement
  !

  CVInfo => self%CVInfo(gp_i, gp_j,rcID)

  if( rcID <= 5 ) then
    CV6GPindex1 => self%rcU_CVGPindex(:,:,1, gp_i,gp_j)
    CV6GPindex2 => self%rcU_CVGPindex(:,:,2, gp_i,gp_j)
	else
	  CV6GPindex1 => self%rcL_CVGPindex(:,:,1, gp_i,gp_j)
	  CV6GPindex2 => self%rcL_CVGPindex(:,:,2, gp_i,gp_j)
	end if

  ! コントロールボリュームの各頂点(点 G)上の物理場の値を, 近傍の格子点(点 P0,P1,P2)上の物理場の値から補間する.
  ! 各格子点の物理場の値は, コントロールボリュームの頂点と格子点が作る三角形(GP0P1, GP1P2, GP2P0)の面積によって重みを付けられる.

   
  do cvID=1, 6

		cv_val = &
      &   CV6Area3(1, cvID) * CV6GPVal9(1, CV6GPindex1(cvID,1), CV6GPindex2(cvID,1)) &
      & + CV6Area3(2, cvID) * CV6GPVal9(1, CV6GPindex1(cvID,2), CV6GPindex2(cvID,2))  &
      & + CV6Area3(3, cvID) * CV6GPVal9(1, CV6GPindex1(cvID,3), CV6GPindex2(cvID,3))

    grad_vec_tmp(:) = 0d0
    
    do k=1, 3
      kup = 1 + mod(k, 3)
      grad_vec_tmp(:) = grad_vec_tmp(:) + &
        & CVInfo%geoArcLength_Tri(k, cvID) * CVInfo%normalUniVec2_Tri(1:3, k,cvID) &
        &  * ( 0.5d0 * ( &
        &         CV6GPVal9(1, CV6GPindex1(cvID,kup), CV6GPindex2(cvID,kup)) &
        &       + CV6GPVal9(1, CV6GPindex1(cvID,k), CV6GPindex2(cvID,k))     & 
         &       ) &
        &      - cv_Val &
         &     )

    end do
    
		CV6_val(:, cvID) = grad_vec_tmp(:) / CVInfo%CVTriArea(cvID)

  end do

end subroutine calc_CV6_gradVec


!
!> \~japanese
!! \brief コントロールボリュームの各頂点を取り囲む近傍の格子点のインデックスを計算する.
!!
!! @param[in,out] self 構造型 Derivate_Field_IcGrid2D の変数.
!!
!! \~english
!! \brief コントロールボリュームの各頂点を取り囲む近傍の格子点のインデックスを計算する.
!!
!! @param[in,out] self The variable of derived type Derivate_Field_IcGrid2D.
!!
subroutine calc_CVGPindex(self)

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self

  ! 作業変数
  ! Work variable
  !
  integer :: i, j
  integer :: glEMin, glEMax
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo

  ! 実行文 ; Executable statement
  !

  glEMin = get_EffSize_Min(self%icgrid)
  glEMax = get_EffSize_Max(self%icgrid)
  lcRCInfo => get_LocalProcRectInfo(self%icgrid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  do j=jEMin,jEMax
    do i=iEMin,iEMax
      if ( ( i==glEMax .and. j==glEMin ) .or. ( i==glEMin .and. j==glEMax )  ) then
        !
        !
        self%rcU_CVGPindex(1:5,:,:, i,j) = &
          & generate_CV5_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=1 )
        self%rcL_CVGPindex(1:5,:,:, i,j) = &
          & generate_CV5_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=6 )
      else

        !
        self%rcU_CVGPindex(1:6,:,:, i,j) = &
          & generate_CV6_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=1 )
        !
        self%rcL_CVGPindex(1:6,:,:, i,j) = &
          & generate_CV6_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=6 )

      end if
    end do
  end do

  if ( iEMax == glEMax .and. jEMax == glEMax ) then
    self%rcU_CVGPindex(1:5,:,:, glEMax,glEMax) = &
      & generate_CV5_GPindex(self%icgrid, GP_i=glEMax, GP_j=glEMax, rcID=1 )
  end if

  if ( iEMin == glEMin .and. jEMin == glEMin ) then
    self%rcL_CVGPindex(1:5,:,:, glEMin,glEMin) = &
      & generate_CV5_GPindex(self%icgrid, GP_i=glEMin, GP_j=glEMin, rcID=6 )
  end if

end subroutine calc_CVGPindex

subroutine set_HDerivate_CVInfo( &
  & self &
  & )

  !
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self

  !
  !
  !
  integer :: rcID, i, j
  integer :: cvID, cvIDup
  integer :: k, kup, kupup
  integer :: glEMin, glEMax
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo

  real(DP) :: ic_radius
  real(DP) :: pt_CV(3,6), GP0(3), Tri_GP3(3,3)
  real(DP) :: normalVec(3)
  integer :: CV_num
  integer, pointer :: CVGPindex_i(:,:)
  integer, pointer :: CVGPindex_j(:,:)
  integer :: TriGP3_i(3), TriGP3_j(3), TriGP3_rcID(3)
  type(FVM_HDerivate_CVInfo), pointer :: cvInfo_rij
  real(DP) :: tmp(3,3)
  !
  !
  !
  glEMin = get_EffSize_Min(self%icgrid)
  glEMax = get_EffSize_Max(self%icgrid)
  ic_radius = get_IcRadius(self%icgrid)
  lcRCInfo => get_LocalProcRectInfo(self%icgrid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
        
        if( rcID <= 5 ) then
          CVGPindex_i => self%rcU_CVGPindex(:,:,1, i,j)
          CVGPindex_j => self%rcU_CVGPindex(:,:,2, i,j)
	      else
	        CVGPindex_i => self%rcL_CVGPindex(:,:,1, i,j)
	        CVGPindex_j => self%rcL_CVGPindex(:,:,2, i,j)
	      end if
        
        if (  ( i==glEMin .and. j==glEMin ) .or. ( i==glEMax .and. j==glEMin ) .or. &
         &    ( i==glEMin .and. j==glEMax ) .or. ( i==glEMax .and. j==glEMax )      & 
         & ) then
         
          CV_num = 5
        else
          CV_num = 6
        end if

        cvInfo_rij => self%cvInfo(i, j, rcID)
        cvInfo_rij%CV_num = CV_num

        GP0(1:3) = self%icgrid%rcs_AGrid(1:3, i,j,rcID)
        pt_CV(1:3, 1:CV_num) = self%icgrid%rcs_CVInfo(i,j,rcID)%cv_pts(1:3, 1:CV_num)

        do cvID=1, CV_num
           cvIDup = 1+mod(cvID, CV_num)

		        ! コントロールボリュームの各辺の長さを求める. 
           cvInfo_rij%geoArcLength_CV(cvID) = geodesic_arc_length(ic_radius, pt_CV(:, cvIDup), pt_CV(:, cvID))
		        ! コントロールボリュームの各辺の法線ベクトルを求める.            
           cvInfo_rij%normalUniVec1_CV(:, cvID) = vec_normalize(cross( pt_CV(:, cvIDup), pt_CV(:, cvID)) )
           cvInfo_rij%normalUniVec2_CV(:, cvID) = vec_normalize(cross( pt_CV(:, cvIDup) - pt_CV(:, cvID), GP0(:)) )
           cvInfo_rij%tangentUniVec_CV(:, cvID) = vec_normalize( pt_CV(:, cvIDup) - pt_CV(:, cvID) )

            !
            !
            
            !
           if (rcID <= 5 .and. i==glEMin .and. j==glEMin ) then
             TriGP3_i(1) = glEMin; TriGP3_i(2) = glEMin+1; TriGP3_i(3) = glEMin; 
             TriGP3_j(1) = glEMin; TriGP3_j(2) = glEMin; TriGP3_j(3) = glEMin+1;
             TriGP3_rcID(:) = cvID                
           else if (rcID > 5 .and. i==glEMax .and. j==glEMax ) then
             TriGP3_i(1) = glEMax; TriGP3_i(2) = glEMax-1; TriGP3_i(3) = glEMax; 
             TriGP3_j(1) = glEMax; TriGP3_j(2) = glEMax; TriGP3_j(3) = glEMax-1; 
             TriGP3_rcID(:) = 11-cvID                
           else
             TriGP3_i(:) = CVGPindex_i(cvID,:); TriGP3_j(:) = CVGPindex_j(cvID,:)
             TriGP3_rcID(:) = rcID
		       end if

             
          do k=1, 3
            Tri_GP3(:,k) = self%icgrid%rcs_AGrid(:, TriGP3_i(k),TriGP3_j(k),TriGP3_rcID(k))
          end do
		       
		       do k=1, 3
		         kup = 1 + mod(k,3); kupup = 1 + mod(kup,3)

  		        ! コントロールボリュームを構成する各三角形の各辺の長さを求める. 
             cvInfo_rij%geoArcLength_Tri(k,cvID) = geodesic_arc_length(ic_radius, Tri_GP3(:,kup), Tri_GP3(:,k))
  		        ! コントロールボリュームを構成する各三角形の各辺の法線ベクトルを求める.
  		       !if ( k ==2 ) then
            !  cvInfo_rij%normalUniVec2_Tri(:, k,cvID) = vec_normalize(cross( Tri_GP3(:,kup)-Tri_GP3(:,k), Tri_Gp3(:,kupup)) )
  		       !else
              cvInfo_rij%normalUniVec2_Tri(:, k,cvID) = vec_normalize(cross( Tri_GP3(:,kup), Tri_GP3(:,k)) )
  	         !end if
		       end do
!write(*,*) '*', cvInfo_rij%normalUniVec2_Tri(:, 2,cvID)
!write(*,*) vec_normalize(cross( Tri_GP3(:,3)-Tri_GP3(:,2), GP0(:)) )

		       cvInfo_rij%CVTriArea(cvID) = spherical_tri_area(ic_radius, Tri_GP3(:,1), Tri_GP3(:,2), Tri_GP3(:,3))

        end do
!if(rcID==3.and.i==33.and.j==33)then
!do cvID=1,5
!write(*,*) cvID,':',cvInfo_rij%CVTriArea(cvID)
!do k=1,3
!tmp(:,k) = cvInfo_rij%normalUniVec2_Tri(:, k,cvID)*cvInfo_rij%geoArcLength_Tri(k,cvID)
!write(*,*) tmp(:,k)
!end do
!write(*,*) sum(tmp, 2)
!end do
!
!end if        
      end do
    end do
  end do

end subroutine set_HDerivate_CVInfo

!
!> \~japanese
!! \brief
!!
!! @param[in,out] self 構造型 Derivate_Field_IcGrid2D の変数.
!!
!! \~english
!! \brief
!!
!! @param[in,out] self The variable of derived type Derivate_Field_IcGrid2D.
!!
subroutine calc_GPP_area_weight( &
  & self &   ! (inout)
  & )

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self

  ! 作業変数
  ! Work variables
  !
  type(IcGrid2D_FVM_CV), pointer :: cv_info
  integer :: rcID, i, j, cvID
  integer :: glEMin, glEMax
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  integer :: glIdMin
  type(LocalProcRectInfo), pointer :: lcRCInfo

  real(DP) :: ic_radius
  integer :: CV5_GPindx(5,3,2), CV6_GPindx(6,3,2)
	real(DP) :: areaGPP(3)
	
  ! 実行文 ; Executable statement
  !

  glEMin = get_EffSize_Min(self%icgrid)
  glEMax = get_EffSize_Max(self%icgrid)
  glIdMin = get_IdMin(self%icgrid)
  ic_radius = get_IcRadius(self%icgrid)
  lcRCInfo => get_LocalProcRectInfo(self%icgrid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
      
        cv_info => self%icgrid%rcs_CVInfo(i,j,rcID)
          
        if ( ( rcID > 5 .and. i==glEMin .and. j==glEMin ) .or. ( i==glEMax .and. j==glEMin ) .or. &
             ( i==glEMin .and. j==glEMax ) .or. ( rcID <= 5 .and. i==glEMax .and. j==glEMax ) ) then
          !
          CV5_GPindx(:,:,:) = generate_CV5_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=rcID )

          do cvID=1,5
            areaGPP(1:3) = calc_usual_area3( CV_GPindx=CV5_GPindx(cvID, :,:), &
              & g=cv_info%cv_pts(:, cvID), &
              & rc_AGrid=self%icgrid%rcs_AGrid(:,:,:,rcID), ic_radius=ic_radius, glIdMin=glIdMin )
              
            self%rcs_areaGPP(1:3, cvID, i,j,rcID) = areaGPP(1:3) / sum(areaGPP(:))
          end do
        else if( rcID <= 5 .and. i==glEMin .and. j==glEMin ) then
          !
          do cvID=1,5
            areaGPP(1:3) = calc_area3_from_4pts( ic_radius=ic_radius, &
              & p1=self%icgrid%rcs_AGrid(:, glEMin,glEMin,cvID), &
              & p2=self%icgrid%rcs_AGrid(:, glEMin+1,glEMin,cvID), &
              & p3=self%icgrid%rcs_AGrid(:, glEMin,glEMin+1,cvID), &
              & g=cv_info%cv_pts(:,cvID) )
            
            self%rcs_areaGPP(1:3, cvID, i,j,rcID) = areaGPP(1:3) / sum(areaGPP(:))
          end do

        else if( rcID > 5 .and. i==glEMax .and. j==glEMax ) then
          !
          do cvID=1,5
            areaGPP(1:3) = calc_area3_from_4pts( ic_radius=ic_radius, &
              & p1=self%icgrid%rcs_AGrid(:,glEMax,glEMax, 11-cvID), &
              & p2=self%icgrid%rcs_AGrid(:,glEMax-1,glEMax, 11-cvID), &
              & p3=self%icgrid%rcs_AGrid(:,glEMax,glEMax-1, 11-cvID), &
              & g=cv_info%cv_pts(:,cvID) )
              
             self%rcs_areaGPP(1:3, cvID, i,j,rcID) = areaGPP(1:3) / sum(areaGPP(:))
          end do

        else
          !
          CV6_GPindx(:,:,:) = generate_CV6_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=rcID )

          do cvID=1,6
            areaGPP(1:3) = calc_usual_area3(     &
              & CV_GPindx=CV6_GPindx(cvID, :,:), &
              & g=cv_info%cv_pts(:,cvID),        &
              & rc_AGrid=self%icgrid%rcs_AGrid(:,:,:,rcID), ic_radius=ic_radius, glIdMin=glIdMin )
            
            self%rcs_areaGPP(1:3, cvID, i,j,rcID) = areaGPP(1:3) / sum(areaGPP(:))

          end do

        end if
      end do
    end do
  end do

end subroutine calc_GPP_area_weight

!
!> \~japanese
!! \brief
!!
!! \param[in] CV_GPindex
!! \param[in] g          重心の位置ベクトル.
!! \param[in] rc_AGrid   対象となる矩形領域の格子座標データを保持する配列.
!! \param[in] ic_rudius  正二十面体を内包する球の半径.
!! \param[in] idMin      配列 rc_AGrid の第 1 (or 2) 次元の添字の下限.
!! \return
!!
!! \~english
!! \brief
!!
!! \param[in] CV_GPindex
!! \param[in] g         Point vector of Gravitational point.
!! \param[in] rc_AGrid
!! \param[in] ic_rudius The radius of a sphere which an icosahedron is embedded into.
!! \param[in] idMin     The lower bound index of array 'rc_AGrid' in the first (or second dimension).
!! \return
!!
function calc_usual_area3( &
  & CV_GPindx, g, rc_AGrid, ic_radius, glIdMin &  ! (in)
  & ) result(area3)

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: glIdMin
  integer, intent(in) :: CV_GPindx(3,2)
  real(DP), intent(in) :: g(3)
  real(DP), intent(in) :: rc_AGrid(:,glIdMin:,glIdMin:)
  real(DP), intent(in) :: ic_radius

  ! 作業変数
  ! Work variable
  !
  real(DP) area3(3)

  ! 実行文 ; Executable statement
  !

  area3(:) = calc_area3_from_4pts( &
    & p1=rc_AGrid(:, CV_GPindx(1,1), CV_GPindx(1,2)), &
    & p2=rc_AGrid(:, CV_GPindx(2,1), CV_GPindx(2,2)), &
    & p3=rc_AGrid(:, CV_GPindx(3,1), CV_GPindx(3,2)), &
    & g=g(:), &
    & ic_radius=ic_radius )

end function calc_usual_area3

!
!> \~japanese
!! \brief
!!
!! \param[in] p1        点 P1 の位置ベクトル.
!! \param[in] p2        点 P2 の位置ベクトル.
!! \param[in] p3        点 P3 の位置ベクトル.
!! \param[in] g         重心の位置ベクトル.
!! \param[in] ic_radius 正二十面体を内包する球の半径.
!! \return
!!
!! \~english
!! \brief
!!
!! \param[in] p1  The position vector of a point P1.
!! \param[in] p2  The position vector of a point P2.
!! \param[in] p3  The position vector of a point P3.
!! \param[in] g         Point vector of Gravitational point.
!! \param[in] ic_radius The radius of a sphere which an icosahedron is embedded into.
!! \return
!!
function calc_area3_from_4pts(p1, p2, p3, g, ic_radius) result(area3)

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: p1(3), p2(3), p3(3)
  real(DP), intent(in) :: g(3)
  real(DP), intent(in) :: ic_radius
  real(DP) area3(3)

  ! 実行文 ; Executable statement
  !

  area3(1) = spherical_tri_area(ic_radius, g, p2, p3)
  area3(2) = spherical_tri_area(ic_radius, g, p3, p1)
  area3(3) = spherical_tri_area(ic_radius, g, p1, p2)

end function calc_area3_from_4pts

end module Derivate_Field_IcGrid2D_Manager
