!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2006. All rights reserved.
!---------------------------------------------------------------------
!
!= 3 次元 (xyz 方向) 等間隔交互格子 格子点設定モジュール
!
!* 履歴
!  * 2010/08/02 (小高正嗣) : コメント追加
!  * 2007/07/15 (小高正嗣) : 3D deepconv へ移植, dc_types を Use.
!  * 2006/12/26 (小高正嗣) : 平均関数の名前を元に戻す
!  * 2006/06/05 (小高正嗣) : 新規作成
!
module axesset
  != 3 次元 (xyz 方向) 等間隔交互格子 格子点設定モジュール
  !
  !== 概要
  !
  ! axesset は, 3 次元 (xyz 方向) 等間隔交互格子を用いた有限差分法に基づく
  ! 数値モデルのための, 基本的な Fortran90 副プログラムおよび関数を提供する. 
  ! 具体的に行っていることは以下の通り.
  !
  ! * 格子点座標配列と格子点間隔配列の設定
  ! * 格子点配列の補間関数の設定
  ! * 格子点配列の積分・平均関数
  ! 
  !== 備考
  !
  ! * 例外処理がない
  ! * その場合のメッセージ表示がない
  ! * 初期化の値はマシンイプシロン値としておくべき
  !

  !モジュール読み込み
  use dc_types,    only: DP, STRING
  use dc_iounit,   only: FileOpen
  use dc_message,  only: MessageNotify
  use mpi_wrapper, only: myrank, nprocs
  use gridset,     only: &
    &                  nx, ny, nz,    & ! 格子点数
    &                  imin, imax, jmin, jmax, kmin, kmax ! 配列の上限値と下限値
  use namelist_util, only: namelist_filename

  !暗黙の型宣言禁止
  implicit none

  !デフォルトは非公開
  private

  !公開要素
  ! 初期化手続き
  public :: axesset_init
!!
  ! 格子点配列の補間関数 
  public :: xyz_avr_pyz, xyr_avr_pyr, xqz_avr_pqz
  public :: pyz_avr_xyz, pyr_avr_xyr, pqz_avr_xqz
  public :: xyz_avr_xqz, pyz_avr_pqz, xyr_avr_xqr
  public :: xqz_avr_xyz, pqz_avr_pyz, xqr_avr_xyr
  public :: xyz_avr_xyr, pyz_avr_pyr, xqz_avr_xqr
  public :: xyr_avr_xyz, pyr_avr_pyz, xqr_avr_xqz
  public :: pqz_avr_xyz, pyr_avr_xyz, xqr_avr_xyz
  public :: xyz_avr_pqz, xyz_avr_pyr, xyz_avr_xqr
