!> \~japanese
!! \brief 3 次元正二十面格子の鉛直座標に地形に沿う座標系を適用するための, 補助手続きを提供するモジュール. 
!!
!! \par TFCoord_IcGrid3D_helper モジュールについて
!!   - 幾何座標系と地形に沿う座標系間の鉛直座標変換や物理場データの補間を手続きを提供する. 
!!   - 地形に沿う座標系におけるメトリック(\f$G^{1/2}, G^{\Dvect{Z}}\f$)を計算し, 管理する. 
!!   - 地形に沿う座標系において, 物理場の全球平均値を求める手続きを提供する. 
!!
!! \~english
!! \brief 3 次元正二十面格子の鉛直座標に地形に沿う座標系を適用するための, 補助手続きを提供するモジュール. 
!!
!! \par About TFCoord_IcGrid3D_helper module
!!   - 幾何座標系と地形に沿う座標系間の鉛直座標変換や物理場データの補間を手続きを提供する. 
!!   - 地形に沿う座標系におけるメトリック(\f$G^{1/2}, G^{\Dvect{Z}}\f$)を計算し, 管理する. 
!!   - 地形に沿う座標系において, 物理場の全球平均値を求める手続きを提供する. 
!!
!! \~
!! \since IGMBaseLib 2.0
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module TFCoord_IcGrid3D_helper

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

  ! Message Passing Interface Library
  !
  !
#ifdef LIB_MPI
  use mpi
#endif  

  !* gtool ****
  !*

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: &
    & DP

  !* IGMBaseLib ****
  !*

  use igmcore_Interpolation, only: &
    & CubicSpline_Interpolation, CubicSpline_eval, &
    & CubicSpline_SData
    
  !
  !
  !
  use TerrainFCoordinate_helper, only: &
    & TF_to_Z, Z_to_TF

  ! 正二十面格子データの管理
  ! Manager of icosaheral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & check_pole, NOT_POLE_FLAG,  &
    & 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
    
  ! 3 次元正二十面格子データの管理
  ! Manager of 3D icosahedral grid data
  !
  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_icgrid2D, get_vertical_level_num, &
    & set_IcGrid3D_VerticalCoord_name, &
    & set_IcGrid3D_VerticalCoord_long_name, &
    & set_IcGrid3D_VerticalCoord_units

  ! 1 次元物理場データの管理
  ! Manager of 1D physical field data.
  !
  use Field_Grid1D_Manager, only: &
    & Field_Grid1D

  ! 正二十面格子上の物理場の管理
  ! Manager of the physical field on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_icGrid2D_Init, Field_IcGrid2D_Final, &
    & paste_field2D_margin

  ! 3 次元正二十面格子上の物理場の管理
  ! Manager of the physical field on 3D icosahedral grid
  !
  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & Field_IcGrid3D_Init, Field_IcGrid3D_Final, &
    & get_icgrid3D, get_field_layer, get_field3D_rank, &
    & paste_field3D_margin


  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  !>\~japanese
  !!\brief IcGrid3D_FVM_Manager モジュールに対して地形に沿う座標系を提供するために必要な情報を管理する構造型. 
  !!
  !!\~english
  !!\brief IcGrid3D_FVM_Manager モジュールに対して地形に沿う座標系を提供するために必要な情報を管理する構造型. 
  !!
  type, public :: TFCoord_IcGrid3D

    !> \~japanese 下部境界の地形高度データを保持する構造型 Field_IcGrid2D の変数を指すポインタ. 
    !! \~english A pointer to a variable of derived type Field_IcGrid2D which contains height data at lower bottom. 
    type(Field_IcGrid2D), pointer :: Z_s
    
    !> \~japanese
    !! \~english
    real(DP) :: Z_T

    !> \~japanese
    !! \~english    
    type(IcGrid3D_FVM), pointer :: icgrid3D

    !> \~japanese
    !! \~english
    type(Field_icGrid2D) :: G_1div2

    !> \~japanese
    !! \~english
    type(Field_icGrid3D) :: G_Z

  end type TFCoord_IcGrid3D

  ! 公開手続き
  ! Public procedures
  !
  public :: TFCoord_IcGrid3D_Init, TFCoord_IcGrid3D_Final

  public :: TF_to_Z_IcGrid3D, Z_to_TF_IcGrid3D, &
    &       calc_metric_G_1div2, calc_metric_G_Z, calc_gamma, &
    &       interpolate_FieldZ1D_into_FieldTF3D,              &
    &       global_mean_Field_TFIcGrid3Dr, global_mean_Field_TFIcGrid3Dz

  public :: get_TFMetric_G_1div2

