!= プランク関数の計算
!
!= Calculate Planck function
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: planck_func.f90,v 1.6 2014/05/07 09:39:21 murashin Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module planck_func
  !
  != プランク関数の計算
  !
  != Calculate Planck function
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! 
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
!!$  ! NAMELIST#radiation_DennouAGCM_nml
  !

  ! モジュール引用 ; USE statements
  !

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 


  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: aaa_PF
  public :: PF
  public :: DPFDT
  public :: Integ_PF_GQ_Array3D
  public :: Integ_PF_GQ_Array2D
  public :: Integ_DPFDT_GQ_Array2D
  public :: Integ_DPFDT_GQ_Array3D

  ! 公開変数
  ! Public variables
  !
!!$  logical, save :: planck_func_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag


  ! 非公開変数
  ! Private variables
  !
  real(DP), parameter ::                &
    & SOL    = 2.99792458e8_DP        , &
    & Planc  = 6.6260755e-34_DP       , &
    & Boltz  = 1.380658e-23_DP

  character(*), parameter:: module_name = 'planck_func'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: planck_func.f90,v 1.6 2014/05/07 09:39:21 murashin Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

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

  function aaa_PF( &
    & is, ie, js, je, ks, ke, &
    & WN, aaa_Temp &
    & ) &
    result( aaa_Res )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! Calculate radiation flux from temperature, specific humidity, and 
    ! air pressure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 宣言文 ; Declaration statements
    !
    integer , intent(in) :: is
    integer , intent(in) :: ie
    integer , intent(in) :: js
    integer , intent(in) :: je
    integer , intent(in) :: ks
    integer , intent(in) :: ke
    real(DP), intent(in) :: WN
    real(DP), intent(in) :: aaa_Temp(is:ie, js:je, ks:ke)
    real(DP)             :: aaa_Res (is:ie, js:je, ks:ke)

    ! 作業変数
    ! Work variables
    !

    ! 実行文 ; Executable statement
    !

    aaa_Res = 2.0_DP * Planc * SOL * SOL * WN * WN * WN &
      / ( exp( Planc * SOL * ( WN+1.0e-10_DP ) / ( Boltz * aaa_Temp ) ) - 1.0_DP )


  end function aaa_PF

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

  function PF( WN, Temp ) result( Res )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! Calculate radiation flux from temperature, specific humidity, and 
    ! air pressure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in) :: WN
    real(DP), intent(in) :: Temp
    real(DP)             :: Res

    ! 作業変数
    ! Work variables
    !
    real(DP) :: aaa_Temp(1,1,1)
    real(DP) :: aaa_Res (1,1,1)

    ! 実行文 ; Executable statement
    !

    aaa_Temp(1,1,1) = Temp
    aaa_Res = &
      & aaa_PF( &
      &         1, 1, 1, 1, 1, 1, &
      &         WN, aaa_Temp &
      &       )

    Res = aaa_Res(1,1,1)


  end function PF

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

  subroutine Integ_PF_GQ_Array3D( &
    & wn1, wn2, num, &
    & is, ie, js, je, ks, ke, &
    & aaa_Temp, &
    & aaa_PFInted &
    & )

    ! ガウス重み, 分点の計算
    ! Calculate Gauss node and Gaussian weight
    !
    use gauss_quad, only : GauLeg

    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: is, ie
    integer , intent(in ) :: js, je
    integer , intent(in ) :: ks, ke
    real(DP), intent(in ) :: aaa_Temp   (is:ie, js:je, ks:ke)
    real(DP), intent(out) :: aaa_PFInted(is:ie, js:je, ks:ke)


    !
    ! local variables
    !
    real(DP):: x( num ), w( num )
    integer :: l


    call GauLeg( wn1, wn2, num, x, w )

    aaa_PFInted(:,:,:) = 0.0_DP

    do l = 1, num
      aaa_PFInted(:,:,:) = aaa_PFInted(:,:,:)       &
        & + aaa_PF(                         &
        &           is, ie, js, je, ks, ke, &
        &           x(l), aaa_Temp          &
        &         )                         &
        & * w( l )
    end do


  end subroutine Integ_PF_GQ_Array3D

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

  subroutine Integ_PF_GQ_Array2D( &
    & wn1, wn2, num, &
    & is, ie, js, je, &
    & xy_Temp, &
    & xy_PFInted &
    & )


    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    integer , intent(in ) :: js
    integer , intent(in ) :: je
    real(DP), intent(in ) :: xy_Temp   (is:ie, js:je)
    real(DP), intent(out) :: xy_PFInted(is:ie, js:je)


    !
    ! local variables
    !
    real(DP) :: xyz_Temp3d   (is:ie, js:je, 1:1)
    real(DP) :: xyz_PFInted3d(is:ie, js:je, 1:1)


    xyz_Temp3d(:,:,1) = xy_Temp(:,:)
    call Integ_PF_GQ_Array3D( &
      & wn1, wn2, num, &
      & is, ie, js, je, 1, 1, &
      & xyz_Temp3d, &
      & xyz_PFInted3d &
      & )
    xy_PFInted(:,:) = xyz_PFInted3d(:,:,1)


  end subroutine Integ_PF_GQ_Array2D

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

  function DPFDT( &
    & WN,    & ! (in )
    & Temp   & ! (in )
    & )      &
    & result( Res )

    ! USE statements
    !

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in ) :: WN
    real(DP), intent(in ) :: Temp
    real(DP)              :: Res


    ! 作業変数
    ! Work variables
    !
    real(DP) :: aaa_Temp(1,1,1)
    real(DP) :: aaa_Res (1,1,1)


    aaa_Temp(1,1,1) = Temp

    aaa_Res = aaa_DPFDT(                  &
      &                 1, 1, 1, 1, 1, 1, & ! (in )
      &                 WN,               & ! (in )
      &                 aaa_Temp          & ! (in )
      & )

    Res = aaa_Res(1,1,1)


  end function DPFDT

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

  function aaa_DPFDT( &
    & is, ie, js, je, ks, ke, & ! (in )
    & WN,                     & ! (in )
    & aaa_Temp                & ! (in )
    & )                       &
    & result( aaa_Res )

    ! USE statements
    !

    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    integer , intent(in ) :: js
    integer , intent(in ) :: je
    integer , intent(in ) :: ks
    integer , intent(in ) :: ke
    real(DP), intent(in ) :: WN
    real(DP), intent(in ) :: aaa_Temp(is:ie, js:je, ks:ke)
    real(DP)              :: aaa_Res (is:ie, js:je, ks:ke)


    real(DP) :: aaa_ExpTerm(is:ie, js:je, ks:ke)
    real(DP) :: aaa_PF     (is:ie, js:je, ks:ke)


    aaa_ExpTerm = exp( Planc * SOL * ( WN + 1.0e-10_DP ) / ( Boltz * aaa_Temp ) )

    aaa_PF = 2.0_DP * Planc * SOL * SOL * WN * WN * WN &
      / ( aaa_ExpTerm - 1.0_DP )

    aaa_Res = &
      & 1.0_DP / ( 2.0_DP * SOL * WN * WN * Boltz ) &
      & * ( aaa_PF / aaa_Temp )**2 &
      & * aaa_ExpTerm


  end function aaa_DPFDT

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

  subroutine Integ_DPFDT_GQ_Array3D(    &
    & WN1, WN2, Num,                    & ! (in )
    & is, ie, js, je, ks, ke, aaa_Temp, & ! (in )
    & aaa_DPFDTInted                    & ! (out)
    & )

    ! USE statements
    !

    ! ガウス重み, 分点の計算
    ! Calculate Gauss node and Gaussian weight
    !
    use gauss_quad, only : GauLeg

    real(DP), intent(in ) :: WN1
    real(DP), intent(in ) :: WN2
    integer , intent(in ) :: Num
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    integer , intent(in ) :: js
    integer , intent(in ) :: je
    integer , intent(in ) :: ks
    integer , intent(in ) :: ke
    real(DP), intent(in ) :: aaa_Temp      (is:ie, js:je, ks:ke)
    real(DP), intent(out) :: aaa_DPFDTInted(is:ie, js:je, ks:ke)


    !
    ! local variables
    !
    real(DP):: GP( Num )
    real(DP):: GW( Num )
    integer :: l


    call GauLeg( WN1, WN2, Num, GP, GW )

    aaa_DPFDTInted = 0.0_DP

    do l = 1, num
      aaa_DPFDTInted = aaa_DPFDTInted &
        & + aaa_DPFDT(                         &
        &              is, ie, js, je, ks, ke, & ! (in )
        &              GP(l),                  & ! (in )
        &              aaa_Temp                & ! (in )
        &           ) &
        & * GW(l)
    end do


  end subroutine Integ_DPFDT_GQ_Array3D

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

  subroutine Integ_DPFDT_GQ_Array2D( &
    & WN1, WN2, Num,                 & ! (in )
    & is, ie, js, je, aa_Temp,       & ! (in )
    & aa_DPFDTInted                  & ! (out)
    & )

    ! USE statements
    !

    real(DP), intent(in ) :: WN1
    real(DP), intent(in ) :: WN2
    integer , intent(in ) :: Num
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    integer , intent(in ) :: js
    integer , intent(in ) :: je
    real(DP), intent(in ) :: aa_Temp      (is:ie, js:je)
    real(DP), intent(out) :: aa_DPFDTInted(is:ie, js:je)


    !
    ! local variables
    !
    real(DP) :: aaa_Temp      (is:ie, js:je, 1:1)
    real(DP) :: aaa_DPFDTInted(is:ie, js:je, 1:1)


    aaa_Temp(:,:,1) = aa_Temp

    call Integ_DPFDT_GQ_Array3D(        &
      & WN1, WN2, Num,                  & ! (in )
      & is, ie, js, je, 1, 1, aaa_Temp, & ! (in )
      & aaa_DPFDTInted                  & ! (out)
      & )

    aa_DPFDTInted = aaa_DPFDTInted(:,:,1)


  end subroutine Integ_DPFDT_GQ_Array2D

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

  subroutine PlanckFuncPrepPFTable(                            &
    & WNs, WNe, NGaussQuad, TableTempMin, TableTempMax, ntmax, &
    & a_TableTemp, a_TableIPF, a_TableIDPFDT                   &
    & )

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    ! ガウス重み, 分点の計算
    ! Calculate Gauss node and Gaussian weight
    !
    use gauss_quad, only : GauLeg


    real(DP), intent(in ) :: WNs
    real(DP), intent(in ) :: WNe
    integer , intent(in ) :: NGaussQuad
    real(DP), intent(in ) :: TableTempMin
    real(DP), intent(in ) :: TableTempMax
    integer , intent(in ) :: ntmax
    real(DP), intent(out) :: a_TableTemp  (1:ntmax)
    real(DP), intent(out) :: a_TableIPF   (1:ntmax)
    real(DP), intent(out) :: a_TableIDPFDT(1:ntmax)


    ! Local variables
    !
    real(DP)              :: TableTempIncrement
    integer               :: nn
    real(DP), allocatable :: aa_TempTMP   (:,:)
    real(DP), allocatable :: aa_PF        (:,:)
    real(DP), allocatable :: aa_DPFDT     (:,:)
    real(DP), allocatable :: aa_PFTable   (:,:)
    real(DP), allocatable :: aa_DPFDTTable(:,:)
    real(DP)              :: ErrorPFInteg
    real(DP), parameter   :: ThresholdErrorPFInteg = 1.0d-3
                              ! Threshold for checking accuracy of calculation of
                              ! integrated Planc function by using a pre-calculated
                              ! table.

    ! Variables for preparation for calculation of Plank function
    !
    real(DP)      , allocatable :: a_GQP(:)
    real(DP)      , allocatable :: a_GQW(:)


    integer:: i
    integer:: j
    integer:: l
    integer:: m


    ! Preparation of tables for calculation of Plank function
    !
    TableTempIncrement = ( TableTempMax - TableTempMin ) / ( ntmax - 1 )

    do m = 1, ntmax
      a_TableTemp(m) = TableTempMin + TableTempIncrement * ( m - 1 )
    end do


    a_TableIPF    = 0.0_DP
    a_TableIDPFDT = 0.0_DP

    allocate( a_GQP(1:NGaussQuad) )
    allocate( a_GQW(1:NGaussQuad) )
    call GauLeg(              &
      & WNs, WNe, NGaussQuad, & ! (in )
      & a_GQP, a_GQW          & ! (out)
      & )
    do m = 1, ntmax
      do l = 1, NGaussQuad
        a_TableIPF   (m) = a_TableIPF   (m) + PF   ( a_GQP(l), a_TableTemp(m) ) * a_GQW(l)
        a_TableIDPFDT(m) = a_TableIDPFDT(m) + DPFDT( a_GQP(l), a_TableTemp(m) ) * a_GQW(l)
      end do
    end do
    deallocate( a_GQP )
    deallocate( a_GQW )


    !----------------------------------------------------
    ! Check accuracy of integration of Planc function by using a pre-calculated table.
    !

    nn = ntmax-1
    allocate( aa_TempTMP   (1:nn, 1:1) )
    allocate( aa_PF        (1:nn, 1:1) )
    allocate( aa_DPFDT     (1:nn, 1:1) )
    allocate( aa_PFTable   (1:nn, 1:1) )
    allocate( aa_DPFDTTable(1:nn, 1:1) )

    j = 1

    do i = 1, nn
      aa_TempTMP(i,j) =                                                &
        &   a_TableTemp(1)                                             &
        & + ( a_TableTemp(2) - a_TableTemp(1) ) * 0.5_DP               &
        & + ( a_TableTemp(2) - a_TableTemp(1) ) * ( i - 1 )
    end do



    call Integ_PF_GQ_Array2D(    &
      & WNs, WNe, NGaussQuad,    &
      & 1, nn, 1, 1, aa_TempTMP, &
      & aa_PF                    &
      & )
    call Integ_DPFDT_GQ_Array2D(     &
      & WNs, WNe, NGaussQuad,        & ! (in )
      & 1, nn, 1, 1, aa_TempTMP,     & ! (in )
      & aa_DPFDT                     & ! (out)
      & )

    call CalcIntegratedPFWithTable2D(   &
      & ntmax, a_TableTemp, a_TableIPF, &
      & 1, nn, 1, 1, aa_TempTMP,        &
      & aa_PFTable,                     &
      & .false.                         &
      & )
    call CalcIntegratedPFWithTable2D(      &
      & ntmax, a_TableTemp, a_TableIDPFDT, &
      & 1, nn, 1, 1, aa_TempTMP,           &
      & aa_DPFDTTable,                     &
      & .true.                             &
      & )

    do i = 1, nn
      ErrorPFInteg = abs( aa_PF   (i,j) - aa_PFTable   (i,j) ) / aa_PF   (i,j)
      if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
        call MessageNotify( 'E', module_name, &
          & 'Error of integrated PF, %f, is greater than threshold, %f.', &
          & d = (/ ErrorPFInteg, ThresholdErrorPFInteg /) )
      end if
      ErrorPFInteg = abs( aa_DPFDT(i,j) - aa_DPFDTTable(i,j) ) / aa_DPFDT(i,j)
      if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
        call MessageNotify( 'E', module_name, &
          & 'Error of integrated DPFDT, %f, is greater than threshold, %f.', &
          & d = (/ ErrorPFInteg, ThresholdErrorPFInteg /) )
      end if
    end do

    deallocate( aa_TempTMP    )
    deallocate( aa_PF         )
    deallocate( aa_DPFDT      )
    deallocate( aa_PFTable    )
    deallocate( aa_DPFDTTable )


  end subroutine PlanckFuncPrepPFTable

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

  subroutine CalcIntegratedPFWithTable2D( &
    & ntmax, a_TableTemp, a_TableIPF,     &
    & is, ie, js, je, xy_Temp,            &
    & xy_IntegPF,                         &
    & flag_DPFDT                          &
    & )

    ! USE statements
    !

    integer , intent(in )           :: ntmax
    real(DP), intent(in )           :: a_TableTemp(1:ntmax)
    real(DP), intent(in )           :: a_TableIPF (1:ntmax)
    integer , intent(in )           :: is
    integer , intent(in )           :: ie
    integer , intent(in )           :: js
    integer , intent(in )           :: je
    real(DP), intent(in )           :: xy_Temp    (is:ie, js:je)
    real(DP), intent(out)           :: xy_IntegPF (is:ie, js:je)
    logical , intent(in ), optional :: flag_DPFDT

    !
    ! local variables
    !
    real(DP) :: xyz_Temp   (is:ie, js:je, 1)
    real(DP) :: xyz_IntegPF(is:ie, js:je, 1)


    xyz_Temp(:,:,1) = xy_Temp

    call CalcIntegratedPFWithTable3D( &
      & ntmax, a_TableTemp, a_TableIPF,     &
      & is, ie, js, je, 1, 1, xyz_Temp,     &
      & xyz_IntegPF,                        &
      & flag_DPFDT                          &
      & )

    xy_IntegPF = xyz_IntegPF(:,:,1)


  end subroutine CalcIntegratedPFWithTable2D

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

  subroutine CalcIntegratedPFWithTable3D( &
    & ntmax, a_TableTemp, a_TableIPF,     &
    & is, ie, js, je, ks, ke, xyz_Temp,   &
    & xyz_IntegPF,                        &
    & flag_DPFDT                          &
    & )

    ! USE statements
    !

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    integer , intent(in )           :: ntmax
    real(DP), intent(in )           :: a_TableTemp(1:ntmax)
    real(DP), intent(in )           :: a_TableIPF (1:ntmax)
    integer , intent(in )           :: is
    integer , intent(in )           :: ie
    integer , intent(in )           :: js
    integer , intent(in )           :: je
    integer , intent(in )           :: ks
    integer , intent(in )           :: ke
    real(DP), intent(in )           :: xyz_Temp   (is:ie, js:je, ks:ke)
    real(DP), intent(out)           :: xyz_IntegPF(is:ie, js:je, ks:ke)
    logical , intent(in ), optional :: flag_DPFDT

    !
    ! local variables
    !
    real(DP)                    :: TableTempMin
    real(DP)                    :: TableTempMax
    real(DP)                    :: TableTempIncrement

    logical                     :: local_flag_DPFDT

    integer                     :: xyz_TempIndex(is:ie, js:je, ks:ke)
    integer                     :: i
    integer                     :: j
    integer                     :: k
    integer                     :: m


    TableTempMin       = a_TableTemp(1)
    TableTempMax       = a_TableTemp(ntmax)
    TableTempIncrement = ( TableTempMax - TableTempMin ) / ( ntmax - 1 )


    do k = ks, ke
      do j = js, je
        do i = is, ie

          if ( ( xyz_Temp(i,j,k) < a_TableTemp(1)     ) .or. &
            &  ( xyz_Temp(i,j,k) > a_TableTemp(ntmax) ) ) then
            call MessageNotify( 'E', module_name, &
              & 'Temperature is not appropriate, Temp(%d,%d,%d) = %f.', &
              & i = (/i, j, k/), d = (/xyz_Temp(i,j,k)/) )
          end if

          xyz_TempIndex(i,j,k) = &
            & int( ( xyz_Temp(i,j,k) - TableTempMin ) / TableTempIncrement ) + 2

          if ( xyz_TempIndex(i,j,k) == 1 ) then
             xyz_TempIndex(i,j,k) = 2
          else if ( xyz_TempIndex(i,j,k) >= ntmax ) then
             xyz_TempIndex(i,j,k) = ntmax - 1
          end if

