3 次元 (xyz 方向) 等間隔交互格子 格子点設定モジュール

  function AvrXYZ_pyz(pyz_Var)
    ! pyz 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: AvrXYZ_pyz                             ! 出力
    
    ! 初期化
    AvrXYZ_pyz = 0.0d0
    
    AvrXYZ_pyz = IntXYZ_pyz(pyz_Var)/ (sum(x_dx(1:nx))*sum(y_dy(1:ny))*sum(z_dz(1:nz)))
    
  end function AvrXYZ_pyz

  function AvrXYZ_xqz(xqz_Var)
    ! xqz 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: AvrXYZ_xqz                             ! 出力
    
    ! 初期化
    AvrXYZ_xqz = 0.0d0
    
    AvrXYZ_xqz = IntXYZ_xqz(xqz_Var)/ (sum(x_dx(1:nx))*sum(y_dy(1:ny))*sum(z_dz(1:nz)))
    
  end function AvrXYZ_xqz

  function AvrXYZ_xyr(xyr_Var)
    ! xyr 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: xyr_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: AvrXYZ_xyr                             ! 出力
    
    ! 初期化
    AvrXYZ_xyr = 0.0d0

    AvrXYZ_xyr = IntXYZ_xyr(xyr_Var)/ (sum(x_dx(1:nx))*sum(y_dy(1:ny))*sum(z_dz(1:nz)))
    
  end function AvrXYZ_xyr

  function AvrXYZ_xyz(xyz_Var)
    ! xyz 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: AvrXYZ_xyz                             ! 出力
    
    ! 初期化
    AvrXYZ_xyz = 0.0d0
    
    AvrXYZ_xyz = IntXYZ_xyz(xyz_Var)/ (sum(x_dx(1:nx))*sum(y_dy(1:ny))*sum(z_dz(1:nz)))
    
  end function AvrXYZ_xyz

  function IntXYZ_pyz(pyz_Var)
    ! pyz 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: IntXYZ_pyz                             ! 出力
    
    ! 初期化
    IntXYZ_pyz = 0.0d0

    IntXYZ_pyz = IntZ_z(z_IntXY_pyz(pyz_Var))
    
  end function IntXYZ_pyz

  function IntXYZ_xqz(xqz_Var)
    ! xqz 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: IntXYZ_xqz                             ! 出力
    
    ! 初期化
    IntXYZ_xqz = 0.0d0
    
    IntXYZ_xqz = IntZ_z(z_IntXY_xqz(xqz_Var))
    
  end function IntXYZ_xqz

  function IntXYZ_xyr(xyr_Var)
    ! xyr 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: xyr_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: IntXYZ_xyr                             ! 出力
    
    ! 初期化
    IntXYZ_xyr = 0.0d0
    
    IntXYZ_xyr = IntZ_r(a_IntXY_xya(xyr_Var))

  end function IntXYZ_xyr

  function IntXYZ_xyz(xyz_Var)
    ! xyz 格子上の配列に対し領域積分を行う
    
    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: IntXYZ_xyz                             ! 出力
    
    ! 初期化
    IntXYZ_xyz = 0.0d0
    
    IntXYZ_xyz = IntZ_z(a_IntXY_xya(xyz_Var))
    
  end function IntXYZ_xyz

  subroutine axesset_init
    ! 格子点座標配列と格子点間隔配列の初期化    

    ! 暗黙の型宣言禁止
    implicit none

    ! 変数定義
    real(DP),allocatable :: xy_X(:,:)! x 座標(半整数格子, 作業配列)
    real(DP),allocatable :: xy_Y(:,:)! y 座標(半整数格子, 作業配列)
    real(DP),allocatable :: yz_Z(:,:)! z 座標(半整数格子, 作業配列)
    integer              :: unit     ! 設定ファイル用装置番号
    

    !設定ファイルから読み込む出力ファイル情報
    NAMELIST /axesset_nml/ xmin, xmax, ymin, ymax, zmin, zmax

    !設定ファイルから出力ファイルに記載する情報を読み込む
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=axesset_nml)
    close(unit)
    
    ! 配列の上下限の値, 座標値と格子点間隔を設定
    ! * 1 次元用のサブルーチンを用いる
    !
    call x_axis_init
    call y_axis_init
    call z_axis_init
    
    ! MPI 対応
    ! * x 方向のみ
    !
    XMin = XMin           ! XMin = 0 を仮定
    XMax = XMax * nprocs  ! CPU の数だけ領域を拡張する

    ! 3 次元格子点座標配列の設定
    ! * 組み込み関数 spread を用いる. 
    ! * 中間配列として 2 次元格子点座標配列を作り, それを 3 次元に拡張する.
    ! 
    allocate(xy_X(imin:imax,jmin:jmax))
    allocate(xy_Y(imin:imax,jmin:jmax))
    allocate(yz_Z(jmin:jmax,kmin:kmax))
    
    allocate(xyz_X(imin:imax,jmin:jmax,kmin:kmax))
    allocate(xyz_Y(imin:imax,jmin:jmax,kmin:kmax))
    allocate(xyz_Z(imin:imax,jmin:jmax,kmin:kmax))
    allocate(xyz_dX(imin:imax,jmin:jmax,kmin:kmax))
    allocate(xyz_dY(imin:imax,jmin:jmax,kmin:kmax))
    allocate(xyz_dZ(imin:imax,jmin:jmax,kmin:kmax))
    
    xy_X  = spread(x_X, 2,size(y_Y))
    xyz_X = spread(xy_X,3,size(z_Z))
    
    xy_X   = spread(x_dX, 2,size(y_dY))
    xyz_dX = spread(xy_X,3,size(z_dZ))
    
    xy_Y  = spread(y_Y, 1,size(x_X))
    xyz_Y = spread(xy_Y,3,size(z_Z))
    
    xy_Y   = spread(y_dY, 1,size(x_dX))
    xyz_dY = spread(xy_Y,3,size(z_dZ))
    
    yz_Z  = spread(z_Z, 1,size(y_Y))
    xyz_Z = spread(yz_Z,1,size(x_X))
    
    yz_Z   = spread(z_dZ, 1,size(y_dY))
    xyz_dZ = spread(yz_Z,1,size(x_dX))
    
    deallocate(xy_X)
    deallocate(xy_Y)
    deallocate(yz_Z)

    !"myrank == 0" に該当する計算ノードが, 読み込んだ情報を出力
    if (myrank == 0) then 
      call MessageNotify( "M", "axesset_init", "XMin = %f", d=(/XMin/)    )
      call MessageNotify( "M", "axesset_init", "XMax = %f", d=(/XMax/)    )
      call MessageNotify( "M", "axesset_init", "YMin = %f", d=(/YMin/)    )
      call MessageNotify( "M", "axesset_init", "YMax = %f", d=(/YMax/)    )
      call MessageNotify( "M", "axesset_init", "ZMin = %f", d=(/ZMin/)    )
      call MessageNotify( "M", "axesset_init", "ZMax = %f", d=(/ZMax/)    )
    end if
  
  end subroutine axesset_init

  function aqa_avr_aya(aya_Var)
    ! 平均操作を行い y 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: aqa_avr_aya(imin:imax,jmin:jmax,kmin:kmax) 
    integer                 :: jy
    
    ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
    if (jmin == jmax) then 
      aqa_avr_aya = aya_Var
      return
    end if
    
    ! 平均操作
    !
    do jy = jmin, jmax-1