!!
  ! 格子点配列の積分・平均関数
  public :: yz_IntX_xyz, yz_IntX_pyz, qz_IntX_xqz, yr_IntX_xyr
  public :: xz_IntY_xyz, xz_IntY_xqz, pz_IntY_pyz, xr_IntY_xyr 
  public :: xy_IntZ_xyz, xy_IntZ_xyr, py_IntZ_pyz, xq_IntZ_xqz 
  public :: z_IntXY_xyz, z_IntXY_pyz, z_IntXY_xqz, r_IntXY_xyr
  public :: IntXYZ_xyz , IntXYZ_pyz , IntXYZ_xqz , IntXYZ_xyr
  public :: yz_AvrX_xyz, yz_AvrX_pyz, qz_AvrX_xqz, yr_AvrX_xyr
  public :: xz_AvrY_xyz, xz_AvrY_xqz, pz_AvrY_pyz, xr_AvrY_xyr 
  public :: xy_AvrZ_xyz, xy_AvrZ_xyr, py_AvrZ_pyz, xq_AvrZ_xqz 
  public :: z_AvrXY_xyz, z_AvrXY_pyz, z_AvrXY_xqz, r_AvrXY_xyr
  public :: AvrXYZ_xyz , AvrXYZ_pyz , AvrXYZ_xqz , AvrXYZ_xyr

  ! 変数定義
  real(DP), public, save :: Xmin = 0.0d0   ! x 座標の始点・終点
  real(DP), public, save :: Xmax = 1.0d4   ! x 座標の始点・終点
  real(DP), public, save :: Ymin = 0.0d0   ! x 座標の始点・終点
  real(DP), public, save :: Ymax = 1.0d4   ! x 座標の始点・終点
  real(DP), public, save :: Zmin = 0.0d0   ! z 座標の始点・終点
  real(DP), public, save :: Zmax = 1.0d4   ! z 座標の始点・終点
  real(DP), allocatable, public, save :: x_X(:)         ! 半整数格子点座標
  real(DP), allocatable, public, save :: p_X(:)         ! 整数格子点座標
  real(DP), allocatable, public, save :: x_dx(:)        ! 半整数格子点間隔
  real(DP), allocatable, public, save :: p_dx(:)        ! 整数格子点間隔
  real(DP), allocatable, public, save :: y_Y(:)         ! 半整数格子点座標
  real(DP), allocatable, public, save :: q_Y(:)         ! 整数格子点座標
  real(DP), allocatable, public, save :: y_dy(:)        ! 半整数格子点間隔
  real(DP), allocatable, public, save :: q_dy(:)        ! 整数格子点間隔
  real(DP), allocatable, public, save :: z_Z(:)         ! 半整数格子点座標
  real(DP), allocatable, public, save :: r_Z(:)         ! 整数格子点座標
  real(DP), allocatable, public, save :: z_dz(:)        ! 半整数格子点間隔
  real(DP), allocatable, public, save :: r_dz(:)        ! 整数格子点間隔
  real(DP), allocatable, public, save :: xyz_X(:,:,:)   ! x 座標(半整数格子)
  real(DP), allocatable, public, save :: xyz_Y(:,:,:)   ! y 座標(半整数格子)
  real(DP), allocatable, public, save :: xyz_Z(:,:,:)   ! z 座標(半整数格子)
  real(DP), allocatable, public, save :: xyz_dX(:,:,:)  ! x 格子間隔(半整数格子)
  real(DP), allocatable, public, save :: xyz_dY(:,:,:)  ! y 格子間隔(半整数格子)
  real(DP), allocatable, public, save :: xyz_dZ(:,:,:)  ! z 格子間隔(半整数格子)


  interface y_IntX_xy
    module procedure a_IntX_xa
  end interface 

  interface q_IntX_xq
    module procedure a_IntX_xa
  end interface 

  interface xyz_avr_pyz
    module procedure xaa_avr_paa
  end interface
  
  interface xyr_avr_pyr
    module procedure xaa_avr_paa
  end interface 

  interface xqz_avr_pqz
    module procedure xaa_avr_paa
  end interface
  
  interface pyz_avr_xyz
    module procedure paa_avr_xaa
  end interface 

  interface pqz_avr_xqz
    module procedure paa_avr_xaa
  end interface 

  interface pyr_avr_xyr
    module procedure paa_avr_xaa
  end interface 

  interface xyz_avr_xqz
    module procedure aya_avr_aqa
  end interface 

  interface pyz_avr_pqz
    module procedure aya_avr_aqa
  end interface 

  interface xyr_avr_xqr
    module procedure aya_avr_aqa
  end interface
  
  interface xqz_avr_xyz
    module procedure aqa_avr_aya
  end interface

  interface pqz_avr_pyz
    module procedure aqa_avr_aya
  end interface 

  interface xqr_avr_xyr
    module procedure aqa_avr_aya
  end interface 

  interface xyz_avr_xyr
    module procedure aaz_avr_aar
  end interface 

  interface pyz_avr_pyr
    module procedure aaz_avr_aar
  end interface 

  interface xqz_avr_xqr
    module procedure aaz_avr_aar
  end interface 

  interface xyr_avr_xyz
    module procedure aar_avr_aaz
  end interface 

  interface pyr_avr_pyz
    module procedure aar_avr_aaz
  end interface 

  interface xqr_avr_xqz
    module procedure aar_avr_aaz
  end interface 

  interface yz_IntX_xyz
    module procedure aa_IntX_xaa
  end interface

  interface qz_IntX_xqz
    module procedure aa_IntX_xaa
  end interface

  interface yr_IntX_xyr
    module procedure aa_IntX_xaa
  end interface

  interface xz_IntY_xyz
    module procedure aa_IntY_aya
  end interface

  interface pz_IntY_pyz
    module procedure aa_IntY_aya
  end interface

  interface xr_IntY_xyr
    module procedure aa_IntY_aya
  end interface

  interface xy_IntZ_xyz
    module procedure aa_IntZ_aaz
  end interface

  interface py_IntZ_pyz
    module procedure aa_IntZ_aaz
  end interface

  interface xq_IntZ_xqz
    module procedure aa_IntZ_aaz
  end interface

  interface z_IntXY_xyz
    module procedure a_IntXY_xya
  end interface

  interface r_IntXY_xyr
    module procedure a_IntXY_xya
  end interface

  interface yz_AvrX_xyz
    module procedure aa_AvrX_xaa
  end interface

  interface qz_AvrX_xqz
    module procedure aa_AvrX_xaa
  end interface

  interface yr_AvrX_xyr
    module procedure aa_AvrX_xaa
  end interface

  interface xz_AvrY_xyz
    module procedure aa_AvrY_aya
  end interface

  interface pz_AvrY_pyz
    module procedure aa_AvrY_aya
  end interface

  interface xr_AvrY_xyr
    module procedure aa_AvrY_aya
  end interface

  interface xy_AvrZ_xyz
    module procedure aa_AvrZ_aaz
  end interface

  interface py_AvrZ_pyz
    module procedure aa_AvrZ_aaz
  end interface

  interface xq_AvrZ_xqz
    module procedure aa_AvrZ_aaz
  end interface

  interface z_AvrXY_xyz
    module procedure a_AvrXY_xya
  end interface

  interface r_AvrXY_xyr
    module procedure a_AvrXY_xya
  end interface

