プランク関数の計算

Calculate Planck function

Note that Japanese and English are described in parallel.

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

Methods

Boltz CalcIntegratedPFWithTable2D CalcIntegratedPFWithTable3D DPFDT Integ_DPFDT_GQ_Array2D Integ_DPFDT_GQ_Array3D Integ_PF_GQ_Array2D Integ_PF_GQ_Array3D PF Planc PlanckFuncPrepPFTable SOL aaa_DPFDT aaa_PF module_name version

Included Modules

dc_types gauss_quad dc_message

Public Instance methods

DPFDT( WN, Temp ) result(Res)

*Function* :
`Res` :	`real(DP)`
`WN` :	`real(DP), intent(in )`
`Temp` :	`real(DP), intent(in )`

[Source]

  function DPFDT( WN, Temp ) 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, WN, aaa_Temp )

    Res = aaa_Res(1,1,1)


  end function DPFDT

Integ_DPFDT_GQ_Array2D( WN1, WN2, Num, is, ie, js, je, aa_Temp, aa_DPFDTInted )

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

[Source]

  subroutine Integ_DPFDT_GQ_Array2D( WN1, WN2, Num, is, ie, js, je, aa_Temp, aa_DPFDTInted )

    ! 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, is, ie, js, je, 1, 1, aaa_Temp, aaa_DPFDTInted )

    aa_DPFDTInted = aaa_DPFDTInted(:,:,1)


  end subroutine Integ_DPFDT_GQ_Array2D

Integ_DPFDT_GQ_Array3D( WN1, WN2, Num, is, ie, js, je, ks, ke, aaa_Temp, aaa_DPFDTInted )

*Subroutine* :
`WN1` :	`real(DP), intent(in )`
`WN2` :	`real(DP), 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` :	`integer , intent(in )`
`aaa_Temp(is:ie, js:je, ks:ke)` :	`real(DP), intent(in )`
`aaa_DPFDTInted(is:ie, js:je, ks:ke)` :	`real(DP), intent(out)`

[Source]

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

    ! 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, GP(l), aaa_Temp ) * GW(l)
    end do


  end subroutine Integ_DPFDT_GQ_Array3D

Integ_PF_GQ_Array2D( wn1, wn2, num, is, ie, js, je, temp, pfinted )

*Subroutine* :
`wn1` :	`real(DP), intent(in )`
`wn2` :	`real(DP), intent(in )`
`num` :	`integer , intent(in )`
`is` :	`integer , intent(in )`
`ie` :	`integer , intent(in )`
`js` :	`integer , intent(in )`
`je` :	`integer , intent(in )`
`temp(is:ie, js:je)` :	`real(DP), intent(in )`
`pfinted(is:ie, js:je)` :	`real(DP), intent(out)`

[Source]

  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

Integ_PF_GQ_Array3D( wn1, wn2, num, is, ie, js, je, ks, ke, aaa_temp, aaa_pfinted )

*Subroutine* :
`wn1` :	`real(DP), intent(in )`
`wn2` :	`real(DP), 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` :	`integer , intent(in )`
`aaa_temp(is:ie, js:je, ks:ke)` :	`real(DP), intent(in )`
`aaa_pfinted(is:ie, js:je, ks:ke)` :	`real(DP), intent(out)`

[Source]

  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

PF( WN, Temp ) result(Res)

*Function* :
`Res` :	`real(DP)`
`WN` :	`real(DP), intent(in)`
`Temp` :	`real(DP), intent(in)`

温度, 比湿, 気圧から, 放射フラックスを計算します.

Calculate radiation flux from temperature, specific humidity, and air pressure.

[Source]

  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

aaa_PF( is, ie, js, je, ks, ke, WN, aaa_Temp ) result(aaa_Res)

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

温度, 比湿, 気圧から, 放射フラックスを計算します.

Calculate radiation flux from temperature, specific humidity, and air pressure.

[Source]

  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

Private Instance methods

Boltz

*Constant* :
`Boltz` `= 1.380658d-23` :	`real(DP), parameter`

CalcIntegratedPFWithTable2D( ntmax, a_TableTemp, a_TableIPF, is, ie, js, je, xy_Temp, xy_IntegPF, [flag_DPFDT] )

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

[Source]

  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

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

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

[Source]

  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

Planc

*Constant* :
`Planc` `= 6.6260755d-34` :	`real(DP), parameter`

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

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

[Source]

  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, a_GQP, a_GQW )
    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, 1, nn, 1, 1, aa_TempTMP, aa_DPFDT )

    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

SOL

*Constant* :
`SOL` `= 2.99792458d8` :	`real(DP), parameter`

aaa_DPFDT( is, ie, js, je, ks, ke, WN, aaa_Temp ) result(aaa_Res)

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

[Source]

  function aaa_DPFDT( is, ie, js, je, ks, ke, WN, aaa_Temp ) 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

module_name

Constant :

module_name = ‘planck_func‘ :

character(*), parameter

:	モジュールの名称. Module name

version

Constant :

version = ’$Name: dcpam5-20121128 $’ // ’$Id: planck_func.f90,v 1.5 2012-01-20 00:24:54 yot Exp $’ :

character(*), parameter

:	モジュールのバージョン Module version