!!$          xyz_TempIndex(i,j,k) = ntmax-1
!!$          search_index: do m = 2, ntmax-1
!!$            if ( a_TableTemp(m) >= xyz_Temp(i,j,k) ) then
!!$              xyz_TempIndex(i,j,k) = m
!!$              exit search_index
!!$            end if
!!$          end do search_index

        end do
      end do
    end do


    local_flag_DPFDT = .false.
    if ( present( flag_DPFDT ) ) then
      if ( flag_DPFDT ) then
        local_flag_DPFDT = .true.
      end if
    end if

    if ( .not. local_flag_DPFDT ) then
      do k = ks, ke
        do j = js, je
          do i = is, ie
            m = xyz_TempIndex(i,j,k)

!!$            xyz_IntegPF(i,j,k) = &
!!$              &   ( aa_TableIPF( m, iband ) - aa_TableIPF( m-1, iband ) ) &
!!$              & / ( a_TableTemp( m )        - a_TableTemp( m-1 )        ) &
!!$              & * ( xyz_Temp(i,j,k)         - a_TableTemp( m-1 )        ) &
!!$              & +   aa_TableIPF( m-1, iband )

            xyz_IntegPF(i,j,k) = &
              &   a_TableIPF(m-1)                                  &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) )    &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) )    &
              & / ( ( a_TableTemp( m-1 ) - a_TableTemp( m   ) )    &
              &   * ( a_TableTemp( m-1 ) - a_TableTemp( m+1 ) ) )  &
              & + a_TableIPF(m  )                                  &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) )    &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) )    &
              & / ( ( a_TableTemp( m   ) - a_TableTemp( m-1 ) )    &
              &   * ( a_TableTemp( m   ) - a_TableTemp( m+1 ) ) )  &
              & + a_TableIPF(m+1)                                  &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) )    &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) )    &
              & / ( ( a_TableTemp( m+1 ) - a_TableTemp( m-1 ) )    &
              &   * ( a_TableTemp( m+1 ) - a_TableTemp( m   ) ) )
          end do
        end do
      end do
    else
      do k = ks, ke
        do j = js, je
          do i = is, ie
            m = xyz_TempIndex(i,j,k)

