! Calculate Planck function

module planck_func_v2

  use vtype_module

  implicit none

  private

  public :: PFV2IntegPFGQArray3D
  public :: PFV2IntegPFGQArray2D
  public :: PFV2IntegPFGQArray1D
  public :: PFV2IntegDPFDTGQArray3D
  public :: PFV2IntegDPFDTGQArray2D
  public :: PFV2IntegDPFDTGQArray1D


  real(DP), parameter ::                &
    & SOL    = 2.99792458e8_DP        , &
    & Planc  = 6.6260755e-34_DP       , &
    & Boltz  = 1.380658e-23_DP

  character(*), parameter:: module_name = 'planck_func'

contains

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

  subroutine CalcPF3D( &
    & imax, jmax, kmax,       &
    & is, ie, js, je, ks, ke, &
    & WN, aaa_Temp,           &
    & aaa_Res                 &
    & )

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

    ! 宣言文 ; Declaration statements
    !
    integer , intent(in ) :: imax
    integer , intent(in ) :: jmax
    integer , intent(in ) :: kmax
    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(imax, jmax, kmax)
    real(DP), intent(out) :: aaa_Res (imax, jmax, kmax)


    ! 作業変数
    ! Work variables
    !
    integer :: i
    integer :: j
    integer :: k


    ! 実行文 ; Executable statement
    !

    ! NOTE:
    !   When wavenumber is large and temperature is low, an argument of
    !   exponential function is very large. This may cause errors.

    do k = ks, ke
      do j = js, je
        do i = is, ie
          aaa_Res(i,j,k) = 2.0_DP * Planc * SOL * SOL * WN * WN * WN &
            / ( exp( Planc * SOL * ( WN+1.0e-10_DP ) / ( Boltz * aaa_Temp(i,j,k) ) ) - 1.0_DP )
        end do
      end do
    end do


  end subroutine CalcPF3D

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

!!$  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
!!$
!!$    call CalcPF3D( &
!!$      &         1, 1, 1,          &
!!$      &         1, 1, 1, 1, 1, 1, &
!!$      &         WN, aaa_Temp,     &
!!$      &         aaa_Res           &
!!$      &          )
!!$
!!$    Res = aaa_Res(1,1,1)
!!$
!!$
!!$  end function PF

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

  subroutine PFV2IntegPFGQArray3D( &
    & wn1, wn2, num,          &
    & imax, jmax, kmax,       &
    & 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 ) :: imax
    integer , intent(in ) :: jmax
    integer , intent(in ) :: kmax
    integer , intent(in ) :: is, ie
    integer , intent(in ) :: js, je
    integer , intent(in ) :: ks, ke
    real(DP), intent(in ) :: aaa_temp   (imax, jmax, kmax)
    real(DP), intent(out) :: aaa_pfinted(imax, jmax, kmax)


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

    real(DP):: aaa_PFEach(imax, jmax, kmax)

    integer :: i
    integer :: j
    integer :: k
    integer :: l


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

    do k = ks, ke
      do j = js, je
        do i = is, ie
          aaa_pfinted(i,j,k) = 0.0_DP
        end do
      end do
    end do

    do l = 1, num

      call CalcPF3D( &
        & imax, jmax, kmax,       &
        & is, ie, js, je, ks, ke, &
        & x(l), aaa_Temp,         &
        & aaa_PFEach              &
        & )
      do k = ks, ke
        do j = js, je
          do i = is, ie
            aaa_pfinted(i,j,k) = aaa_pfinted(i,j,k)       &
              & + aaa_PFEach(i,j,k) * w( l )
          end do
        end do
      end do

    end do


  end subroutine PFV2IntegPFGQArray3D

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

  subroutine PFV2IntegPFGQArray2D( &
    & wn1, wn2, num, &
    & imax, jmax, &
    & is, ie, js, je, &
    & aa_temp, &
    & aa_pfinted &
    & )


    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: imax
    integer , intent(in ) :: jmax
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    integer , intent(in ) :: js
    integer , intent(in ) :: je
    real(DP), intent(in ) :: aa_temp   (imax, jmax)
    real(DP), intent(out) :: aa_pfinted(imax, jmax)


    !
    ! local variables
    !
    real(DP) :: aaa_temp   (imax, jmax, 1:1)
    real(DP) :: aaa_pfinted(imax, jmax, 1:1)

    integer  :: i
    integer  :: j


    do j = js, je
      do i = is, ie
        aaa_temp(i,j,1) = aa_temp(i,j)
      end do
    end do

    call PFV2IntegPFGQArray3D( &
      & wn1, wn2, num,        &
      & imax, jmax, 1,        &
      & is, ie, js, je, 1, 1, &
      & aaa_temp,             &
      & aaa_pfinted           &
      & )

    do j = js, je
      do i = is, ie
        aa_pfinted(i,j) = aaa_pfinted(i,j,1)
      end do
    end do


  end subroutine PFV2IntegPFGQArray2D

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

  subroutine PFV2IntegPFGQArray1D( &
    & wn1, wn2, num, &
    & imax, &
    & is, ie, &
    & a_temp, &
    & a_pfinted &
    & )


    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: imax
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    real(DP), intent(in ) :: a_temp   (imax)
    real(DP), intent(out) :: a_pfinted(imax)


    !
    ! local variables
    !
    real(DP) :: aaa_temp   (imax, 1:1, 1:1)
    real(DP) :: aaa_pfinted(imax, 1:1, 1:1)

    integer  :: i


    do i = is, ie
      aaa_temp(i,1,1) = a_temp(i)
    end do

    call PFV2IntegPFGQArray3D( &
      & wn1, wn2, num,        &
      & imax, 1, 1,        &
      & is, ie, 1, 1, 1, 1, &
      & aaa_temp,             &
      & aaa_pfinted           &
      & )

    do i = is, ie
      a_pfinted(i) = aaa_pfinted(i,1,1)
    end do


  end subroutine PFV2IntegPFGQArray1D

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

