!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2006. All rights reserved.
!---------------------------------------------------------------------
!
!= 3 次元 (xyz 方向) 不等間隔交互格子 微分演算モジュール
!
!* 履歴
!  * 2007/07/15 (小高正嗣) : 3D deepconv へ移植, dc_types を Use.
!  * 2007/06/27 (小高正嗣) : 1 次元微分演算モジュールに依存しないよう修正
!  * 2006/06/07 (小高正嗣) : 新規作成
!

module xyz_deriv_c4_module
  != 3 次元 (xyz 方向) 不等間隔交互格子 微分演算モジュール
  !
  !== 機能
  !
  ! xyz_deriv_module は, 3 次元 (xyz 方向) 不等間隔交互格子を用いた有限差分
  ! 法に基づく数値モデルのための, 微分演算を行う Fortran90 関数を提供する.
  ! 微分演算は 4 次精度中心差分を用いて行う.
  !
  ! このモジュールは下請けモジュールとして xyz_base_module, data_type
  ! モジュールを用いている.
  !
  ! モジュール内の変数と手続きの命名方法については xyz_module を参照のこと.
  !
  !
  !== 変数, 手続きの要約
  !
  !=== 1 階微分
  !
  ! xyz_dx_pyz, pyz_dx_xyz :: x 方向 1 階微分を計算する
  ! xyz_dy_xqz, xqz_dy_xyz :: y 方向 1 階微分を計算する
  ! xyz_dz_xyr, xyr_dz_xyz :: z 方向 1 階微分を計算する
  !
  !
  !== 備考
  !
  ! 差分式は格子間隔を一定とした
  !
  !    du/dx(i) = 9/8*[u(i)-u(i-1)]/dx(i) - 1/24*[u(i+1)-u(i-2)]/dx(i)
  !
  ! を用いている. 不等間隔格子の場合は精度が落ちることに注意.
  !

  use dc_types,        only : DBKIND => DP
  use dc_message,      only: MessageNotify
  use xyz_base_module, only : imin, imax, jmin, jmax, kmin, kmax, &
    &                         x_dx, p_dx, y_dy, q_dy, z_dz, r_dz

  implicit none
  private
  public :: xyz_dx_pyz, pyz_dx_xyz, pyr_dx_xyr, pqz_dx_xqz, &
    &       xyz_dy_xqz, xqz_dy_xyz, pqz_dy_pyz, xqr_dy_xyr, &
    &       xyz_dz_xyr, xyr_dz_xyz, pyr_dz_pyz, xqr_dz_xqz

  interface pyz_dx_xyz
    module procedure paa_dx_xaa
  end interface

  interface pyr_dx_xyr
    module procedure paa_dx_xaa
  end interface

  interface pqz_dx_xqz
    module procedure paa_dx_xaa
  end interface

  interface xqz_dy_xyz
    module procedure aqa_dy_aya
  end interface

  interface pqz_dy_pyz
    module procedure aqa_dy_aya
  end interface

  interface xqr_dy_xyr
    module procedure aqa_dy_aya
  end interface

  interface xyr_dz_xyz
    module procedure aar_dz_aaz
  end interface

  interface pyr_dz_pyz
    module procedure aar_dz_aaz
  end interface

  interface xqr_dz_xqz
    module procedure aar_dz_aaz
  end interface

  contains
!--------------------------------------------------------------------
    function xyz_dx_pyz(pyz_Var)
      ! 半整数格子点上の 1 階微分を計算する

      real(DBKIND),intent(in) :: pyz_Var(imin:imax,jmin:jmax,kmin:kmax)
      real(DBKIND)            :: xyz_dx_pyz(imin:imax,jmin:jmax,kmin:kmax)
      integer                 :: ix

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

      ! 1 階微分の計算
      !
      do ix = imin+2, imax-1
        xyz_dx_pyz(ix,:,:) = &
          &   ((pyz_Var(ix,:,:)   - pyz_Var(ix-1,:,:))*9.0d0/(8.0d0*x_dx(ix)) &
          & - ((pyz_Var(ix+1,:,:) - pyz_Var(ix-2,:,:))/(24.0d0*x_dx(ix))))
      end do

    end function xyz_dx_pyz