contains

!
!> \~japanese
!! \brief 構造型 TFCoord_IcGrid3D の変数を初期化する. 
!!
!! @param[in,out] self      構造型 TFCoord_IcGrid3D の変数. 
!! @param[in,out] icgrid3D  地形に沿う座標系を関連付ける 3 次元正二十面体格子の管理オブジェクト. 
!! @param[in]     Z_s       地形の高度データを保持する水平正二十面格子上の物理場データの管理オブジェクト. 
!! @param[in]     Z_T       モデル上端の高度. 
!!
!! \~english
!! \brief Initialize a variable of derived type TFCoord_IcGrid3D.  
!!
!! @param[in,out] self      構造型 TFCoord_IcGrid3D の変数. 
!! @param[in,out] icgrid3D  地形に沿う座標系を関連付ける 3 次元正二十面体格子の管理オブジェクト. 
!! @param[in]     Z_s       地形の高度データを保持する水平正二十面格子上の物理場データの管理オブジェクト. 
!! @param[in]     Z_T       モデル上端の高度. 
!!
subroutine TFCoord_IcGrid3D_Init( &
	& self, icgrid3D, &                ! (inout)
	& Z_s, Z_T        &                ! (in)
	& )

  ! 宣言文 ; Declaration statements
  !
  type(TFCoord_IcGrid3D), intent(inout) :: self
  type(IcGrid3D_FVM), intent(inout), target :: icgrid3D
  type(Field_IcGrid2D), intent(in), target :: Z_s
  real(DP), intent(in) :: Z_T

  ! 実行文 ; Executable statements
  !
  self%icgrid3D => icgrid3D
  self%Z_s => Z_S
  self%Z_T = Z_T

  call Field_IcGrid2D_Init(self=self%G_1div2, icgrid=get_icgrid2D(icgrid3D), name='G_1div2', rank=1)
  call Field_IcGrid3D_Init(self=self%G_Z, icgrid3D=icgrid3D, name='G_Z', rank=3)

  ! IcGrid3D_FVM の成分として関係づけられている鉛直座標のメタデータを, 地形に沿う座標系用に設定する.
  !

  call set_IcGrid3D_VerticalCoord_name(icgrid3D, 'xi')
  call set_IcGrid3D_VerticalCoord_long_name(icgrid3D, 'Terrain following coordinate')
  call set_IcGrid3D_VerticalCoord_units(icgrid3D, '1')

end subroutine TFCoord_IcGrid3D_Init

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

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

  ! 実行文 ; Executable statements
  !

  call Field_IcGrid2D_Final(self%G_1div2)
  call Field_IcGrid3D_Final(self%G_Z)

end subroutine TFCoord_IcGrid3D_Final

!
!> \~japanese
!! \brief
!!
!! \~english
!! \brief 
!!
function get_TFMetric_G_1div2( &
  & self &  ! (inout)
  & ) result(xy_G1div2)

  ! 宣言文 ; Declaration statements
  !
  type(TFCoord_IcGrid3D), intent(in), target :: self
  type(Field_IcGrid2D), pointer :: xy_G1div2

  ! 実行文 ; Executable statements
  !
  xy_G1div2 => self%G_1div2

end function get_TFMetric_G_1div2

!
!> \~japanese
!! \brief
!!
!! \~english
!! \brief 
!!
function Z_to_TF_IcGrid3D( &
  & self, z, rcID, i, j    & ! (in)
  & ) result(xi)

  ! 宣言文 ; Declaration statements
  !
  type(TFCoord_IcGrid3D), intent(in) :: self
  real(DP), intent(in) :: z
  integer, intent(in) :: rcID
  integer, intent(in) :: i
  integer, intent(in) :: j
  real(DP) xi

  ! 実行文 ; Executable statements
  !
  xi = Z_to_TF(z, self%Z_s%val(1,i,j,rcID), self%Z_T)

