subroutine physics_interpolate_geopot( xyz_Press, xyr_Press, xyz_GeoPot, xyr_GeoPot, xy_Press, xyz_Temp, xyr_Temp, z_Sigma, r_Sigma                               )  ! (in)

    !==== Dependency
    use type_mod,     only : REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod,   only: im, jm, km
    use constants_mod, only: constants_init, RAir, Grav !!$ , Cp
    use dc_trace,     only : SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Input
    !
    real(DBKIND), intent(in) :: xy_Press(im*jm), xyz_Temp(im*jm,km), xyr_Temp(im*jm,km+1), z_Sigma(:), r_Sigma(:)              ! intent(in): σレベル(半整数)座標

    !==== Output
    !
    real(DBKIND), intent(out) :: xyz_Press(im*jm,km), xyr_Press(im*jm,km+1), xyz_GeoPot(im*jm,km), xyr_GeoPot(im*jm,km+1)   ! intent(out): ジオポテンシャル(半整数)座標

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_interpolate_geopot"
    real(DBKIND)        :: z_DelSigma(km), r_DelSigma(km+1)        ! Δσ(半整数レベル)
    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: ij, k

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   気圧演算 (units:[Pa])
    !----------------------------------------------------------------

    !----- 整数座標 -----
    do k = 1, km
       do ij = 1, im*jm
          xyz_Press(ij,k) = xy_Press(ij) * z_Sigma(k) 
       end do
    end do

    !----- 半整数座標 -----
    do k = 1, km+1
       do ij = 1, im*jm
          xyr_Press(ij,k) = xy_Press(ij) * r_Sigma(k) 
       end do
    end do

    !----------------------------------------------------------------
    !   高度演算 (units:[m]???, よ, よく分かんない…)
    !----------------------------------------------------------------

    !----- Δσの計算 -----
    do k = 1, km
       z_DelSigma(k) = r_Sigma(k) - r_Sigma(k+1)
    end do

    do k = 2, km
       r_DelSigma(k) = z_Sigma(k-1) - z_Sigma(k)
    end do
    r_DelSigma(1)    = r_Sigma(1)  - z_Sigma(1)
    r_DelSigma(km+1) = z_Sigma(km) - r_Sigma(km+1)

    !----- 整数座標 -----
    do ij = 1, im*jm
       xyz_GeoPot(ij,1) = RAir * xyz_Temp(ij,1) / Grav * ( 1.0d0 - z_Sigma(1) )
    end do

    do k = 2, km
       do ij = 1, im*jm
          xyz_GeoPot(ij,k) = xyz_GeoPot(ij,k-1) + RAir / Grav * xyr_Temp(ij,k) * r_DelSigma(k) / r_Sigma(k)
       end do
    end do

    !----- 半整数座標 -----
    xyr_GeoPot(:,:) = 0.0d0
    do k = 2, km+1
       do ij = 1, im*jm
          xyr_GeoPot(ij,k) = xyr_GeoPot(ij,k-1) + RAir / Grav * xyz_Temp(ij,k-1) * z_DelSigma(k-1) / z_Sigma(k-1)
       end do
    end do


!!!!! AFES での算出方法 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
!!!!!!    real(DBKIND)        :: APHI(km), BPHI(km), CPHI(km)
!    do k = 1, km
!       APHI(k) = ( r_Sigma(k)/z_Sigma(k) )**(RAir/Cp) - 1.D0
!       BPHI(k) = 1.D0 - ( r_Sigma(k+1)/z_Sigma(k) )**(RAir/Cp)
!       CPHI(k) = (r_Sigma(k)/z_Sigma(k))**(RAir/Cp) - (r_Sigma(k+1)/z_Sigma(k))**(RAir/Cp)
!    end do
!    BPHI(km) = 0.0d0
!
!    do ij = 1, im*jm
!       xyz_GeoPot(ij,1) = Cp/Grav*APHI(1)*xyz_Temp(ij,1) 
!       xyr_GeoPot(ij,1) = 0.0d0
!    end do
!
!    do k = 2, km
!       do ij = 1, im*jm
!       xyz_GeoPot(ij,k) = xyz_GeoPot(ij,k-1) &
!            &            + Cp/Grav * (APHI(k)*xyz_Temp(ij,k) &
!            &                       + BPHI(k-1)*xyz_Temp(ij,k-1))
!       end do
!    end do
!
!    do k = 2, km+1
!       do ij = 1, im*jm
!          xyr_GeoPot(ij,k) = xyr_GeoPot(ij,k-1) &
!               &            + Cp/Grav * CPHI(k-1) * xyz_Temp(ij,k-1)
!       end do
!    end do

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_interpolate_geopot

  subroutine physics_interpolate_temp(xyr_Temp, xyz_Temp, z_Sigma, r_Sigma)

    !==== Dependency
    use type_mod,    only : REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use dc_trace,    only : SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Input
    !
    real(DBKIND), intent(in) :: z_Sigma(:)        , r_Sigma(:)        , xyz_Temp(im*jm,km)    ! intent(in): 温度 (整数)座標

    !==== Output
    !
    real(DBKIND), intent(out) :: xyr_Temp(im*jm,km+1) ! intent(out): 温度 (半整数)座標

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_interpolate_temp"
    real(DBKIND)        :: z_TempInpolM(km), z_TempInpol(km)            ! 補間係数
    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: ij, k

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

    !----------------------------------------------------------------
    !   温度補間演算
    !----------------------------------------------------------------
    
    !----- 1. 温度補間係数 -----
    do k = 2, km
       z_TempInpolM(k) = (  log( r_Sigma( k )  ) - log( z_Sigma( k )  )  ) / (  log( z_Sigma( k-1 ) ) - log( z_Sigma( k )  )  )
       z_TempInpol(k)  = (  log( z_Sigma( k-1 ) ) - log( r_Sigma( k ) )  ) / (  log( z_Sigma( k-1 ) ) - log( z_Sigma( k )  )  )
    end  do

    !----- 2. 温度の補間 -----
    do k = 2, km
       do ij = 1, im*jm
          xyr_Temp(ij,k) = z_TempInpolM(k) * xyz_Temp(ij,k-1) + z_TempInpol(k) * xyz_Temp(ij,k) 
       end do
    end do
    
    do ij = 1, im*jm
       xyr_Temp(ij,km+1) = xyz_Temp(ij,km)
       xyr_Temp(ij,   1) = xyz_Temp(ij, 1)
    end do

    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_interpolate_temp