!--------------------------------------------------------------------
    function paa_dx_xaa(xaa_Var)
      ! 半整数格子点上の 1 階微分を計算する

      real(DBKIND),intent(in) :: xaa_Var(imin:imax,jmin:jmax,kmin:kmax)
      real(DBKIND)            :: paa_dx_xaa(imin:imax,jmin:jmax,kmin:kmax)
      integer                 :: ix

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

      ! 1 階微分の計算
      !
      do ix = imin+1, imax-2
        paa_dx_xaa(ix,:,:) = &
          &   ((xaa_Var(ix+1,:,:) - xaa_Var(ix,:,:))*9.0d0/(8.0d0*p_dx(ix)) &
          & - ((xaa_Var(ix+2,:,:) - xaa_Var(ix-1,:,:))/(24.0d0*p_dx(ix))))
      end do

    end function paa_dx_xaa
!--------------------------------------------------------------------
    function xyz_dy_xqz(xqz_Var)
      ! 半整数格子点上の 1 階微分を計算する

      real(DBKIND),intent(in) :: xqz_Var(imin:imax,jmin:jmax,kmin:kmax)
      real(DBKIND)            :: xyz_dy_xqz(imin:imax,jmin:jmax,kmin:kmax)
      integer                 :: jy

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

      ! 1 階微分の計算
      !
      do jy = jmin+2, jmax-1
        xyz_dy_xqz(:,jy,:) = &
          &   ((xqz_Var(:,jy,:)   - xqz_Var(:,jy-1,:))*9.0d0/(8.0d0*y_dy(jy)) &
          & - ((xqz_Var(:,jy+1,:) - xqz_Var(:,jy-2,:))/(24.0d0*y_dy(jy))))
      end do

    end function xyz_dy_xqz
!--------------------------------------------------------------------
    function aqa_dy_aya(aya_Var)
      ! 半整数格子点上の 1 階微分を計算する

      real(DBKIND),intent(in) :: aya_Var(imin:imax,jmin:jmax,kmin:kmax)
      real(DBKIND)            :: aqa_dy_aya(imin:imax,jmin:jmax,kmin:kmax)
      integer                 :: jy

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

      ! 1 階微分の計算
      !
      do jy = jmin+1, jmax-2
        aqa_dy_aya(:,jy,:) = &
          &   ((aya_Var(:,jy+1,:) - aya_Var(:,jy,:))*9.0d0/(8.0d0*q_dy(jy)) &
          & - ((aya_Var(:,jy+2,:) - aya_Var(:,jy-1,:))/(24.0d0*q_dy(jy))))
      end do

    end function aqa_dy_aya
!--------------------------------------------------------------------
    function xyz_dz_xyr(xyr_Var)
      ! 半整数格子点上の 1 階微分を計算する

      real(DBKIND),intent(in) :: xyr_Var(imin:imax,jmin:jmax,kmin:kmax)
      real(DBKIND)            :: xyz_dz_xyr(imin:imax,jmin:jmax,kmin:kmax)
      integer                 :: kz

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

      ! 1 階微分の計算
      !
      do kz = kmin+2, kmax-1
        xyz_dz_xyr(:,:,kz) = &
          &   ((xyr_Var(:,:,kz)   - xyr_Var(:,:,kz-1))*9.0d0/(8.0d0*z_dz(kz)) &
          & - ((xyr_Var(:,:,kz+1) - xyr_Var(:,:,kz-2))/(24.0d0*z_dz(kz))))
      end do

    end function xyz_dz_xyr
!--------------------------------------------------------------------
    function aar_dz_aaz(aaz_Var)
      ! 半整数格子点上の 1 階微分を計算する

      real(DBKIND),intent(in) :: aaz_Var(imin:imax,jmin:jmax,kmin:kmax)
      real(DBKIND)            :: aar_dz_aaz(imin:imax,jmin:jmax,kmin:kmax)
      integer                 :: kz

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

      ! 1 階微分の計算
      ! * r_dz_z を用いて計算
      !
      do kz = kmin+1, kmax-2
        aar_dz_aaz(:,:,kz) = &
          &   ((aaz_Var(:,:,kz+1) - aaz_Var(:,:,kz))*9.0d0/(8.0d0*r_dz(kz)) &
          & - ((aaz_Var(:,:,kz+2) - aaz_Var(:,:,kz-1))/(24.0d0*r_dz(kz))))
      end do

    end function aar_dz_aaz
!--------------------------------------------------------------------
end module xyz_deriv_c4_module