!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2006. All rights reserved.
!---------------------------------------------------------------------
!
!= 3 次元 (xyz 方向) 等間隔交互格子 格子点設定モジュール
!
!* 履歴
!  * 2007/07/15 (小高正嗣) : 3D deepconv へ移植, dc_types を Use.
!  * 2006/12/26 (小高正嗣) : 平均関数の名前を元に戻す
!  * 2006/06/05 (小高正嗣) : 新規作成
!
module xyz_base_module
  != 3 次元 (xyz 方向) 等間隔交互格子 格子点設定モジュール
  !
  !== 概要
  !
  ! xyz_base_module は, 3 次元 (xyz 方向) 等間隔交互格子を用いた有限差分
  ! 法に基づく数値モデルのための, 基本的な Fortran90 副プログラムおよび
  ! 関数を提供する. 
  !
  ! このモジュールは xyz_module の下位モジュールである. 下請けモジュール
  ! として data_type, x_base_module, y_base_module z_base_module 
  ! モジュールを用いている. 
  ! 
  !
  !== 備考
  !
  ! * 例外処理がない
  ! * その場合のメッセージ表示がない
  ! * 初期化の値はマシンイプシロン値としておくべき
  !
  use dc_types,  only : DP
  use gridset,   only : imin, imax, &
    &                   jmin, jmax, &
    &                   kmin, kmax
  use axesset,   only : dx, dy, dz 
  
  implicit none

  private

  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


  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 

  contains

    !--------------------------------------------------------------------
    !
    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                 :: jy, kz
      integer                 :: ix

      ! 初期化
      ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
      !
      xaa_avr_paa = 0.0d0

      ! 平均操作
      ! * 関数 x_avr_p を用いて計算
      !
!      do kz = kmin, kmax
!        do jy = jmin, jmax
!          xaa_avr_paa(:,jy,kz) = x_avr_p(paa_Var(:,jy,kz))
!        end do
!      end do

      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,:,:) = xaa_avr_paa(imin+1,:,:)


    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                 :: jy, kz                    
      integer                 :: ix

      ! 初期化
      ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
      !
!      paa_avr_xaa = 0.0d0

      ! 平均操作
      ! * 関数 p_avr_x を用いて計算
      !
!      do kz = kmin, kmax
!        do jy = jmin, jmax
!          paa_avr_xaa(:,jy,kz) = p_avr_x(xaa_Var(:,jy,kz))
!        end do
!      end do

      do ix = imin, imax-1
        paa_avr_xaa(ix,:,:) = &
          &  ( xaa_Var(ix+1,:,:) + xaa_Var(ix,:,:)) &
          &  *0.5d0
      end do

      paa_avr_xaa(imax,:,:) = paa_avr_xaa(imax-1,:,:) 

    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                 :: ix, kz
      integer                 :: jy

      ! 初期化
      ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
      !
!      aya_avr_aqa = 0.0d0

      ! 平均操作
      ! * 関数 y_q を用いて計算
      !
!      do kz = kmin, kmax
!        do ix = imin, imax
!          aya_avr_aqa(ix,:,kz) = y_avr_q(aqa_Var(ix,:,kz))
!        end do
!      end do

      do jy = jmin+1, jmax
        aya_avr_aqa(:,jy,:) = (aqa_Var(:,jy,:) + aqa_Var(:,jy-1,:))*0.5d0 
      end do
      
      aya_avr_aqa(:,jmin,:) = aya_avr_aqa(:,jmin+1,:)

    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                 :: ix, kz
      integer                 :: jy

      ! 初期化
      ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
      !
      aqa_avr_aya = 0.0d0

      ! 平均操作
      ! * 関数 q_avr_y を用いて計算
      !
!      do kz = kmin, kmax
!        do ix = imin, imax
!          aqa_avr_aya(ix,:,kz) = q_avr_y(aya_Var(ix,:,kz))
!        end do
!      end do

      do jy = jmin, jmax-1
        aqa_avr_aya(:,jy,:) = &
          &  ( aya_Var(:,jy+1,:) + aya_Var(:,jy,:)) &
          &  * 0.5d0
      end do
      
      aqa_avr_aya(:,jmax,:) =  aqa_avr_aya(:,jmax-1,:) 

    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                 :: ix, jy
      integer                 :: kz

      ! 初期化
      ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
      !
1      aaz_avr_aar = 0.0d0

      ! 平均操作
      ! * 関数 z_avr_r を用いて計算
      !
!      do jy = jmin, jmax
!        do ix = imin, imax
!          aaz_avr_aar(ix,jy,:) = z_avr_r(aar_Var(ix,jy,:))
!        end do
!      end do

      do kz = kmin+1, kmax
        aaz_avr_aar(:,:,kz) = (aar_Var(:,:,kz) + aar_Var(:,:,kz-1))*0.5d0 
      end do

      aaz_avr_aar(:,:,kmin) = aaz_avr_aar(:,:,kmin+1)

    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                 :: ix, jy
      integer                 :: kz

      ! 初期化
      ! * 0 割りを防ぐためにはマシンイプシロン値を用いるべき
      !