contains
  !--------------------------------------------------------------------
  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
  
  !--------------------------------------------------------------------
  subroutine x_axis_init()
    ! 配列の上下限の値と格子点座標値, 格子点間隔を設定する
    
    integer                 :: ix   ! do ループ添字
    real(DP)                :: dx   ! 格子間隔

    allocate(x_X(imin:imax))
    allocate(p_X(imin:imax))
    allocate(x_dx(imin:imax))
    allocate(p_dx(imin:imax))

    ! 等間隔格子
    dx = (xmax - xmin) / nx
    
    do ix = imin, imax
      p_X(ix) = dx * ix  + myrank * xmax
      x_X(ix) = dx * (ix - 0.5) + myrank * xmax
!      p_X(ix) = dx * ix  
!      x_X(ix) = dx * (ix - 0.5)
      x_dx(ix) = dx
      p_dx(ix) = dx
    end do
    
  end subroutine x_axis_init
  !--------------------------------------------------------------------
  subroutine y_axis_init()
    ! y 方向の座標値と格子点間隔を設定する
    
    integer                 :: jy   ! ループ添字
    real(DP)                :: dy
    
    allocate(y_Y(jmin:jmax))
    allocate(q_Y(jmin:jmax))
    allocate(y_dy(jmin:jmax))
    allocate(q_dy(jmin:jmax))
    
    ! 等間隔格子
    dy = (ymax - ymin) / ny
    
    do jy = jmin, jmax
      q_Y(jy) = dy * jy
      y_Y(jy) = dy * (jy - 0.5)
      y_dy(jy) = dy
      q_dy(jy) = dy
    end do
    
  end subroutine y_axis_init
  !--------------------------------------------------------------------
  subroutine z_axis_init()
    ! z 方向の座標値と格子点間隔を設定する
    
    integer                 :: kz   ! ループ添字
    real(DP)            :: dz
    
    allocate(z_Z(kmin:kmax))
    allocate(r_Z(kmin:kmax))
    allocate(z_dz(kmin:kmax))
    allocate(r_dz(kmin:kmax))
    
    ! 等間隔格子
    dz = (zmax - zmin) / nz
    
    do kz = kmin, kmax
      r_Z(kz) = dz * kz
      z_Z(kz) = dz * (kz - 0.5)
      r_dz(kz) = dz
      z_dz(kz) = dz
    end do
    
  end subroutine z_axis_init
  
  !--------------------------------------------------------------------
  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 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)
    end do
    
    ! imax 格子上の値
    paa_avr_xaa(imax,:,:) = 1.0d10
    
  end function paa_avr_xaa
  !--------------------------------------------------------------------
  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 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)
    end do
    
    ! jmax 格子上の値
    aqa_avr_aya(:,jmax,:) = 1.0d10
    
  end function aqa_avr_aya
  !--------------------------------------------------------------------
  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 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)
    end do
    
    ! kmax 格子上の値
    aar_avr_aaz(:,:,kmax) = 1.0d10
    
  end function aar_avr_aaz
  !--------------------------------------------------------------------
  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 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 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 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 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 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_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_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 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_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 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 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 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 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
  !--------------------------------------------------------------------
  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 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_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_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_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_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_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 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 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 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_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 IntX_p(p_Var)
    ! 整数格子上の配列に対し x 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: p_Var(imin:imax)  ! 入力
    real(DP)             :: IntX_p            ! 出力
    
    ! 初期化
    IntX_p = 0.0d0
    
    ! 積分
    IntX_p = sum(p_Var(1:nx)*p_dx(1:nx))
    
  end function IntX_p
  !--------------------------------------------------------------------
  function IntX_x(x_Var)
    ! 半整数格子上の配列に対し x 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: x_Var(imin:imax)  ! 入力
    real(DP)             :: IntX_x            ! 出力
    real(DP), allocatable:: work(:)
    
    ! 初期化
    IntX_x = 0.0d0
    allocate(work(imin:imax))
    
    ! 積分
    work = x_Var*x_dx
    
    IntX_x = sum(work(1:nx))
    
    deallocate(work)
    
  end function IntX_x