!      aqa_avr_aya(:,jy,:) = &
!        &  (y_dy(jy)*aya_Var(:,jy+1,:) + y_dy(jy+1)*aya_Var(:,jy,:)) &
!        &  * 0.5d0/q_dy(jy)
      aqa_avr_aya(:,jy,:) = (aya_Var(:,jy+1,:) + aya_Var(:,jy,:)) * 0.5d0
    end do
    
    ! jmax 格子上の値
    aqa_avr_aya(:,jmax,:) = 1.0d10
    
  end function aqa_avr_aya

  function paa_avr_xaa(xaa_Var)
    ! 半整数格子点の配列値を整数点上へ返す
    
    real(DP),intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: paa_avr_xaa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: ix

    ! 平均操作: x 座標から p 座標へ返す. imax の値は陽に代入する.
    !     
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !
    do ix = imin, imax-1
!      paa_avr_xaa(ix,:,:) = &
!        &  (x_dx(ix)*xaa_Var(ix+1,:,:) + x_dx(ix+1)*xaa_Var(ix,:,:)) &
!        &  *0.5d0/p_dx(ix)
      paa_avr_xaa(ix,:,:) = (xaa_Var(ix+1,:,:) + xaa_Var(ix,:,:)) * 0.5d0
    end do
    
    ! imax 格子上の値
    paa_avr_xaa(imax,:,:) = 1.0d10
    
  end function paa_avr_xaa

  function pqz_avr_xyz(xyz_Var)
    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax)      ! 入力
    real(DP)             :: pqz_avr_xyz(imin:imax,jmin:jmax,kmin:kmax)  ! 出力
    integer              :: ix, jy                                      ! ループ添字

    ! 2 次元計算(y 方向の配列要素数が 1)の場合, pz_avr_xz と同じ計算を行う.
    if (jmin == jmax) then 
      do ix = imin, imax-1
        pqz_avr_xyz(ix,:,:) = ( xyz_Var(ix,:,:) + xyz_Var(ix+1,:,:) ) * 0.5d0 
      end do

      ! imax 格子上の値
      pqz_avr_xyz(imax,:,:) = 1.0d10

    ! 3 次元計算の場合
    else
      do jy = jmin, jmax-1
        do ix = imin, imax-1
          pqz_avr_xyz(ix,jy,:) = ( xyz_Var(ix,jy,:)   + xyz_Var(ix+1,jy,:) + xyz_Var(ix,jy+1,:) + xyz_Var(ix+1,jy+1,:) ) * 0.25d0 
        end do
      end do

      ! imax 格子上の値
      pqz_avr_xyz(imax,:,:) = 1.0d10
      ! jmax 格子上の値
      pqz_avr_xyz(:,jmax,:) = 1.0d10
    end if
        
  end function pqz_avr_xyz

  function aa_AvrZ_aaz(aaz_Var)
    ! aaz 格子上の配列に対し z 方向平均を行う
    
    real(DP), intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrZ_aaz(imin:imax,jmin:jmax)       ! 出力
    
    ! 初期化
    aa_AvrZ_aaz = 0.0d0
    
    ! 平均
    aa_AvrZ_aaz = aa_IntZ_aaz(aaz_Var)/sum(z_dz(1:nz))
    
  end function aa_AvrZ_aaz

  function aa_IntZ_aaz(aaz_Var)
    ! aaz 格子上の配列に対し z 方向積分を行う
    
    real(DP), intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntZ_aaz(imin:imax,jmin:jmax)       ! 出力
    integer                  :: ix, jy                           ! ループ添字
    
    ! 初期化
    aa_IntZ_aaz = 0.0d0

    ! 積分
    do jy = jmin, jmax
      do ix = imin, imax
        aa_IntZ_aaz(ix,jy) = IntZ_z(aaz_Var(ix,jy,:))
      end do
    end do
    
  end function aa_IntZ_aaz

  function aar_avr_aaz(aaz_Var)
    ! 平均操作を行い z 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP)            :: aar_avr_aaz(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: kz
    
    ! 平均操作
    !
    do kz = kmin, kmax-1
!      aar_avr_aaz(:,:,kz) = &
!        &  (z_dz(kz)*aaz_Var(:,:,kz+1) + z_dz(kz+1)*aaz_Var(:,:,kz)) &
!        &  *0.5d0/r_dz(kz)
      aar_avr_aaz(:,:,kz) = (aaz_Var(:,:,kz+1) + aaz_Var(:,:,kz)) * 0.5d0
    end do
    
    ! kmax 格子上の値
    aar_avr_aaz(:,:,kmax) = 1.0d10
    
  end function aar_avr_aaz

  function paa_avr_xaa(xaa_Var)
    ! 半整数格子点の配列値を整数点上へ返す
    
    real(DP),intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: paa_avr_xaa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: ix

    ! 平均操作: x 座標から p 座標へ返す. imax の値は陽に代入する.
    !     
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !
    do ix = imin, imax-1
!      paa_avr_xaa(ix,:,:) = &
!        &  (x_dx(ix)*xaa_Var(ix+1,:,:) + x_dx(ix+1)*xaa_Var(ix,:,:)) &
!        &  *0.5d0/p_dx(ix)
      paa_avr_xaa(ix,:,:) = (xaa_Var(ix+1,:,:) + xaa_Var(ix,:,:)) * 0.5d0
    end do
    
    ! imax 格子上の値
    paa_avr_xaa(imax,:,:) = 1.0d10
    
  end function paa_avr_xaa

  function pyr_avr_xyz(xyz_Var)
    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax)     ! 入力
    real(DP)             :: pyr_avr_xyz(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: ix, kz                                     ! ループ添字
    
    do kz = kmin, kmax-1
      do ix = imin, imax-1
        pyr_avr_xyz(ix,:,kz) = ( xyz_Var(ix,:,kz)   + xyz_Var(ix+1,:,kz) + xyz_Var(ix,:,kz+1) + xyz_Var(ix+1,:,kz+1) ) * 0.25d0 
      end do
    end do

    ! imax 格子上の値
    pyr_avr_xyz(imax,:,:) = 1.0d10
    ! kmax 格子上の値
    pyr_avr_xyz(:,:,kmax) = 1.0d10

    
  end function pyr_avr_xyz

  function aya_avr_aqa(aqa_Var)
    ! 平均操作を行い y 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aqa_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aya_avr_aqa(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: jy
    
    ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
    if (jmin == jmax) then 
      aya_avr_aqa = aqa_Var
      return
    end if
    
    ! 平均操作
    !
    do jy = jmin+1, jmax
      aya_avr_aqa(:,jy,:) = (aqa_Var(:,jy,:) + aqa_Var(:,jy-1,:))*0.5d0 
    end do
    
    ! jmin 格子上の値
    aya_avr_aqa(:,jmin,:) = 1.0d10
    
  end function aya_avr_aqa

  function aaz_avr_aar(aar_Var)
    ! 平均操作を行い z 方向整数格子点の配列値を半整数格子点上へ返す
    
    real(DP),intent(in) :: aar_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aaz_avr_aar(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: kz
    
    ! 平均操作
    !
    do kz = kmin+1, kmax
      aaz_avr_aar(:,:,kz) = (aar_Var(:,:,kz) + aar_Var(:,:,kz-1))*0.5d0 
    end do
    
    ! kmin 格子上の値
    aaz_avr_aar(:,:,kmin) = 1.0d10
    
  end function aaz_avr_aar

  function paa_avr_xaa(xaa_Var)
    ! 半整数格子点の配列値を整数点上へ返す
    
    real(DP),intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: paa_avr_xaa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: ix

    ! 平均操作: x 座標から p 座標へ返す. imax の値は陽に代入する.
    !     
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !
    do ix = imin, imax-1
!      paa_avr_xaa(ix,:,:) = &
!        &  (x_dx(ix)*xaa_Var(ix+1,:,:) + x_dx(ix+1)*xaa_Var(ix,:,:)) &
!        &  *0.5d0/p_dx(ix)
      paa_avr_xaa(ix,:,:) = (xaa_Var(ix+1,:,:) + xaa_Var(ix,:,:)) * 0.5d0
    end do
    
    ! imax 格子上の値
    paa_avr_xaa(imax,:,:) = 1.0d10
    
  end function paa_avr_xaa

  function aa_AvrY_aya(aya_Var)
    ! aya 格子上の配列に対し y 方向平均を行う
    
    real(DP), intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrY_aya(imin:imax,kmin:kmax)       ! 出力
    
    ! 初期化
    aa_AvrY_aya = 0.0d0
    
    ! 平均
    aa_AvrY_aya = aa_IntY_aya(aya_Var)/sum(y_dy(1:ny))
    
  end function aa_AvrY_aya

  function aa_IntY_aya(aya_Var)
    ! aya 格子上の配列に対し y 方向積分を行う
    
    real(DP), intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntY_aya(imin:imax,kmin:kmax)       ! 出力
    integer                  :: ix, kz                           ! ループ添字
    
    ! 初期化
    aa_IntY_aya = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do ix = imin, imax
        aa_IntY_aya(ix,kz) = IntY_y(aya_Var(ix,:,kz))
      end do
    end do
    
  end function aa_IntY_aya

  function aa_AvrX_xaa(xaa_Var)
    ! xaa 格子上の配列に対し x 方向平均を行う

    real(DP), intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrX_xaa(jmin:jmax,kmin:kmax)       ! 出力
    
    ! 初期化
    aa_AvrX_xaa = 0.0d0
    
    ! 平均
    aa_AvrX_xaa = aa_IntX_xaa(xaa_Var)/sum(x_dx(1:nx))
    
  end function aa_AvrX_xaa

  function aa_IntX_xaa(xaa_Var)
    ! xaa 格子上の配列に対し x 方向積分を行う
    
    real(DP), intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntX_xaa(jmin:jmax,kmin:kmax)       ! 出力
    integer                  :: jy, kz                           ! ループ添字
    
    ! 初期化
    aa_IntX_xaa = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do jy = jmin, jmax
        aa_IntX_xaa(jy,kz) = IntX_x(xaa_Var(:,jy,kz))
      end do
    end do
    
  end function aa_IntX_xaa

  function a_AvrXY_xya(xya_Var)
    ! xya 格子上の配列に対し xy 方向平均を行う
    
    real(DP), intent(in) :: xya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: a_AvrXY_xya(kmin:kmax)       ! 出力

    ! 初期化
    a_AvrXY_xya = 0.0d0
    
    ! 平均
    a_AvrXY_xya = a_IntXY_xya(xya_Var)/(sum(x_dx(1:nx))*sum(y_dy(1:ny)))
    
  end function a_AvrXY_xya

  function a_IntXY_xya(xya_Var)
    ! xya 格子上の配列に対し xy 方向積分を行う
    
    real(DP), intent(in) :: xya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: a_IntXY_xya(kmin:kmax)       ! 出力
    integer                  :: kz                           ! ループ添字
    
    ! 初期化
    a_IntXY_xya = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      a_IntXY_xya(kz) = IntY_y(y_IntX_xy(xya_Var(:,:,kz)))
    end do
    
  end function a_IntXY_xya

  function aa_AvrZ_aaz(aaz_Var)
    ! aaz 格子上の配列に対し z 方向平均を行う
    
    real(DP), intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrZ_aaz(imin:imax,jmin:jmax)       ! 出力
    
    ! 初期化
    aa_AvrZ_aaz = 0.0d0
    
    ! 平均
    aa_AvrZ_aaz = aa_IntZ_aaz(aaz_Var)/sum(z_dz(1:nz))
    
  end function aa_AvrZ_aaz

  function aa_IntZ_aaz(aaz_Var)
    ! aaz 格子上の配列に対し z 方向積分を行う
    
    real(DP), intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntZ_aaz(imin:imax,jmin:jmax)       ! 出力
    integer                  :: ix, jy                           ! ループ添字
    
    ! 初期化
    aa_IntZ_aaz = 0.0d0

    ! 積分
    do jy = jmin, jmax
      do ix = imin, imax
        aa_IntZ_aaz(ix,jy) = IntZ_z(aaz_Var(ix,jy,:))
      end do
    end do
    
  end function aa_IntZ_aaz

  function aar_avr_aaz(aaz_Var)
    ! 平均操作を行い z 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP)            :: aar_avr_aaz(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: kz
    
    ! 平均操作
    !
    do kz = kmin, kmax-1
!      aar_avr_aaz(:,:,kz) = &
!        &  (z_dz(kz)*aaz_Var(:,:,kz+1) + z_dz(kz+1)*aaz_Var(:,:,kz)) &
!        &  *0.5d0/r_dz(kz)
      aar_avr_aaz(:,:,kz) = (aaz_Var(:,:,kz+1) + aaz_Var(:,:,kz)) * 0.5d0
    end do
    
    ! kmax 格子上の値
    aar_avr_aaz(:,:,kmax) = 1.0d10
    
  end function aar_avr_aaz

  function aqa_avr_aya(aya_Var)
    ! 平均操作を行い y 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: aqa_avr_aya(imin:imax,jmin:jmax,kmin:kmax) 
    integer                 :: jy
    
    ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
    if (jmin == jmax) then 
      aqa_avr_aya = aya_Var
      return
    end if
    
    ! 平均操作
    !
    do jy = jmin, jmax-1
!      aqa_avr_aya(:,jy,:) = &
!        &  (y_dy(jy)*aya_Var(:,jy+1,:) + y_dy(jy+1)*aya_Var(:,jy,:)) &
!        &  * 0.5d0/q_dy(jy)
      aqa_avr_aya(:,jy,:) = (aya_Var(:,jy+1,:) + aya_Var(:,jy,:)) * 0.5d0
    end do
    
    ! jmax 格子上の値
    aqa_avr_aya(:,jmax,:) = 1.0d10
    
  end function aqa_avr_aya

  function xqr_avr_xyz(xyz_Var)
    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax)       ! 入力
    real(DP)             :: xqr_avr_xyz(imin:imax,jmin:jmax,kmin:kmax)   ! 出力
    integer              :: jy, kz                                       ! ループ添字
    
    if (jmin == jmax) then 
      ! xr_avr_xz と同じになる
      do kz = kmin, kmax-1
        xqr_avr_xyz(:,:,kz) = ( xyz_Var(:,:,kz) + xyz_Var(:,:,kz+1) ) * 0.5d0 
      end do

      ! kmax 格子上の値
      xqr_avr_xyz(:,:,kmax) = 1.0d10
      
    else
      do kz = kmin, kmax-1
        do jy = jmin, jmax-1
          xqr_avr_xyz(:,jy,kz) = ( xyz_Var(:,jy,kz)   + xyz_Var(:,jy+1,kz) + xyz_Var(:,jy,kz+1) + xyz_Var(:,jy+1,kz+1) ) * 0.25d0 
        end do
      end do

      ! jmax 格子上の値
      xqr_avr_xyz(:,jmax,:) = 1.0d10
      ! kmax 格子上の値
      xqr_avr_xyz(:,:,kmax) = 1.0d10

    end if
    
  end function xqr_avr_xyz

  function xaa_avr_paa(paa_Var)
    ! 平均操作を行い x 方向整数格子点の配列値を半整数点上へ返す
    
    real(DP),intent(in) :: paa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xaa_avr_paa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: ix

    ! 平均操作: p 座標から x 座標へ返す. imin の値は陽に代入する.
    !
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !
    do ix = imin+1, imax
      xaa_avr_paa(ix,:,:) = (paa_Var(ix,:,:) + paa_Var(ix-1,:,:))*0.5d0 
    end do
    
    ! imin 格子上の値
    xaa_avr_paa(imin,:,:) = 1.0d10 
    
  end function xaa_avr_paa

  function aaz_avr_aar(aar_Var)
    ! 平均操作を行い z 方向整数格子点の配列値を半整数格子点上へ返す
    
    real(DP),intent(in) :: aar_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aaz_avr_aar(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: kz
    
    ! 平均操作
    !
    do kz = kmin+1, kmax
      aaz_avr_aar(:,:,kz) = (aar_Var(:,:,kz) + aar_Var(:,:,kz-1))*0.5d0 
    end do
    
    ! kmin 格子上の値
    aaz_avr_aar(:,:,kmin) = 1.0d10
    
  end function aaz_avr_aar

  function aqa_avr_aya(aya_Var)
    ! 平均操作を行い y 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: aqa_avr_aya(imin:imax,jmin:jmax,kmin:kmax) 
    integer                 :: jy
    
    ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
    if (jmin == jmax) then 
      aqa_avr_aya = aya_Var
      return
    end if
    
    ! 平均操作
    !
    do jy = jmin, jmax-1
!      aqa_avr_aya(:,jy,:) = &
!        &  (y_dy(jy)*aya_Var(:,jy+1,:) + y_dy(jy+1)*aya_Var(:,jy,:)) &
!        &  * 0.5d0/q_dy(jy)
      aqa_avr_aya(:,jy,:) = (aya_Var(:,jy+1,:) + aya_Var(:,jy,:)) * 0.5d0
    end do
    
    ! jmax 格子上の値
    aqa_avr_aya(:,jmax,:) = 1.0d10
    
  end function aqa_avr_aya

  function aa_AvrY_aya(aya_Var)
    ! aya 格子上の配列に対し y 方向平均を行う
    
    real(DP), intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrY_aya(imin:imax,kmin:kmax)       ! 出力
    
    ! 初期化
    aa_AvrY_aya = 0.0d0
    
    ! 平均
    aa_AvrY_aya = aa_IntY_aya(aya_Var)/sum(y_dy(1:ny))
    
  end function aa_AvrY_aya

  function aa_IntY_aya(aya_Var)
    ! aya 格子上の配列に対し y 方向積分を行う
    
    real(DP), intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntY_aya(imin:imax,kmin:kmax)       ! 出力
    integer                  :: ix, kz                           ! ループ添字
    
    ! 初期化
    aa_IntY_aya = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do ix = imin, imax
        aa_IntY_aya(ix,kz) = IntY_y(aya_Var(ix,:,kz))
      end do
    end do
    
  end function aa_IntY_aya

  function xy_AvrZ_xyr(xyr_Var)
    ! xyr 格子上の配列に対し z 方向平均を行う
    
    real(DP), intent(in) :: xyr_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: xy_AvrZ_xyr(imin:imax,jmin:jmax)       ! 出力
      
    ! 初期化
    xy_AvrZ_xyr = 0.0d0
    
    ! 平均
    xy_AvrZ_xyr = xy_IntZ_xyr(xyr_Var)/sum(z_dz(1:nz))
    
  end function xy_AvrZ_xyr

  function aa_AvrZ_aaz(aaz_Var)
    ! aaz 格子上の配列に対し z 方向平均を行う
    
    real(DP), intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrZ_aaz(imin:imax,jmin:jmax)       ! 出力
    
    ! 初期化
    aa_AvrZ_aaz = 0.0d0
    
    ! 平均
    aa_AvrZ_aaz = aa_IntZ_aaz(aaz_Var)/sum(z_dz(1:nz))
    
  end function aa_AvrZ_aaz

  function xy_IntZ_xyr(xyr_Var)
    ! xyr 格子上の配列に対し z 方向積分を行う
    
    real(DP), intent(in) :: xyr_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: xy_IntZ_xyr(imin:imax,jmin:jmax)       ! 出力
    integer                  :: ix, jy                           ! ループ添字
    
    ! 初期化
    xy_IntZ_xyr = 0.0d0
    
    ! 積分
    do jy = jmin, jmax
      do ix = imin, imax
        xy_IntZ_xyr(ix,jy) = IntZ_r(xyr_Var(ix,jy,:))
      end do
    end do
    
  end function xy_IntZ_xyr

  function aa_IntZ_aaz(aaz_Var)
    ! aaz 格子上の配列に対し z 方向積分を行う
    
    real(DP), intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntZ_aaz(imin:imax,jmin:jmax)       ! 出力
    integer                  :: ix, jy                           ! ループ添字
    
    ! 初期化
    aa_IntZ_aaz = 0.0d0

    ! 積分
    do jy = jmin, jmax
      do ix = imin, imax
        aa_IntZ_aaz(ix,jy) = IntZ_z(aaz_Var(ix,jy,:))
      end do
    end do
    
  end function aa_IntZ_aaz

  function xaa_avr_paa(paa_Var)
    ! 平均操作を行い x 方向整数格子点の配列値を半整数点上へ返す
    
    real(DP),intent(in) :: paa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xaa_avr_paa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: ix

    ! 平均操作: p 座標から x 座標へ返す. imin の値は陽に代入する.
    !
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !
    do ix = imin+1, imax
      xaa_avr_paa(ix,:,:) = (paa_Var(ix,:,:) + paa_Var(ix-1,:,:))*0.5d0 
    end do
    
    ! imin 格子上の値
    xaa_avr_paa(imin,:,:) = 1.0d10 
    
  end function xaa_avr_paa

  function aya_avr_aqa(aqa_Var)
    ! 平均操作を行い y 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aqa_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aya_avr_aqa(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: jy
    
    ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
    if (jmin == jmax) then 
      aya_avr_aqa = aqa_Var
      return
    end if
    
    ! 平均操作
    !
    do jy = jmin+1, jmax
      aya_avr_aqa(:,jy,:) = (aqa_Var(:,jy,:) + aqa_Var(:,jy-1,:))*0.5d0 
    end do
    
    ! jmin 格子上の値
    aya_avr_aqa(:,jmin,:) = 1.0d10
    
  end function aya_avr_aqa

  function aar_avr_aaz(aaz_Var)
    ! 平均操作を行い z 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP)            :: aar_avr_aaz(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: kz
    
    ! 平均操作
    !
    do kz = kmin, kmax-1
!      aar_avr_aaz(:,:,kz) = &
!        &  (z_dz(kz)*aaz_Var(:,:,kz+1) + z_dz(kz+1)*aaz_Var(:,:,kz)) &
!        &  *0.5d0/r_dz(kz)
      aar_avr_aaz(:,:,kz) = (aaz_Var(:,:,kz+1) + aaz_Var(:,:,kz)) * 0.5d0
    end do
    
    ! kmax 格子上の値
    aar_avr_aaz(:,:,kmax) = 1.0d10
    
  end function aar_avr_aaz

  function xyz_avr_pqz(pqz_Var)
    real(DP), intent(in) :: pqz_Var(imin:imax,jmin:jmax,kmin:kmax)     ! 入力
    real(DP)             :: xyz_avr_pqz(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: ix, jy                                     ! ループ添字
    
    if (jmin == jmax) then 
      ! xz_avr_pz と同じになる
      do ix = imin+1, imax
        xyz_avr_pqz(ix,:,:) = ( pqz_Var(ix-1,:,:) + pqz_Var(ix,:,:) ) * 0.5d0 
      end do

      ! imin 格子上の値
      xyz_avr_pqz(imin,:,:) = 1.0d10

    else
      do jy = jmin+1, jmax
        do ix = imin+1, imax
          xyz_avr_pqz(ix,jy,:) = ( pqz_Var(ix-1,jy-1,:) + pqz_Var(ix,jy-1,:) + pqz_Var(ix-1,jy,:)   + pqz_Var(ix,jy,:) ) * 0.25d0 
        end do
      end do

      ! imin 格子上の値
      xyz_avr_pqz(imin,:,:) = 1.0d10
      ! jmin 格子上の値
      xyz_avr_pqz(:,jmin,:) = 1.0d10
    end if
    
  end function xyz_avr_pqz

  function xyz_avr_pyr(pyr_Var)
    real(DP), intent(in) :: pyr_Var(imin:imax,jmin:jmax,kmin:kmax)     ! 入力
    real(DP)             :: xyz_avr_pyr(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: ix, kz                                     ! ループ添字
        
    do kz = kmin+1, kmax
      do ix = imin+1, imax
        xyz_avr_pyr(ix,:,kz) = ( pyr_Var(ix-1,:,kz-1) + pyr_Var(ix,:,kz-1) + pyr_Var(ix-1,:,kz)   + pyr_Var(ix,:,kz)     ) * 0.25d0 
      end do
    end do

    ! imin 格子上の値
    xyz_avr_pyr(imin,:,:) = 1.0d10
    ! kmin 格子上の値
    xyz_avr_pyr(:,:,kmin) = 1.0d10

  end function xyz_avr_pyr

  function xaa_avr_paa(paa_Var)
    ! 平均操作を行い x 方向整数格子点の配列値を半整数点上へ返す
    
    real(DP),intent(in) :: paa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xaa_avr_paa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: ix

    ! 平均操作: p 座標から x 座標へ返す. imin の値は陽に代入する.
    !
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !
    do ix = imin+1, imax
      xaa_avr_paa(ix,:,:) = (paa_Var(ix,:,:) + paa_Var(ix-1,:,:))*0.5d0 
    end do
    
    ! imin 格子上の値
    xaa_avr_paa(imin,:,:) = 1.0d10 
    
  end function xaa_avr_paa

  function xyz_avr_xqr(xqr_Var)
    real(DP), intent(in) :: xqr_Var(imin:imax,jmin:jmax,kmin:kmax)     ! 入力
    real(DP)             :: xyz_avr_xqr(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: jy, kz                                     ! ループ添字
        
    if (jmin == jmax) then  
      !xz_avr_xr と同じ
      do kz = kmin+1, kmax
        xyz_avr_xqr(:,:,kz) = ( xqr_Var(:,:,kz-1) + xqr_Var(:,:,kz) ) * 0.5d0 
      end do

      ! kmin 格子上の値
      xyz_avr_xqr(:,:,kmin) = 1.0d10

    else
      do kz = kmin+1, kmax
        do jy = jmin+1, jmax
          xyz_avr_xqr(:,jy,kz) = ( xqr_Var(:,jy-1,kz-1) + xqr_Var(:,jy,kz-1) + xqr_Var(:,jy-1,kz)   + xqr_Var(:,jy,kz)     ) * 0.25d0 
        end do
      end do

      ! jmin 格子上の値
      xyz_avr_xqr(:,jmin,:) = 1.0d10
      ! kmin 格子上の値
      xyz_avr_xqr(:,:,kmin) = 1.0d10
    end if
    
  end function xyz_avr_xqr

  function aya_avr_aqa(aqa_Var)
    ! 平均操作を行い y 方向半整数格子点の配列値を整数格子点上へ返す
    
    real(DP),intent(in) :: aqa_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aya_avr_aqa(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: jy
    
    ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
    if (jmin == jmax) then 
      aya_avr_aqa = aqa_Var
      return
    end if
    
    ! 平均操作
    !
    do jy = jmin+1, jmax
      aya_avr_aqa(:,jy,:) = (aqa_Var(:,jy,:) + aqa_Var(:,jy-1,:))*0.5d0 
    end do
    
    ! jmin 格子上の値
    aya_avr_aqa(:,jmin,:) = 1.0d10
    
  end function aya_avr_aqa

  function aaz_avr_aar(aar_Var)
    ! 平均操作を行い z 方向整数格子点の配列値を半整数格子点上へ返す
    
    real(DP),intent(in) :: aar_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aaz_avr_aar(imin:imax,jmin:jmax,kmin:kmax)
    integer                 :: kz
    
    ! 平均操作
    !
    do kz = kmin+1, kmax
      aaz_avr_aar(:,:,kz) = (aar_Var(:,:,kz) + aar_Var(:,:,kz-1))*0.5d0 
    end do
    
    ! kmin 格子上の値
    aaz_avr_aar(:,:,kmin) = 1.0d10
    
  end function aaz_avr_aar

  function xz_AvrY_xqz(xqz_Var)
    ! xqz 格子上の配列に対し y 方向平均を行う
    
    real(DP), intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: xz_AvrY_xqz(imin:imax,kmin:kmax)       ! 出力
    
    ! 初期化
    xz_AvrY_xqz = 0.0d0
    
    ! 平均
    xz_AvrY_xqz = xz_IntY_xqz(xqz_Var)/sum(y_dy(1:ny))
    
  end function xz_AvrY_xqz

  function aa_AvrY_aya(aya_Var)
    ! aya 格子上の配列に対し y 方向平均を行う
    
    real(DP), intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrY_aya(imin:imax,kmin:kmax)       ! 出力
    
    ! 初期化
    aa_AvrY_aya = 0.0d0
    
    ! 平均
    aa_AvrY_aya = aa_IntY_aya(aya_Var)/sum(y_dy(1:ny))
    
  end function aa_AvrY_aya

  function xz_IntY_xqz(xqz_Var)
    ! xqz 格子上の配列に対し y 方向積分を行う
    
    real(DP), intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: xz_IntY_xqz(imin:imax,kmin:kmax)       ! 出力
    integer                  :: ix, kz                           ! ループ添字
    
    ! 初期化
    xz_IntY_xqz = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do ix = imin, imax
        xz_IntY_xqz(ix,kz) = IntY_q(xqz_Var(ix,:,kz))
      end do
    end do
    
  end function xz_IntY_xqz

  function aa_IntY_aya(aya_Var)
    ! aya 格子上の配列に対し y 方向積分を行う
    
    real(DP), intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntY_aya(imin:imax,kmin:kmax)       ! 出力
    integer                  :: ix, kz                           ! ループ添字
    
    ! 初期化
    aa_IntY_aya = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do ix = imin, imax
        aa_IntY_aya(ix,kz) = IntY_y(aya_Var(ix,:,kz))
      end do
    end do
    
  end function aa_IntY_aya

  function aa_AvrX_xaa(xaa_Var)
    ! xaa 格子上の配列に対し x 方向平均を行う

    real(DP), intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrX_xaa(jmin:jmax,kmin:kmax)       ! 出力
    
    ! 初期化
    aa_AvrX_xaa = 0.0d0
    
    ! 平均
    aa_AvrX_xaa = aa_IntX_xaa(xaa_Var)/sum(x_dx(1:nx))
    
  end function aa_AvrX_xaa

  function aa_IntX_xaa(xaa_Var)
    ! xaa 格子上の配列に対し x 方向積分を行う
    
    real(DP), intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntX_xaa(jmin:jmax,kmin:kmax)       ! 出力
    integer                  :: jy, kz                           ! ループ添字
    
    ! 初期化
    aa_IntX_xaa = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do jy = jmin, jmax
        aa_IntX_xaa(jy,kz) = IntX_x(xaa_Var(:,jy,kz))
      end do
    end do
    
  end function aa_IntX_xaa

  function yz_AvrX_pyz(pyz_Var)
    ! pyz 格子上の配列に対し x 方向平均を行う
    
    real(DP), intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: yz_AvrX_pyz(jmin:jmax,kmin:kmax)       ! 出力
    
    ! 初期化
    yz_AvrX_pyz = 0.0d0
    
    ! 平均
    yz_AvrX_pyz = yz_IntX_pyz(pyz_Var)/sum(x_dx(1:nx))
    
  end function yz_AvrX_pyz

  function aa_AvrX_xaa(xaa_Var)
    ! xaa 格子上の配列に対し x 方向平均を行う

    real(DP), intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_AvrX_xaa(jmin:jmax,kmin:kmax)       ! 出力
    
    ! 初期化
    aa_AvrX_xaa = 0.0d0
    
    ! 平均
    aa_AvrX_xaa = aa_IntX_xaa(xaa_Var)/sum(x_dx(1:nx))
    
  end function aa_AvrX_xaa

  function yz_IntX_pyz(pyz_Var)
    ! pyz 格子上の配列に対し x 方向積分を行う
    
    real(DP), intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: yz_IntX_pyz(jmin:jmax,kmin:kmax)       ! 出力
    integer                  :: jy, kz                           ! ループ添字
    
    ! 初期化
    yz_IntX_pyz = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do jy = jmin, jmax
        yz_IntX_pyz(jy,kz) = IntX_p(pyz_Var(:,jy,kz))
      end do
    end do
    
  end function yz_IntX_pyz

  function aa_IntX_xaa(xaa_Var)
    ! xaa 格子上の配列に対し x 方向積分を行う
    
    real(DP), intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: aa_IntX_xaa(jmin:jmax,kmin:kmax)       ! 出力
    integer                  :: jy, kz                           ! ループ添字
    
    ! 初期化
    aa_IntX_xaa = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      do jy = jmin, jmax
        aa_IntX_xaa(jy,kz) = IntX_x(xaa_Var(:,jy,kz))
      end do
    end do
    
  end function aa_IntX_xaa

  function z_AvrXY_pyz(pyz_Var)
    ! pyz 格子上の配列に対し xy 方向平均を行う
    
    real(DP), intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: z_AvrXY_pyz(kmin:kmax)       ! 出力
    
    ! 初期化
    z_AvrXY_pyz = 0.0d0
    
    ! 平均
    z_AvrXY_pyz = z_IntXY_pyz(pyz_Var)/(sum(x_dx(1:nx))*sum(y_dy(1:ny)))
    
  end function z_AvrXY_pyz

  function z_AvrXY_xqz(xqz_Var)
    ! xqz 格子上の配列に対し xy 方向積分を行う
    
    real(DP), intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: z_AvrXY_xqz(kmin:kmax)       ! 出力
    
    ! 初期化
    z_AvrXY_xqz = 0.0d0
    
    ! 平均
    z_AvrXY_xqz = z_IntXY_xqz(xqz_Var)/(sum(x_dx(1:nx))*sum(y_dy(1:ny)))
    
  end function z_AvrXY_xqz

  function a_AvrXY_xya(xya_Var)
    ! xya 格子上の配列に対し xy 方向平均を行う
    
    real(DP), intent(in) :: xya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: a_AvrXY_xya(kmin:kmax)       ! 出力

    ! 初期化
    a_AvrXY_xya = 0.0d0
    
    ! 平均
    a_AvrXY_xya = a_IntXY_xya(xya_Var)/(sum(x_dx(1:nx))*sum(y_dy(1:ny)))
    
  end function a_AvrXY_xya

  function z_IntXY_pyz(pyz_Var)
    ! pyz 格子上の配列に対し xy 方向積分を行う
    
    real(DP), intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: z_IntXY_pyz(kmin:kmax)       ! 出力
    integer                  :: kz                           ! ループ添字
    
    ! 初期化
    z_IntXY_pyz = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      z_IntXY_pyz(kz) = IntY_y(y_IntX_py(pyz_Var(:,:,kz)))
    end do
    
  end function z_IntXY_pyz

  function z_IntXY_xqz(xqz_Var)
    ! xqz 格子上の配列に対し xy 方向積分を行う
    
    real(DP), intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: z_IntXY_xqz(kmin:kmax)       ! 出力
    integer                  :: kz                           ! ループ添字
    
    ! 初期化
    z_IntXY_xqz = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      z_IntXY_xqz(kz) = IntY_q(q_IntX_xq(xqz_Var(:,:,kz)))
    end do
    
  end function z_IntXY_xqz

  function a_IntXY_xya(xya_Var)
    ! xya 格子上の配列に対し xy 方向積分を行う
    
    real(DP), intent(in) :: xya_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: a_IntXY_xya(kmin:kmax)       ! 出力
    integer                  :: kz                           ! ループ添字
    
    ! 初期化
    a_IntXY_xya = 0.0d0
    
    ! 積分
    do kz = kmin, kmax
      a_IntXY_xya(kz) = IntY_y(y_IntX_xy(xya_Var(:,:,kz)))
    end do
    
  end function a_IntXY_xya