!      aar_avr_aaz = 0.0d0

      ! 平均操作
      ! * 関数 r_avr_z を用いて計算
      !
!      do jy = jmin, jmax
!        do ix = imin, imax
!          aar_avr_aaz(ix,jy,:) = r_avr_z(aaz_Var(ix,jy,:))
!        end do
!      end do

      do kz = kmin, kmax-1
        aar_avr_aaz(:,:,kz) = &
          &  (aaz_Var(:,:,kz+1) + aaz_Var(:,:,kz)) &
          &  *0.5d0
      end do

      aar_avr_aaz(:,:,kmax) = aar_avr_aaz(:,:,kmax-1)

    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                          ! ループ添字

      ! 初期化
      pqz_avr_xyz = 0.0d0

      do jy = jmin, jmax-1
        do ix = imin, imax-1
          pqz_avr_xyz(ix,jy,:) =        &
            &  (                        &
            &    + xyz_Var(ix+1,jy+1,:) &
            &    + xyz_Var(ix+1,jy,:)   &
            &    + xyz_Var(ix,jy+1,:)   &
            &    + xyz_Var(ix,jy,:)     &
            &  ) * 0.25d0 
        end do
      end do

!            &  ( xyz_Var(ix,jy,:)   + xyz_Var(ix+1,jy,:) +  &
!            &    xyz_Var(ix,jy+1,:) + xyz_Var(ix+1,jy+1,:) ) * 0.25d0 

    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                          ! ループ添字

      ! 初期化
      pyr_avr_xyz = 0.0d0

      do kz = kmin, kmax-1
        do ix = imin, imax-1
          pyr_avr_xyz(ix,:,kz) =        &
            &  (                        &
            &    + xyz_Var(ix+1,:,kz+1) &
            &    + xyz_Var(ix+1,:,kz)   &
            &    + xyz_Var(ix,:,kz+1)   &
            &    + xyz_Var(ix,:,kz)     &
            &  ) * 0.25d0 
        end do
      end do

!            &  ( xyz_Var(ix,:,kz)   + xyz_Var(ix+1,:,kz) +  &
!            &    xyz_Var(ix,:,kz+1) + xyz_Var(ix+1,:,kz+1) ) * 0.25d0 


    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                          ! ループ添字

      ! 初期化
      xqr_avr_xyz = 0.0d0

      do kz = kmin, kmax-1
        do jy = jmin, jmax-1
          xqr_avr_xyz(:,jy,kz) =          &
            &  (                          &
            &    + xyz_Var(:,jy+1,kz+1)   &
            &    + xyz_Var(:,jy+1,kz)     &
            &    + xyz_Var(:,jy,kz+1)     &
            &    + xyz_Var(:,jy,kz)       &
            &  ) * 0.25d0 
        end do
      end do

!            &  ( xyz_Var(:,jy,kz)   + xyz_Var(:,jy+1,kz) +  &
!            &    xyz_Var(:,jy,kz+1) + xyz_Var(:,jy+1,kz+1) ) * 0.25d0 

    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                          ! ループ添字

      ! 初期化
      xyz_avr_pqz = 0.0d0

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

!            &  ( pqz_Var(ix-1,jy-1,:) + pqz_Var(ix,jy-1,:) +  &
!            &    pqz_Var(ix-1,jy,:)   + pqz_Var(ix,jy,:)     ) * 0.25d0 

    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                          ! ループ添字

      ! 初期化
      xyz_avr_pyr = 0.0d0

      do kz = kmin+1, kmax
        do ix = imin+1, imax
          xyz_avr_pyr(ix,:,kz) =          &
            &  (                          &
            &    + pyr_Var(ix,  :,kz)     &
            &    + pyr_Var(ix,  :,kz-1)   &
            &    + pyr_Var(ix-1,:,kz)     &
            &    + pyr_Var(ix-1,:,kz-1)   &
            &  ) * 0.25d0 
        end do
      end do

!            &  ( pyr_Var(ix-1,:,kz-1) + pyr_Var(ix,:,kz-1) +  &
!            &    pyr_Var(ix-1,:,kz)   + pyr_Var(ix,:,kz)     ) * 0.25d0 

    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                          ! ループ添字

      ! 初期化
      xyz_avr_xqr = 0.0d0

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

!            &  ( xqr_Var(:,jy-1,kz-1) + xqr_Var(:,jy,kz-1) +  &
!            &    xqr_Var(:,jy-1,kz)   + xqr_Var(:,jy,kz)     ) * 0.25d0 

    end function xyz_avr_xqr
!--------------------------------------------------------------------
end module xyz_base_module
