!= 3 次元 (xyz 方向) 等間隔交互格子 格子点設定モジュール
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: axesset.f90,v 1.16 2014/07/08 01:05:32 sugiyama Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]

module axesset
  != 3 次元 (xyz 方向) 等間隔交互格子 格子点設定モジュール
  !
  !== 概要
  !
  ! axesset は, 3 次元 (xyz 方向) 等間隔交互格子を用いた有限差分法に基づく
  ! 数値モデルのための, 基本的な Fortran90 副プログラムおよび関数を提供する. 
  ! 具体的に行っていることは以下の通り.
  !
  ! * 格子点座標配列と格子点間隔配列の設定
  ! * 格子点配列の平均関数
  ! 

  !モジュール読み込み
  !
  use dc_types,      only: DP
  use dc_iounit,     only: FileOpen
  use dc_message,    only: MessageNotify
  use mpi_wrapper,   only: myrank,               &
    &                      MPIWrapperCartCreate, &
    &                      MPIWrapperCartShift,  &
    &                      MPIWrapperCommFree
  use gridset,       only: FlagCalc3D,    &
    &                      xsub, ysub,    & ! 
    &                      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 :: xyz_pyz, xyr_pyr, xqz_pqz
  public :: pyz_xyz, pyr_xyr, pqz_xqz
  public :: xyz_xqz, pyz_pqz, xyr_xqr
  public :: xqz_xyz, pqz_pyz, xqr_xyr
  public :: xyz_xyr, pyz_pyr, xqz_xqr
  public :: xyr_xyz, pyr_pyz, xqr_xqz
  public :: pqz_xyz, pyr_xyz, xqr_xyz
  public :: xyz_pqz, xyz_pyr, xyz_xqr

  ! 変数定義
  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), public, save :: DX  
  real(DP), public, save :: DY           
  real(DP), public, save :: DZ           
  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 pqz_avr_xyz
    module procedure pqz_xyz
  end interface

  interface pyr_avr_xyz
    module procedure pyr_xyz
  end interface

  interface xqr_avr_xyz
    module procedure xqr_xyz
  end interface

  interface xyz_avr_pqz
    module procedure xyz_pqz
  end interface
  
  interface xyz_avr_pyr
    module procedure xyz_pyr
  end interface

  interface xyz_avr_xqr
    module procedure xyz_xqr
  end interface

!------------------------------------

  interface xyz_avr_pyz
    module procedure xaa_paa
  end interface
  
  interface xyr_avr_pyr
    module procedure xaa_paa
  end interface 

  interface xqz_avr_pqz
    module procedure xaa_paa
  end interface
  
  interface pyz_avr_xyz
    module procedure paa_xaa
  end interface 

  interface pqz_avr_xqz
    module procedure paa_xaa
  end interface 

  interface pyr_avr_xyr
    module procedure paa_xaa
  end interface 

  interface xyz_avr_xqz
    module procedure aya_aqa
  end interface 

  interface pyz_avr_pqz
    module procedure aya_aqa
  end interface 

  interface xyr_avr_xqr
    module procedure aya_aqa
  end interface
  
  interface xqz_avr_xyz
    module procedure aqa_aya
  end interface

  interface pqz_avr_pyz
    module procedure aqa_aya
  end interface 

  interface xqr_avr_xyr
    module procedure aqa_aya
  end interface 

  interface xyz_avr_xyr
    module procedure aaz_aar
  end interface 

  interface pyz_avr_pyr
    module procedure aaz_aar
  end interface 

  interface xqz_avr_xqr
    module procedure aaz_aar
  end interface 

  interface xyr_avr_xyz
    module procedure aar_aaz
  end interface 

  interface pyr_avr_pyz
    module procedure aar_aaz
  end interface 

  interface xqr_avr_xqz
    module procedure aar_aaz
  end interface 