end function Z_to_TF_IcGrid3D
!
!> \~japanese
!! \brief
!!
!! \~english
!! \brief 
!!
function TF_to_Z_IcGrid3D( &
  & self, xi, rcID, i, j   &  ! (in)
  & ) result(z)

  ! 宣言文 ; Declaration statements
  !
  type(TFCoord_IcGrid3D), intent(in) :: self
  real(DP), intent(in) :: xi
  integer, intent(in) :: rcID
  integer, intent(in) :: i
  integer, intent(in) :: j
  real(DP) z

  ! 実行文 ; Executable statements
  !
  z = TF_to_Z(xi, self%Z_s%val(1,i,j,rcID), self%Z_T)

end function TF_to_Z_IcGrid3D

!
!> \~japanese
!! \brief
!!
!! \~english
!! \brief 
!!
function calc_gamma( &
  & z, a             & ! (in)
  & ) result(val)
  
  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: z
  real(DP), intent(in) :: a
  real(DP) val

  ! 実行文 ; Executable statements
  !
  val = 1.0d0 + z / a

end function calc_gamma

!
!> \~japanese
!! \brief  地形に沿う座標系におけるメトリックの一つである\f$ G^{1/2} \f$を計算する.
!!
!! @param[in,out] self
!!
!! \~english
!! \brief
!!
!! @param[in,out] self
!!
subroutine calc_metric_G_1div2( & 
  & self                        &  ! (inout)
  & )

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

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

  ! 実行文 ; Declaration statements
  !

  lcRCInfo => get_LocalProcRectInfo(get_icgrid2D(self%icgrid3D))

  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
        self%G_1div2%val(1, i,j,rcID) = ( self%Z_T - self%Z_s%val(1, i,j,rcID) ) / self%Z_T
      end do
    end do
  end do

  !
  call paste_field2D_margin(self%G_1div2)

end subroutine calc_metric_G_1div2

!
!> \~japanese
!! \brief 地形に沿う座標系におけるメトリックの一つである\f$ G^\Dvect{Z} \f$を計算する.
!!
!! @param[in,out] self
!! @param[in,out] derivate
!!
!! \~english
!! \brief
!!
!! @param[in,out] self
!! @param[in,out] derivate
!!
subroutine calc_metric_G_Z( & 
  & self, derivate          &  ! (inout)
  & )

  ! モジュール引用 ; Use statements
  !
  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & gradient_op

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

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid3D) :: xyT_Height
  type(Field_IcGrid2D) :: HDiv_Zk
  integer :: rcID, i,j,verticalLvID
  integer :: iEMin, iEMax, jEMin, jEMax, KMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  real(DP) :: xi

  ! 実行文 ; Executable statements
  !

  ! 水平方向, 鉛直方向のデータを保持する配列境界インデックスを取得する.
  !

  KMax = get_vertical_level_num(self%icgrid3D)
  lcRCInfo => get_LocalProcRectInfo(get_icgrid2D(self%icgrid3D))

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

  ! 作業変数を初期化し, メモリを確保する.
  !

  call Field_IcGrid3D_Init(xyT_Height, self%icgrid3D, 'xyT_Height', rank=1)
  call Field_IcGrid2D_Init(HDiv_Zk, get_icgrid2D(self%icgrid3D), 'HDiv_Zk', rank=3)

  ! 各鉛直レベルに対して, G_Z を求める.
  ! 始めに, 地形に沿う座標系における幾何高度場 z(lon,lat,xi) の水平勾配場を求めるために,
  ! xi 系の各鉛直レベルにおいて, z の場を計算する.
  !
  do verticalLvID=1, KMax

    ! verticalLvID に対応する鉛直高度 xi を取得する.
    xi = self%icgrid3D%verticalGrid(verticalLvID)

    ! 等 xi 面における幾何高度 z の水平分布を求める.
    ! (下端地形が変化する場合には, 等 xi 面と等 z 面は一致しない)

    do rcID=rcIDMin, rcIDMax
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          xyT_Height%field2Ds(verticalLvID)%val(1, i,j,rcID) = TF_to_Z_IcGrid3D(self, xi, rcID, i, j)

        end do
      end do
    end do
  end do

  ! 袖領域の物理場データを貼り合わせる.
  call paste_field3D_margin(xyT_Height)

  ! 上で求めた, 幾何高度場 z(lon,lat,xi) に対して, 水平勾配演算子を作用させる.
  ! そして, その結果を G_1/2 で割ることで, 最終的な G_Z の場を求める.
  !

  do verticalLvID=1, KMax-1

    call gradient_op(derivate, get_field_layer(xyT_Height, verticalLvID), HDiv_Zk)

    do rcID=rcIDMin, rcIDMax
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          self%G_Z%field2Ds(verticalLvID)%val(1:3, i,j,rcID) = &
            & - HDiv_Zk%val(1:3, i,j,rcID) / self%G_1div2%val(1, i,j,rcID)

        end do
      end do
    end do
  end do

  ! 領域上端では, G_Z は恒等的にゼロとなる.
  self%G_Z%field2Ds(KMax)%val = 0.0d0

  ! 袖領域の物理場データを貼り合わせる.
  call paste_field3D_margin(self%G_z)

  ! 作業変数のメモリを解放する.
  !

  call Field_IcGrid3D_Final(xyT_Height)
  call Field_IcGrid2D_Final(HDiv_Zk)