!!$  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

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

  subroutine CalcDPFDT3D( &
    & imax, jmax, kmax,       &
    & is, ie, js, je, ks, ke, & ! (in )
    & WN,                     & ! (in )
    & aaa_Temp,               & ! (in )
    & aaa_Res                 &
    & )

    ! USE statements
    !

    integer , intent(in ) :: imax
    integer , intent(in ) :: jmax
    integer , intent(in ) :: kmax
    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(imax, jmax, kmax)
    real(DP), intent(out) :: aaa_Res (imax, jmax, kmax)


    real(DP) :: aaa_ExpTerm(imax, jmax, kmax)
    real(DP) :: aaa_PF     (imax, jmax, kmax)
    integer  :: i
    integer  :: j
    integer  :: k


    do k = ks, ke
      do j = js, je
        do i = is, ie
          aaa_ExpTerm(i,j,k) = exp( Planc * SOL * ( WN + 1.0e-10_DP ) / ( Boltz * aaa_Temp(i,j,k) ) )

          aaa_PF(i,j,k) = 2.0_DP * Planc * SOL * SOL * WN * WN * WN &
            / ( aaa_ExpTerm(i,j,k) - 1.0_DP )

          aaa_Res(i,j,k) = &
            & 1.0_DP / ( 2.0_DP * SOL * WN * WN * Boltz ) &
            & * ( aaa_PF(i,j,k) / aaa_Temp(i,j,k) )**2 &
            & * aaa_ExpTerm(i,j,k)
        end do
      end do
    end do


  end subroutine CalcDPFDT3D

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

  subroutine PFV2IntegDPFDTGQArray3D(    &
    & WN1, WN2, Num,                    & ! (in )
    & imax, jmax, kmax,                 &
    & 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 ) :: imax
    integer , intent(in ) :: jmax
    integer , intent(in ) :: kmax
    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      (imax, jmax, kmax)
    real(DP), intent(out) :: aaa_DPFDTInted(imax, jmax, kmax)


    !
    ! local variables
    !
    real(DP):: GP( Num )
    real(DP):: GW( Num )
    real(DP):: aaa_DPFDTEach(imax, jmax, kmax)

    integer :: i
    integer :: j
    integer :: k
    integer :: l


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

    do k = ks, ke
      do j = js, je
        do i = is, ie
          aaa_DPFDTInted(i,j,k) = 0.0_DP
        end do
      end do
    end do

    do l = 1, num
      call CalcDPFDT3D( &
        &              imax, jmax, kmax,       &
        &              is, ie, js, je, ks, ke, & ! (in )
        &              GP(l),                  & ! (in )
        &              aaa_Temp,               & ! (in )
        &              aaa_DPFDTEach           &
        &             )
      do k = ks, ke
        do j = js, je
          do i = is, ie
            aaa_DPFDTInted(i,j,k) = aaa_DPFDTInted(i,j,k) &
              & + aaa_DPFDTEach(i,j,k) * GW(l)
          end do
        end do
      end do
    end do


  end subroutine PFV2IntegDPFDTGQArray3D

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

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


    !
    ! local variables
    !
    real(DP) :: aaa_Temp      (imax, jmax, 1:1)
    real(DP) :: aaa_DPFDTInted(imax, jmax, 1:1)

    integer  :: i
    integer  :: j


    do j = js, je
      do i = is, ie
        aaa_Temp(i,j,1) = aa_Temp(i,j)
      end do
    end do

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

    do j = js, je
      do i = is, ie
        aa_DPFDTInted(i,j) = aaa_DPFDTInted(i,j,1)
      end do
    end do


  end subroutine PFV2IntegDPFDTGQArray2D

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

  subroutine PFV2IntegDPFDTGQArray1D( &
    & WN1, WN2, Num,                 & ! (in )
    & imax,                          &
    & is, ie, a_Temp,                & ! (in )
    & a_DPFDTInted                   & ! (out)
    & )

    ! USE statements
    !

    real(DP), intent(in ) :: WN1
    real(DP), intent(in ) :: WN2
    integer , intent(in ) :: Num
    integer , intent(in ) :: imax
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    real(DP), intent(in ) :: a_Temp      (imax)
    real(DP), intent(out) :: a_DPFDTInted(imax)


    !
    ! local variables
    !
    real(DP) :: aaa_Temp      (imax, 1:1, 1:1)
    real(DP) :: aaa_DPFDTInted(imax, 1:1, 1:1)

    integer  :: i


    do i = is, ie
      aaa_Temp(i,1,1) = a_Temp(i)
    end do

    call PFV2IntegDPFDTGQArray3D(        &
      & WN1, WN2, Num,                  & ! (in )
      & imax, 1, 1,                     &
      & is, ie, 1, 1, 1, 1, aaa_Temp,   & ! (in )
      & aaa_DPFDTInted                  & ! (out)
      & )

    do i = is, ie
      a_DPFDTInted(i) = aaa_DPFDTInted(i,1,1)
    end do


  end subroutine PFV2IntegDPFDTGQArray1D

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

end module planck_func_v2
