!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_interpolate.f90 
!
! History
!   2005/09/18 Yamada Yukiko     create
!

module physics_interpolate_mod

  use type_mod,    only : REKIND, DBKIND, INTKIND, TOKEN, STRING

  implicit none

  private
  public :: physics_interpolate_temp
  public :: physics_interpolate_geopot

contains

  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(:)        , & ! intent(in): σレベル(整数)座標
         & r_Sigma(:)        , & ! intent(in): σレベル(半整数)座標
         & 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

  subroutine physics_interpolate_geopot( &
       &  xyz_Press, xyr_Press, xyz_GeoPot, xyr_GeoPot, &   ! (out)
       &  xy_Press, xyz_Temp, xyr_Temp,                 &   ! (in)
       &  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),      & ! intent(in): 地表面気圧
         & xyz_Temp(im*jm,km),   & ! intent(in): 温度 (整数)座標
         & xyr_Temp(im*jm,km+1), & ! intent(in): 温度 (半整数)座標
         & z_Sigma(:),           & ! intent(in): σレベル(整数)座標
         & r_Sigma(:)              ! intent(in): σレベル(半整数)座標

    !==== Output
    !
    real(DBKIND), intent(out) ::  &
         & xyz_Press(im*jm,km),   & ! intent(out): 気圧 (整数)座標
         & xyr_Press(im*jm,km+1), & ! intent(out): 気圧 (半整数)座標
         & xyz_GeoPot(im*jm,km),  & ! intent(out): ジオポテンシャル(整数)座標
         & 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


end module physics_interpolate_mod