!!$            xyz_IntegPF(i,j,k) = &
!!$              &   ( aa_TableIDPFDT( m, iband ) - aa_TableIDPFDT( m-1, iband ) ) &
!!$              & / ( a_TableTemp   ( m )        - a_TableTemp   ( m-1 )        ) &
!!$              & * ( xyz_Temp(i,j,k)         - a_TableTemp( m-1 )        ) &
!!$              & +   aa_TableIPF( m-1, iband )

            xyz_IntegPF(i,j,k) = &
              &   a_TableIPF(m-1)                                  &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) )    &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) )    &
              & / ( ( a_TableTemp( m-1 ) - a_TableTemp( m   ) )    &
              &   * ( a_TableTemp( m-1 ) - a_TableTemp( m+1 ) ) )  &
              & + a_TableIPF(m  )                                  &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) )    &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) )    &
              & / ( ( a_TableTemp( m   ) - a_TableTemp( m-1 ) )    &
              &   * ( a_TableTemp( m   ) - a_TableTemp( m+1 ) ) )  &
              & + a_TableIPF(m+1)                                  &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) )    &
              &   * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) )    &
              & / ( ( a_TableTemp( m+1 ) - a_TableTemp( m-1 ) )    &
              &   * ( a_TableTemp( m+1 ) - a_TableTemp( m   ) ) )
          end do
        end do
      end do
    end if


  end subroutine CalcIntegratedPFWithTable3D

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

end module planck_func