end subroutine calc_metric_G_Z

subroutine interpolate_FieldZ1D_into_FieldTF3D( & 
  & self, Z_field1D,                            &  ! (in)
  & TF_field3D                                  &  ! (inout) 
  & )

  ! 宣言文 ; Declaration statements
  !
  type(TFCoord_IcGrid3D), intent(in) :: self
  type(Field_Grid1D), intent(in) :: Z_field1D
  type(Field_IcGrid3D), intent(inout) :: TF_field3D

  ! 作業変数
  ! Work variables
  !
  integer :: Z_field_size, fieldRank
  
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid2D), pointer :: TF_field3D_k
  integer :: i, j, rcID, kIntr, rankID, S_id
  
  integer :: EMin, EMax, KMax
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  
  real(DP) :: xi, z
  real(DP) :: gridZ(size(Z_field1D%grid1D(:)))
  type(CubicSpline_SData) :: csSData(size(gridZ)-1)
  
  
  ! 実行文 ; Executable statements
  !
  

  icgrid2D => get_icgrid2D(self%icgrid3D)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)

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

  KMax = get_vertical_level_num(self%icgrid3D)
  Z_field_size = size(Z_field1D%grid1D)

  gridZ(:) = Z_field1D%grid1D(:)
  fieldRank = size(Z_field1D%val, 2)
  
  do rankID=1, fieldRank

    ! 3 次スプライン補間を用いて, 対象となる一次元物理場データから補間多項式 csSData を生成する. 
    ! 
    csSData(:) = CubicSpline_Interpolation(gridZ(:), Z_field1D%val(:,rankID))
  
    do kIntr=1, KMax
      TF_field3D_k => get_field_layer(TF_field3D, kIntr)
      xi = self%icgrid3D%verticalGrid(kIntr)

      do rcID=rcIDMin, rcIDMax
        !$omp parallel do private(i, S_id, z)
        do j=jEMin, jEMax
          do i=iEMin, iEMax
             
             ! 地形に沿う座標系における一格子点 (rcID,i,j,k) の鉛直高度 xi を, 幾何座標系の高度に変換する. 
            z = TF_to_Z_IcGrid3D(self, xi, rcID, i, j)
              
             ! 格子位置 z　が, 補間多項式の配列 csSData の何番目の要素に対応するかを求める.  
            do S_id=1, Z_field_Size-2
              if (gridZ(S_id) <= z .and. z < gridZ(S_id+1) ) then
                 exit
              end if
            end do
            
             ! 補間多項式から, xi における物理場データの補間値を求める. 
            TF_field3D_k%val(rankID, i,j,rcID) = CubicSpline_eval(csSData(S_id), z)
            
          end do
        end do
      end do
    end do
    
  end do
  
end subroutine interpolate_FieldZ1D_into_FieldTF3D