!--------------------------------------------------------------------
  function IntY_q(q_Var)
    ! 整数格子上の配列に対し y 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: q_Var(jmin:jmax)  ! 入力
    real(DP)             :: IntY_q            ! 出力
    
    ! 初期化
    IntY_q = 0.0d0
    
    ! 積分
    IntY_q = sum(q_Var(1:ny)*q_dy(1:ny))
    
  end function IntY_q
!--------------------------------------------------------------------
  function IntY_y(y_Var)
    ! 半整数格子上の配列に対し y 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: y_Var(jmin:jmax)  ! 入力
    real(DP)             :: IntY_y            ! 出力
    
    ! 初期化
    IntY_y = 0.0d0
    
    ! 積分
    IntY_y = sum(y_Var(1:ny)*y_dy(1:ny))
    
  end function IntY_y
  !--------------------------------------------------------------------
  function IntZ_r(r_Var)
    ! 整数格子上の配列に対し z 方向に重み付きの積分を行う

    real(DP), intent(in) :: r_Var(kmin:kmax)  ! 入力
    real(DP)             :: IntZ_r            ! 出力
    
    ! 初期化
    IntZ_r = 0.0d0
    
    ! 積分
    IntZ_r = sum(r_Var(1:nz)*r_dz(1:nz))
    
  end function IntZ_r
!--------------------------------------------------------------------
  function IntZ_z(z_Var)
    ! 半整数格子上の配列に対し z 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: z_Var(kmin:kmax)  ! 入力
    real(DP)             :: IntZ_z            ! 出力

    ! 初期化
    IntZ_z = 0.0d0
    
    ! 積分
    IntZ_z = sum(z_Var(1:nz)*z_dz(1:nz))
    
  end function IntZ_z
!--------------------------------------------------------------------
  function a_IntX_xa(xa_Var)
    ! xa 格子上の配列に対し x 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: xa_Var(imin:imax,jmin:jmax) ! 入力
    real(DP)             :: a_IntX_xa(jmin:jmax)        ! 出力
    integer                  :: jy                          ! ループ添字

    ! 初期化
    a_IntX_xa = 0.0d0
    
    ! 積分
    do jy = jmin, jmax
      a_IntX_xa(jy) = IntX_x(xa_Var(:,jy))
    end do
    
  end function a_IntX_xa
!--------------------------------------------------------------------
  function y_IntX_py(py_Var)
    ! py 格子上の配列に対し x 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: py_Var(imin:imax,jmin:jmax) ! 入力
    real(DP)             :: y_IntX_py(jmin:jmax)        ! 出力
    integer                  :: jy                          ! ループ添字
    
    ! 初期化
    y_IntX_py = 0.0d0

    ! 積分
    do jy = jmin, jmax
      y_IntX_py(jy) = IntX_p(py_Var(:,jy))
    end do
    
  end function y_IntX_py
  !--------------------------------------------------------------------
  function a_IntY_ay(ay_Var)
    ! ay 格子上の配列に対し y 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: ay_Var(imin:imax,jmin:jmax) ! 入力
    real(DP)             :: a_IntY_ay(imin:imax)        ! 出力
    integer                  :: ix                          ! ループ添字
    
    ! 初期化
    a_IntY_ay = 0.0d0
    
    ! 積分
    do ix = imin, imax
      a_IntY_ay(ix) = IntY_y(ay_Var(ix,:))
    end do
    
  end function a_IntY_ay
  !--------------------------------------------------------------------
  function x_IntY_xq(xq_Var)
    ! ay 格子上の配列に対し y 方向に重み付きの積分を行う
    
    real(DP), intent(in) :: xq_Var(imin:imax,jmin:jmax) ! 入力
    real(DP)             :: x_IntY_xq(imin:imax)        ! 出力
    integer                  :: ix                          ! ループ添字
    
    ! 初期化
    x_IntY_xq = 0.0d0
    
    ! 積分
    do ix = imin, imax
      x_IntY_xq(ix) = IntY_q(xq_Var(ix,:))
    end do
    
  end function x_IntY_xq

end module axesset

