!= プランク関数の計算
!
!= Calculate Planck function
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: planck_func.f90,v 1.2 2010-09-18 01:38:22 yot Exp $ 
! Tag Name::  $Name: dcpam5-20110225-4 $
! 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, public:: planck_func_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag


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

  character(*), parameter:: module_name = 'planck_func'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20110225-4 $' // &
    & '$Id: planck_func.f90,v 1.2 2010-09-18 01:38:22 yot 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.0d0 * Planc * SOL * SOL * WN * WN * WN &
      / ( exp( Planc * SOL * ( WN+1.0d-10 ) / ( Boltz * aaa_Temp ) ) - 1.0d0 )


  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.0d0

    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, &
    & temp, &
    & 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 ) :: temp   (is:ie, js:je)
    real(DP), intent(out) :: pfinted(is:ie, js:je)


    !
    ! local variables
    !
    real(DP) :: temp3d   (is:ie, js:je, 1:1)
    real(DP) :: pfinted3d(is:ie, js:je, 1:1)


    temp3d(:,:,1) = temp(:,:)
    call Integ_PF_GQ_Array3D( &
      & wn1, wn2, num, &
      & is, ie, js, je, 1, 1, &
      & temp3d, &
      & pfinted3d &
      & )
    pfinted(:,:) = 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.0d-10 ) / ( Boltz * aaa_Temp ) )

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

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


  end function aaa_DPFDT

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

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

    ! USE statements
    !

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

    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 ) :: WN1
    real(DP), intent(in ) :: WN2
    integer , intent(in ) :: Num
    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( &
    & is, ie, js, je,         & ! (in )
    & WN1, WN2, Num, aa_Temp, & ! (in )
    & aa_DPFDTInted           & ! (out)
    & )

    ! USE statements
    !

    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    integer , intent(in ) :: js
    integer , intent(in ) :: je
    real(DP), intent(in ) :: WN1
    real(DP), intent(in ) :: WN2
    integer , intent(in ) :: Num
    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( &
      & is, ie, js, je, 1, 1,    & ! (in )
      & WN1, WN2, Num, aaa_Temp, & ! (in )
      & aaa_DPFDTInted           & ! (out)
      & )

    aa_DPFDTInted = aaa_DPFDTInted(:,:,1)


  end subroutine Integ_DPFDT_GQ_Array2D

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

  end module planck_func