!------------------------------------

  interface xyz_pyz
    module procedure xaa_paa
  end interface
  
  interface xyr_pyr
    module procedure xaa_paa
  end interface 

  interface xqz_pqz
    module procedure xaa_paa
  end interface
  
  interface pyz_xyz
    module procedure paa_xaa
  end interface 

  interface pqz_xqz
    module procedure paa_xaa
  end interface 

  interface pyr_xyr
    module procedure paa_xaa
  end interface 

  interface xyz_xqz
    module procedure aya_aqa
  end interface 

  interface pyz_pqz
    module procedure aya_aqa
  end interface 

  interface xyr_xqr
    module procedure aya_aqa
  end interface
  
  interface xqz_xyz
    module procedure aqa_aya
  end interface

  interface pqz_pyz
    module procedure aqa_aya
  end interface 

  interface xqr_xyr
    module procedure aqa_aya
  end interface 

  interface xyz_xyr
    module procedure aaz_aar
  end interface 

  interface pyz_pyr
    module procedure aaz_aar
  end interface 

  interface xqz_xqr
    module procedure aaz_aar
  end interface 

  interface xyr_xyz
    module procedure aar_aaz
  end interface 

  interface pyr_pyz
    module procedure aar_aaz
  end interface 

  interface xqr_xqz
    module procedure aar_aaz
  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     ! 設定ファイル用装置番号
    integer              :: comm_cart
    logical, parameter   :: periodic = .false.

    !設定ファイルから読み込む出力ファイル情報
    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 MPIWrapperCartCreate(xsub, ysub, periodic, comm_cart)
    call x_axis_init(comm_cart)
    call y_axis_init(comm_cart)
    call z_axis_init
    call MPIWrapperCommFree(comm_cart)
    
    ! 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/)    )
      call MessageNotify( "M", "axesset_init", "dx = %f", d=(/dx/)    )
      call MessageNotify( "M", "axesset_init", "dy = %f", d=(/dy/)    )
      call MessageNotify( "M", "axesset_init", "dz = %f", d=(/dz/)    )
    end if
  
  end subroutine axesset_init
  
  !--------------------------------------------------------------------
  subroutine x_axis_init(comm_cart)
    !
    != 格子点座標配列と格子点間隔配列の初期化    
    !
    
    ! 暗黙の型宣言禁止
    implicit none

    integer, intent(in) :: comm_cart

    ! 作業変数
    integer             :: ix   ! do ループ添字
    integer             :: rank, tmp
    integer             :: lp
    integer, parameter  :: direction = 1

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

    ! 初期化
    lp = -999

    ! 等間隔格子
    dx = (xmax - xmin) / ( nx * xsub )
    
    ! 隣り合う CPU の個数を得る
    sx: do ix = 1, xsub
      call MPIWrapperCartShift( comm_cart, direction, ix, rank, tmp )
      if ( rank < 0 ) then 
        lp = ix - 1
        exit sx
      end if
    end do sx

    !確認
    if ( lp == -999 ) then 
      call MessageNotify( "E", "axesset_init", "lp (x) is undifined" )
    end if

    ! 配列の作成
    do ix = imin, imax
      p_X(ix) = XMin + dx * nx * lp + dx * ix 
      x_X(ix) = XMin + dx * nx * lp + dx * ix - dx * 0.5d0
      x_dx(ix) = dx
      p_dx(ix) = dx
    end do
    
  end subroutine x_axis_init

  !--------------------------------------------------------------------
  subroutine y_axis_init(comm_cart)
    !
    != y 方向の座標値と格子点間隔を設定する
    !
    
    ! 暗黙の型宣言禁止
    implicit none

    integer, intent(in) :: comm_cart

    ! 作業変数
    integer             :: jy   ! ループ添字
    integer             :: lp
    integer             :: rank, tmp
    integer, parameter  :: direction = 0
    
    allocate(y_Y(jmin:jmax))
    allocate(q_Y(jmin:jmax))
    allocate(y_dy(jmin:jmax))
    allocate(q_dy(jmin:jmax))

    ! 初期化
    lp = -999

    ! 2 次元の場合は ymax = dx とする. 
    if (.NOT. FlagCalc3D ) then 
      ymax = dx
    end if
    
    ! 等間隔格子
    dy = (ymax - ymin) / ( ny * ysub )
    
    ! 隣り合う CPU の個数を得る
    sy: do jy = 1, ysub
      call MPIWrapperCartShift( comm_cart, direction, jy, tmp, rank )
      if ( rank < 0 ) then 
        lp = jy - 1
        exit sy
      end if
    end do sy

    !確認
    if ( lp == -999 ) then 
      call MessageNotify( "E", "axesset_init", "lp (y) is undifined" )
    end if
    
    ! 配列の作成
    do jy = jmin, jmax
      q_Y(jy) = YMin + dy * ny * lp  + dy * jy
      y_Y(jy) = YMin + dy * ny * lp  + dy * jy - dy * 0.5d0
      y_dy(jy) = dy
      q_dy(jy) = dy
    end do
    
  end subroutine y_axis_init
  !--------------------------------------------------------------------
  subroutine z_axis_init()
    !
    != z 方向の座標値と格子点間隔を設定する
    !
    
    ! 暗黙の型宣言禁止
    implicit none

    ! 作業変数
    integer                 :: kz   ! ループ添字
    
    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) = ZMin + dz * kz
      z_Z(kz) = ZMin + dz * (kz - 0.5)
      r_dz(kz) = dz
      z_dz(kz) = dz
    end do
    
  end subroutine z_axis_init
  
  !--------------------------------------------------------------------
  function xaa_paa(paa_Var)
    !
    != 平均操作を行い x 方向整数格子点の配列値を半整数点上へ返す 
    ! 
    ! 平均操作: p 座標から x 座標へ返す. imin の値は陽に代入する.
    !
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !

    ! 暗黙の型宣言禁止
    implicit none
        
    real(DP),intent(in) :: paa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xaa_paa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: i, j, k

    do k = kmin, kmax
      do j = jmin, jmax    
        do i = imin+1, imax
          xaa_paa(i,j,k) = ( paa_Var(i,j,k) + paa_Var(i-1,j,k) ) * 0.5d0 
        end do
      end do
    end do

    ! imin 格子上の値
    xaa_paa(imin,:,:) = 1.0d10 
    
  end function xaa_paa
  !--------------------------------------------------------------------
  function paa_xaa(xaa_Var)
    !
    != 半整数格子点の配列値を整数点上へ返す
    !
    ! 平均操作: x 座標から p 座標へ返す. imax の値は陽に代入する.
    !     
    !    p     imin    0     1         imax-1  imax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    x  imin   0      1         imax-1  imax
    !
    
    ! 暗黙の型宣言禁止
    implicit none

    real(DP),intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: paa_xaa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: i, j, k

    do k = kmin, kmax
      do j = jmin, jmax    
        do i = imin, imax-1
          paa_xaa(i,j,k) = ( xaa_Var(i+1,j,k) + xaa_Var(i,j,k) ) * 0.5d0
        end do
      end do
    end do
    
    ! imax 格子上の値
    paa_xaa(imax,:,:) = 1.0d10
    
  end function paa_xaa
  !--------------------------------------------------------------------
  function aya_aqa(aqa_Var)
    !
    != 平均操作を行い y 方向半整数格子点の配列値を整数格子点上へ返す
    !
    ! 平均操作: q 座標から y 座標へ返す. jmin の値は陽に代入する.
    !
    !    q     jmin    0     1         jmax-1  jmax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    y  jmin   0      1         jmax-1  jmax
    !
    
    ! 暗黙の型宣言禁止
    implicit none

    real(DP),intent(in) :: aqa_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aya_aqa(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: i, j, k 
 
    if ( FlagCalc3D ) then 
      do k = kmin, kmax
        do j = jmin+1, jmax    
          do i = imin,   imax
            aya_aqa(i,j,k) = ( aqa_Var(i,j,k) + aqa_Var(i,j-1,k) ) * 0.5d0 
          end do
        end do
      end do
      
      ! jmin 格子上の値
      aya_aqa(:,jmin,:) = 1.0d10

    else      
      ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
      aya_aqa = aqa_Var
      return
      
    end if   
    
  end function aya_aqa
  !--------------------------------------------------------------------
  function aqa_aya(aya_Var)
    !
    != 平均操作を行い y 方向の整数格子点の配列値を半整数格子点上へ返す
    !
    ! 平均操作: y 座標から q 座標へ返す. jmax の値は陽に代入する.
    !
    !    q     jmin    0     1         jmax-1  jmax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    y  jmin   0      1         jmax-1  jmax
    !
    
    ! 暗黙の型宣言禁止
    implicit none

    real(DP),intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: aqa_aya(imin:imax,jmin:jmax,kmin:kmax) 
    integer             :: i, j, k

    if ( FlagCalc3D ) then 
      do k = kmin, kmax
        do j = jmin, jmax-1  
          do i = imin, imax
            aqa_aya(i,j,k) = ( aya_Var(i,j+1,k) + aya_Var(i,j,k) ) * 0.5d0
          end do
        end do
      end do
      
      ! jmax 格子上の値
      aqa_aya(:,jmax,:) = 1.0d10

    else
    
      ! 2 次元計算(y 方向の配列要素数が 1)の場合, 代入値をそのまま返す
      aqa_aya = aya_Var
      return

    end if
    
    
  end function aqa_aya
  !--------------------------------------------------------------------
  function aaz_aar(aar_Var)
    !
    != 平均操作を行い z 方向整数格子点の配列値を半整数格子点上へ返す
    !
    ! 平均操作: r 座標から z 座標へ返す. zmin の値は陽に代入する.
    !
    !    r     kmin    0     1         kmax-1  kmax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    z  kmin   0      1         kmax-1  kmax
    !
    
    ! 暗黙の型宣言禁止
    implicit none

    real(DP),intent(in) :: aar_Var(imin:imax,jmin:jmax,kmin:kmax)   
    real(DP)            :: aaz_aar(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: i, j, k
    
    do k = kmin+1, kmax
      do j = jmin, jmax    
        do i = imin, imax
          aaz_aar(i,j,k) = ( aar_Var(i,j,k) + aar_Var(i,j,k-1) ) * 0.5d0 
        end do
      end do
    end do
    
    ! kmin 格子上の値
    aaz_aar(:,:,kmin) = 1.0d10
    
  end function aaz_aar
  !--------------------------------------------------------------------
  function aar_aaz(aaz_Var)
    !
    != 平均操作を行い z 方向半整数格子点の配列値を整数格子点上へ返す
    !
    ! 平均操作: z 座標から r 座標へ返す. zmax の値は陽に代入する.
    !
    !    r     kmin    0     1         kmax-1  kmax
    !     |--*--|--*--||--*--|--*--|--*--||--*--|
    !    z  kmin   0      1         kmax-1  kmax
    !
    
    ! 暗黙の型宣言禁止
    implicit none
    
    real(DP),intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP)            :: aar_aaz(imin:imax,jmin:jmax,kmin:kmax)
    integer             :: i, j, k
    
    do k = kmin, kmax-1
      do j = jmin, jmax    
        do i = imin, imax
          aar_aaz(i,j,k) = ( aaz_Var(i,j,k+1) + aaz_Var(i,j,k) ) * 0.5d0
        end do
      end do
    end do
    
    ! kmax 格子上の値
    aar_aaz(:,:,kmax) = 1.0d10
    
  end function aar_aaz
  !--------------------------------------------------------------------
  function pqz_xyz(xyz_Var)
    !
    != XY 平面において, 半整数格子 (xy) から整数格子 (pq) へ変換
    !

    ! 暗黙の型宣言禁止
    implicit none

    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax)  ! 入力
    real(DP)             :: pqz_xyz(imin:imax,jmin:jmax,kmin:kmax)  ! 出力
    integer              :: i, j, k                                 ! ループ添字

    if ( FlagCalc3D ) then 
      ! 3 次元の計算
      !
      do k = kmin, kmax
        do j = jmin, jmax-1
          do i = imin, imax-1
            pqz_xyz(i,j,k) =             &
              &  (                       &
              &     xyz_Var(i,   j,   k) &
              &   + xyz_Var(i+1, j,   k) &
              &   + xyz_Var(i,   j+1, k) &
              &   + xyz_Var(i+1, j+1, k) &
              &  ) * 0.25d0 
          end do
        end do
      end do
      
      ! imax 格子上の値
      pqz_xyz(imax,:,:) = 1.0d10
      ! jmax 格子上の値
      pqz_xyz(:,jmax,:) = 1.0d10
      
    else
      ! 2 次元計算 (y 方向の配列要素数が 1) の場合
      !
      do k = kmin, kmax
        do j = jmin, jmax    
          do i = imin, imax-1
            pqz_xyz(i,j,k) =           &
              &  (                     &
              &     xyz_Var(i,   j, k) &
              &   + xyz_Var(i+1, j, k) &
              &  ) * 0.5d0 
          end do
        end do
      end do

      ! imax 格子上の値
      pqz_xyz(imax,:,:) = 1.0d10
      
    end if
        
  end function pqz_xyz
  !--------------------------------------------------------------------
  function pyr_xyz(xyz_Var)
    !
    != XZ 平面において半整数格子 (xz) から整数格子 (pr) への変換を行う
    !

    ! 暗黙の型宣言禁止
    implicit none

    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: pyr_xyz(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: i, j, k                                ! ループ添字
    
    do k = kmin, kmax-1
      do j = jmin, jmax
        do i = imin, imax-1
          pyr_xyz(i,j,k) =                &
            &  (                          &
            &    + xyz_Var(i,   j, k  )   &
            &    + xyz_Var(i+1, j, k  )   &
            &    + xyz_Var(i,   j, k+1)   &
            &    + xyz_Var(i+1, j, k+1)   &
            &  ) * 0.25d0 
        end do
      end do
    end do

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

    
  end function pyr_xyz
  !--------------------------------------------------------------------
  function xqr_xyz(xyz_Var)
    !
    != YZ 平面において, 整数格子 (qr) から半整数格子 (yz) への変換を行う
    !

    ! 暗黙の型宣言禁止
    implicit none

    real(DP), intent(in) :: xyz_Var(imin:imax,jmin:jmax,kmin:kmax)   ! 入力
    real(DP)             :: xqr_xyz(imin:imax,jmin:jmax,kmin:kmax)   ! 出力
    integer              :: i, j, k                                  ! ループ添字
    
    if ( FlagCalc3D ) then 
      
      do k = kmin, kmax-1
        do j = jmin, jmax-1
          do i = imin, imax
            xqr_xyz(i,j,k) =              &
              &  (                        &
              &   + xyz_Var(i, j,   k  )  &
              &   + xyz_Var(i, j+1, k  )  &
              &   + xyz_Var(i, j,   k+1)  &
              &   + xyz_Var(i, j+1, k+1)  &
              &  ) * 0.25d0 
          end do
        end do
      end do

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

    else
      
      do k = kmin, kmax-1
        do j = jmin, jmax
          do i = imin, imax
            xqr_xyz(i,j,k) =            &
              &  (                      &
              &    + xyz_Var(i, j, k)   &
              &    + xyz_Var(i, j, k+1) &
              &  ) * 0.5d0 
          end do
        end do
      end do

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

    end if
    
  end function xqr_xyz
  !--------------------------------------------------------------------
  function xyz_pqz(pqz_Var)
    !
    != XY 平面において, 整数格子 (pq) から半整数格子 (xz) へ変換
    !

    ! 暗黙の型宣言禁止
    implicit none

    real(DP), intent(in) :: pqz_Var(imin:imax,jmin:jmax,kmin:kmax) ! 入力
    real(DP)             :: xyz_pqz(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: i, j, k                                ! ループ添字
    
    if ( FlagCalc3D ) then 
      ! 3 次元の計算
      !
      do k = kmin, kmax
        do j = jmin+1, jmax
          do i = imin+1, imax
            xyz_pqz(i,j,k) =               &
              &  (                         &
              &    + pqz_Var(i-1, j-1, k)  &
              &    + pqz_Var(i,   j-1, k)  &
              &    + pqz_Var(i-1, j,   k)  &
              &    + pqz_Var(i,   j,   k)  &
              &  ) * 0.25d0 
          end do
        end do
      end do

      ! imin 格子上の値
      xyz_pqz(imin,:,:) = 1.0d10
      ! jmin 格子上の値
      xyz_pqz(:,jmin,:) = 1.0d10
      
    else
      ! 2 次元計算 (y 方向の配列要素数が 1) の場合
      !
      do k = kmin, kmax
        do j = jmin, jmax
          do i = imin+1, imax
            xyz_pqz(i,j,k) =            &
              &  (                      &
              &    + pqz_Var(i-1, j, k) &
              &    + pqz_Var(i,   j, k) &
              &  ) * 0.5d0 
          end do
        end do
      end do
      
      ! imin 格子上の値
      xyz_pqz(imin,:,:) = 1.0d10

    end if
    
  end function xyz_pqz
  !--------------------------------------------------------------------
  function xyz_pyr(pyr_Var)
    !
    != XZ 平面において整数格子 (pr) から半整数格子 (xz) への変換を行う
    !

    ! 暗黙の型宣言禁止
    implicit none

    real(DP), intent(in) :: pyr_Var(imin:imax,jmin:jmax,kmin:kmax)     ! 入力
    real(DP)             :: xyz_pyr(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: i, j, k                                    ! ループ添字
        
    do k = kmin+1, kmax
      do j = jmin,   jmax
        do i = imin+1, imax
          xyz_pyr(i,j,k) =               &
            &  (                         &
            &    + pyr_Var(i-1, j, k-1)  &
            &    + pyr_Var(i,   j, k-1)  &
            &    + pyr_Var(i-1, j, k  )  &
            &    + pyr_Var(i,   j, k  )  &
            &  ) * 0.25d0 
        end do
      end do
    end do
    
    ! imin 格子上の値
    xyz_pyr(imin,:,:) = 1.0d10
    ! kmin 格子上の値
    xyz_pyr(:,:,kmin) = 1.0d10
    
  end function xyz_pyr
  !--------------------------------------------------------------------
  function xyz_xqr(xqr_Var)
    !
    != YZ 平面において, 半整数格子 (yz) から整数格子 (qr) への変換を行う
    !

    ! 暗黙の型宣言禁止
    implicit none

    real(DP), intent(in) :: xqr_Var(imin:imax,jmin:jmax,kmin:kmax)     ! 入力
    real(DP)             :: xyz_xqr(imin:imax,jmin:jmax,kmin:kmax) ! 出力
    integer              :: i, j, k                                    ! ループ添字

    if ( FlagCalc3D ) then 

      do k = kmin+1, kmax
        do j = jmin+1, jmax
          do i = imin, imax
            xyz_xqr(i,j,k) =               &
              &  (                         &
              &    + xqr_Var(i, j-1, k-1)  &
              &    + xqr_Var(i, j,   k-1)  &
              &    + xqr_Var(i, j-1, k  )  &
              &    + xqr_Var(i, j,   k  )  &
              &  ) * 0.25d0 
          end do
        end do
      end do

      ! jmin 格子上の値
      xyz_xqr(:,jmin,:) = 1.0d10
      ! kmin 格子上の値
      xyz_xqr(:,:,kmin) = 1.0d10
      
    else

      do k = kmin+1, kmax
        do j = jmin, jmax
          do i = imin, imax
            xyz_xqr(i,j,k) =              &
              &  (                        &
              &    + xqr_Var(i, j, k-1)   &
              &    + xqr_Var(i, j, k  )   &
              &  ) * 0.5d0 
          end do
        end do
      end do
      
      ! kmin 格子上の値
      xyz_xqr(:,:,kmin) = 1.0d10
      
    end if
    
  end function xyz_xqr

end module axesset