function global_mean_Field_TFIcGrid3Dr(                             &
  & self, fieldrTF3D, field_rank,         &  ! (in)
  & verticalVirtualRegionWidth            &
  & ) result(gl_mean)
  
  ! 宣言文 ; Declaration statements
  !
  type(TFCoord_IcGrid3D), intent(in) :: self
  type(Field_IcGrid3D), intent(in)  :: fieldrTF3D
  integer, intent(in) :: field_rank
  real(DP) gl_mean(field_rank)
  integer, optional :: verticalVirtualRegionWidth
  
  ! 作業変数 
  ! Work variables
  !
  type(Field_IcGrid2D), pointer :: fieldrTF3D_kU, fieldrTF3D_kL
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dr
  
  integer :: rcID, i, j, kr
  integer :: glEMin, glEMax, krMax, KrMin
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  logical :: isPole
  
  real(DP) :: xiL, xiU
  real(DP) :: Zr, ZLr, ZUr
  real(DP) :: radius_h0, dArea_kr, dvol
  real(DP) :: regionVol, global_regionVol
  real(DP) :: tmp_integ(field_rank), global_tmp_integ(field_rank)

  integer :: ierr
  
  ! 実行文 ; Executable statements
  !
  
  icgrid3Dr => self%icgrid3D
  icgrid2D => get_icgrid2D(icgrid3Dr)

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

  glEMin = get_EffSize_Min(icgrid2D)
  glEMax = get_EffSize_Max(icgrid2D)
  KrMin = 1
  krMax = get_vertical_level_num(icgrid3Dr)
  
  if (present(verticalVirtualRegionWidth)) then
    KrMax = KrMax - verticalVirtualRegionWidth
    KrMin = KrMin + verticalVirtualRegionWidth
  end if

  radius_h0 = get_IcRadius(icgrid2D)

  tmp_integ = 0.0d0
  regionVol = 0.0d0

  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
      
       if ( check_pole(icgrid2D, rcID, i, j) /= NOT_POLE_FLAG ) isPole = .true.
	     
	      if (     ( j < glEMax .and. i > glEMin )                      &
	      	 & .or. ( ( rcID == 1 .or. rcID == 10 ) .and. ( check_pole(icgrid2D, rcID, i, j) /= NOT_POLE_FLAG ) )  &
	      	 & ) then


						do kr=KrMin, krMax-1
         	 fieldrTF3D_kL => get_field_layer(fieldrTF3D, kr)
         	 fieldrTF3D_kU => get_field_layer(fieldrTF3D, kr+1)

         	 xiL = icgrid3Dr%verticalGrid(kr)
         	 xiU = icgrid3Dr%verticalGrid(kr+1)

         	 ZLr = TF_to_Z_IcGrid3D(self, xiL, rcID, i, j)
         	 ZUr = TF_to_Z_IcGrid3D(self, xiU, rcID, i, j)
         	 Zr = 0.5d0 * ( ZLr + ZUr )

         	 dArea_kr = ( 1.0d0 + Zr / radius_h0 )**2 * icgrid2D%rcs_CVInfo(i,j, rcID)%cv_area
         	 dvol = dArea_kr * ( ZUr - ZLr )
         	 
         	 tmp_integ(:) =   tmp_integ(:) &
         	    &            + 0.5d0 * ( fieldrTF3D_kU%val(:, i,j,rcID) + fieldrTF3D_kL%val(:, i,j,rcID) ) * dvol

         	 regionVol = regionVol + dvol

          end do
     
       end if 
     
     end do
  	end do
	end do

  if ( lcRCInfo%mpi_flag ) then
  
    ! MPI による並列化が有効である場合には, 各プロセスから tmp_integ を収集して総和する. 
    !
    
#ifdef LIB_MPI      
    call MPI_AllReduce( tmp_integ, global_tmp_integ, field_rank, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )      

    call MPI_AllReduce( regionVol, global_regionVol, 1, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )      
#endif

  else
    global_tmp_integ(:) = tmp_integ(:)
    global_regionVol = regionVol
  end if

  gl_mean(:) = global_tmp_integ(:) / global_regionVol
 ! write(*,*)'TMP--',  tmp_integ, regionVol, abs(regionVol-4.0d0*3.14d0*radius_h0**2*1.0d4)/(4.0d0*3.14d0*radius_h0**2*1.0d4),  & 
 ! & gl_mean

end function global_mean_Field_TFIcGrid3Dr

function global_mean_Field_TFIcGrid3Dz(    &
  & self, fieldzTF3D, field_rank, tfcoordr, &  ! (in)
  & verticalVirtualRegionWidth            &
  & ) result(gl_mean)

  ! 宣言文 ; Declaration statements
  !
  type(TFCoord_IcGrid3D), intent(in) :: self
  type(Field_IcGrid3D), intent(in)  :: fieldzTF3D
  integer, intent(in) :: field_rank
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordr
  integer, optional :: verticalVirtualRegionWidth
  real(DP) gl_mean(field_rank)

  ! 作業変数 
  ! Work variables
  !
  type(Field_IcGrid2D), pointer :: fieldzTF3D_k
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dr, icgrid3Dz
  integer :: rcID, i, j, k
  integer :: glEMin, glEMax, KzMin, kzMax
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  logical :: isPole
  
  real(DP) :: xiL, xiU
  real(DP) :: Zz, ZL, ZU
  real(DP) :: radius_h0, dArea_kr, dvol
  real(DP) :: regionVol, global_regionVol
  real(DP) :: tmp_integ(field_rank), global_tmp_integ(field_rank)

  integer :: ierr
  
  ! 実行文 ; Executable statements
  !

  icgrid3Dz => self%icgrid3D
  icgrid3Dr => tfcoordr%icgrid3D
  icgrid2D => get_icgrid2D(icgrid3Dz)
  glEMin = get_EffSize_Min(icgrid2D)
  glEMax = get_EffSize_Max(icgrid2D)
 
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  kzMax = get_vertical_level_num(icgrid3Dz)
  KzMin = 1
  
  if (present(verticalVirtualRegionWidth)) then
    KzMax = KzMax - verticalVirtualRegionWidth
    KzMin = KzMin + verticalVirtualRegionWidth
  end if


  radius_h0 = get_IcRadius(icgrid2D)

  tmp_integ = 0.0d0
  regionVol = 0.0d0

  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
      
	      if (     ( j < glEMax .and. i > glEMin )                      &
	      	 & .or. ( ( rcID == 1 .or. rcID == 10 ) .and. ( check_pole(icgrid2D, rcID, i, j) /= NOT_POLE_FLAG ) )  &
	      	 & ) then
	      	
						do k=KzMin, kzMax
         	 fieldzTF3D_k => get_field_layer(fieldzTF3D, k)

         	 xiL = icgrid3Dr%verticalGrid(k)
         	 xiU = icgrid3Dr%verticalGrid(k+1)

         	 Zz = Z_to_TF_IcGrid3D(self, icgrid3Dz%verticalGrid(k), rcID, i, j)
						 ZU	= Z_to_TF_IcGrid3D(tfcoordr, xiU, rcID, i, j)						
						 ZL	= Z_to_TF_IcGrid3D(tfcoordr, xiL, rcID, i, j)						
							
         	 dArea_kr = ( 1.0d0 + Zz / radius_h0 )**2 * icgrid2D%rcs_CVInfo(i,j, rcID)%cv_area
         	 dvol = dArea_kr * ( ZU - ZL )

         	 tmp_integ(:) =   tmp_integ(:) + fieldzTF3D_k%val(:, i,j,rcID) * dvol
         	 regionVol = regionVol + dvol

          end do
     
       end if 
     
     end do
  	end do
	end do

  if ( lcRCInfo%mpi_flag ) then
  
    ! MPI による並列化が有効である場合には, 各プロセスから tmp_integ を収集して総和する. 
    !
    
#ifdef LIB_MPI      
    call MPI_AllReduce( tmp_integ, global_tmp_integ, field_rank, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )      

    call MPI_AllReduce( regionVol, global_regionVol, 1, MPI_REAL8, &
      &                 MPI_SUM, MPI_COMM_WORLD, ierr )      
#endif

  else
    global_tmp_integ(:) = tmp_integ(:)
    global_regionVol = regionVol
  end if

  gl_mean(:) = global_tmp_integ(:) / global_regionVol
  !write(*,*)'TMP--',  tmp_integ, regionVol, gl_mean

end function global_mean_Field_TFIcGrid3Dz

end module TFCoord_IcGrid3D_helper
