!= 地球大気向け長波放射モデル Ver. 2.4
!
!= long wave radiation model for the Earth's atmosphere Ver. 2.4
!
! Authors::   Yoshiyuki O. TAKAHASHI
! Version::   $Id: rad_Earth_LW_V2_4.f90,v 1.3 2015/01/29 12:06:43 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module rad_Earth_LW_V2_4
  !
  != 地球大気向け長波放射モデル Ver. 2.4
  !
  != long wave radiation model for the Earth's atmosphere Ver. 2.4
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 長波放射モデル.
  !
  ! This is a model of long wave radiation for the Earth's atmospehre. 
  ! Radiation in the wavenumber range from    0 to  3000 cm-1 
  ! is calculated following the scheme by Chou et al. (2001). 
  ! But absorptions by CFC and weak bands desinated as band 10 in Chou et al. (2001) are neglected.
  !
  !== References
  !
  !  Chou, M.-D., M. J. Suarez, X.-Z. Liang, and M. M.-H. Yan, 
  !    A thermal infrared radiation parameterization for atmospheric studies, 
  !    NASA Technical Report Series on Global Modeling and Data Assimilation, 
  !    19, NASA/TM-2001-104606, 2001.
  !
  !== Procedures List
  !
  ! RadEarthLWV24Flux :: 放射フラックスの計算
  ! ------------      :: ------------
  ! RadEarthLWV24Flux :: Calculate radiation flux
  !
  !== NAMELIST
  !
  ! NAMELIST#rad_Earth_LW_V2_4_nml
  !

  ! USE statements
  !

  ! 
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

  ! 物理・数学定数設定
  ! Physical and mathematical constants settings
  !
  use constants0, only: &
    & PI                    ! $ \pi $ .
                            ! 円周率.  Circular constant

  ! 
  ! Grid points settings
  !
  use gridset, only: imax, & ! 
                             ! Number of grid points in longitude
    &                jmax, & ! 
                             ! Number of grid points in latitude
    &                kmax    ! 
                             ! Number of vertical level


  ! Declaration statements
  !
  implicit none
  private

  !
  ! Public procedure
  !
  public :: RadEarthLWV24Flux
  public :: RadEarthLWV24Init


  character(*), parameter:: module_name = 'rad_Earth_LW_V2_4'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_Earth_LW_V2_4.f90,v 1.3 2015/01/29 12:06:43 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version


  logical , save :: FlagHighAlt


  integer , parameter :: nbmax = 10
  real(DP), save      :: aa_BandParam(1:2, 1:nbmax)

  real(DP), allocatable, save :: xyrra_TransSaved (:,:,:,:,:)
  real(DP), allocatable, save :: xyrr_TransMASaved(:,:,:,:)

  real(DP), parameter :: DiffFactor = 1.66d0



  ! MEMO:
  ! Bands range from 0 to 3000 cm-1.
  !
  data aa_BandParam &
    & / &
    &    0.0d0,  340.0d0, & ! 1:H2O
    &  340.0d0,  540.0d0, & ! 2:H2O
    &  540.0d0,  800.0d0, & ! 3:H2O + CO2
    &  800.0d0,  980.0d0, & ! 4:H2O (+ CO2)
    &  980.0d0, 1100.0d0, & ! 5:H2O (+ CO2) + O3
    & 1100.0d0, 1215.0d0, & ! 6:H2O + N2O + CH4
    & 1215.0d0, 1380.0d0, & ! 7:H2O + N2O + CH4
    & 1380.0d0, 1900.0d0, & ! 8:H2O
    & 1900.0d0, 3000.0d0, & ! 9:H2O
    &  540.0d0,  620.0d0  & !10:H2O + CO2, N2O  ! This band is not used, now.
    & /


  logical , save:: flag_save_time


  real(DP), save:: IntTimeSave
                              ! 長波フラックスを計算する時間間隔.
                              ! Interval time of long wave flux calculation
  real(DP), save:: PrevTimeSave
                              ! 前回長波フラックスを計算した時刻.
                              ! Time when long wave flux is calculated


  logical              , save:: FlagTransSaved
  data FlagTransSaved / .false. /



  ! Variables for integration of Planc function by using a pre-calculated table.
  !
  integer , save              :: ntmax
  real(DP), save, allocatable :: a_TableTemp   (:)
  real(DP), save              :: TableTempMin
  real(DP), save              :: TableTempMax
  real(DP), save              :: TableTempIncrement
  real(DP), save, allocatable :: aa_TableIPF   (:,:)
  real(DP), save, allocatable :: aa_TableIDPFDT(:,:)



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



contains

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

  subroutine RadEarthLWV24Flux(                                &
    & xyz_DelCO2Mass,                                          & ! (in )
    & xyz_DelH2OVapMass, xyz_DelH2OLiqMass, xyz_DelH2OSolMass, & ! (in )
    & xyz_DelO3Mass,                                           & ! (in )
    & xyz_DelN2OMass, xyz_DelCH4Mass,                          & ! (in )
    & xyz_Press, xyz_Temp, xy_SurfTemp,                        & ! (in )
    & xyz_QCO2, xyz_QH2OVap,                                   & ! (in )
    & xyz_QN2O, xyz_QCH4,                                      & ! (in )
    & xyz_CloudCover,                                          & ! (in )
    & xyz_CloudWatREff, xyz_CloudIceREff,                      & ! (in )
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                          & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   & ! (out)
    & )


    ! USE statements
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻.
                              ! Time of step $ t $.
      & TimesetClockStart, TimesetClockStop


!!$    ! Chou et al (1991) による長波放射モデル
!!$    ! Long radiation model described by Chou et al (1991)
!!$    !
!!$    use rad_C1991, only :               &
!!$      & RadC1991CalcTransMAH2O

    ! Chou and Kouvaris (1991) による長波放射モデル
    ! Long radiation model described by Chou and Kouvaris (1991)
    !
    use rad_CK1991, only : RadCK1991CalcTrans

    ! Chou et al. (2001) による長波放射モデル
    ! Long radiation model described by Chou et al. (2001)
    !
    use rad_C2001, only :          &
      & RadC2001CalcTransBand3CO2, &
      & RadC2001CalcTransBand3H2O, &
      & RadC2001CalcTrans,         &
      & RadC2001ReduceCloudOptDep, &
      & RadC2001CalcCloudOptProp , &
      & RadC2001CalcIntegratedPF2D, &
      & RadC2001CalcIntegratedPF3D

    ! 散乱を無視した放射伝達方程式
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScat !, RadRTENonScatAnotherForm

    ! 雲関系ルーチン
    ! Cloud-related routines
    !
    use cloud_utils, only :         &
      & CloudUtilsLocalizeCloud,    &
      & CloudUtilsCalcOverlapCloudTrans


    real(DP), intent(in ):: xyz_DelCO2Mass   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelH2OVapMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelH2OLiqMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelH2OSolMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelO3Mass    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelN2OMass   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelCH4Mass   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_Press        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_Temp         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xy_SurfTemp      (0:imax-1, 1:jmax)
    real(DP), intent(in ):: xyz_QCO2         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QH2OVap      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QN2O         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QCH4         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudCover   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudWatREff (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudIceREff (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyr_RadLUwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyr_RadLDwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP), intent(out):: xyra_DelRadLDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


    !
    ! Work variables
    !
    real(DP) :: xy_SurfEmis       (0:imax-1, 1:jmax)

    real(DP) :: xyz_CloudREff     (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudExtCoef  (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_CloudWatSSA      (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudIceSSA      (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudWatAF       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudIceAF       (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_DelCloudWatOptDep(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelCloudIceOptDep(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_TransCloudOneLayer (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyrr_TransCloud        (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    real(DP) :: xyrr_Trans            (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    real(DP) :: xyrr_TransEach        (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    real(DP) :: xyz_IntPF   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_SurfIntPF(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT0(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT1(0:imax-1, 1:jmax)


!!$    real(DP) :: xyr_RadFluxMA    (0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP) :: xyra_DelRadFluxMA(0:imax-1, 1:jmax, 0:kmax, 0:1)


    real(DP) :: xyr_RadUwFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_RadDwFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyra_DelRadUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP) :: xyra_DelRadDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)


    integer  :: n


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_Earth_LW_V2_4_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! Check whether the transmittance is saved or not
    !
    if (  .not. FlagTransSaved ) then
      if ( TimeN - PrevTimeSave >= IntTimeSave ) then
        call MessageNotify( 'M', module_name, &
          & 'Transmittance is not saved, but criterion for transmittance calculation is met.' )
      else
        call MessageNotify( 'M', module_name, &
          & 'Transmittance is not saved, and criterion for transmittance calculation ' &
          & // 'is not met. However, transmittance will be calculated.' )
      end if
    end if


    if ( ( TimeN - PrevTimeSave >= IntTimeSave ) .or. ( .not. FlagTransSaved ) ) then

!!$      write( 6, * ) 'CalcTrans'

      if ( .not. FlagTransSaved ) then
        PrevTimeSave = TimeN
      else
        PrevTimeSave = PrevTimeSave + IntTimeSave
      end if

      FlagTransSaved = .true.


      LOOP_band_trans: do n = 1, nbmax

        xyrr_Trans = 1.0_DP

        if ( n == nbmax ) then

          ! Now, nothing is done when n = nbmax.

        else if ( n == 3 ) then
          ! 540-800 cm-1

          !   Calculation of H2O line and continuum transmittance
          call RadC2001CalcTransBand3H2O(                          &
            & xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
            & xyrr_TransEach                                       & ! (out)
            & )
          xyrr_Trans = xyrr_Trans * xyrr_TransEach
          !   Calculation of CO2 transmittance
          if ( FlagHighAlt ) then
            ! Transmittance calculation for middle atmospehre as well as lower atmosphere
            call RadCK1991CalcTrans(                 &
              & xyz_DelCO2Mass, xyz_Press, xyz_Temp, & ! (in)
              & 'CO2',                               & ! (in)
              & xyrr_TransEach                       & ! (out)
              & )
          else
            ! Transmittance calculation for lower atmoshere
            call RadC2001CalcTransBand3CO2(                       &
              & xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2,    & ! (in)
              & xyrr_TransEach                                    & ! (out)
              & )
          end if
          xyrr_Trans = xyrr_Trans * xyrr_TransEach

        else

          !   Calculation of H2O continuum transmittance
          call RadC2001CalcTrans(                                  &
            & 'H2OCont', n,                                        & ! (in)
            & xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
            & xyrr_TransEach                                       & ! (out)
            & )
          xyrr_Trans = xyrr_Trans * xyrr_TransEach
          !   Calculation of H2O line transmittance
          call RadC2001CalcTrans(                                  &
            & 'H2OLine', n,                                        & ! (in)
            & xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
            & xyrr_TransEach                                       & ! (out)
            & )
          xyrr_Trans = xyrr_Trans * xyrr_TransEach

          ! Calculation of O3 transmittance
          if ( n == 5 ) then
            ! 980-1100 cm-1
            call RadCK1991CalcTrans(                &
              & xyz_DelO3Mass, xyz_Press, xyz_Temp, & ! (in)
              & 'O3',                               & ! (in)
              & xyrr_TransEach                      & ! (out)
              & )
            xyrr_Trans = xyrr_Trans * xyrr_TransEach
          end if

          ! Calculation of N2O transmittance
          if ( ( n == 6 ) .or. ( n == 7 ) .or. ( n == 10 ) ) then
            if ( maxval( xyz_DelN2OMass ) > 0.0_DP ) then
              call RadC2001CalcTrans(                                  &
                & 'N2OLine', n,                                        & ! (in)
                & xyz_DelN2OMass, xyz_Press, xyz_Temp, xyz_QN2O,       & ! (in)
                & xyrr_TransEach                                       & ! (out)
                & )
              xyrr_Trans = xyrr_Trans * xyrr_TransEach
            end if
          end if

          ! Calculation of CH4 transmittance
          if ( ( n == 6 ) .or. ( n == 7 ) ) then
            if ( maxval( xyz_DelCH4Mass ) > 0.0_DP ) then
              call RadC2001CalcTrans(                                  &
                & 'CH4Line', n,                                        & ! (in)
                & xyz_DelCH4Mass, xyz_Press, xyz_Temp, xyz_QCH4,       & ! (in)
                & xyrr_TransEach                                       & ! (out)
                & )
              xyrr_Trans = xyrr_Trans * xyrr_TransEach
            end if
          end if

          ! Calculation of CO2 transmittance in weak bands
          if ( ( n == 4 ) .or. ( n == 5 ) .or. ( n == 10 ) ) then
            call RadC2001CalcTrans(                                  &
              & 'WeakBandCO2Line', n,                                & ! (in)
              & xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2,       & ! (in)
              & xyrr_TransEach                                       & ! (out)
              & )
            xyrr_Trans = xyrr_Trans * xyrr_TransEach
          end if

        end if

        xyrra_TransSaved(:,:,:,:,n) = xyrr_Trans

      end do LOOP_band_trans

      !
      !   Calculation of transmittance of water vapor by using a method for middle 
      !   atmosphere
      !
!!$      call RadC1991CalcTransMAH2O(      &
!!$        & xyr_Press, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
!!$        & xyrr_Trans                                & ! (out)
!!$        & )
      xyrr_Trans = -1.0d100

      xyrr_TransMASaved = xyrr_Trans

    end if



    !
    ! Calculate radiative flux
    !

    xy_SurfEmis = 1.0_DP

    xyr_RadLUwFlux     = 0.0_DP
    xyr_RadLDwFlux     = 0.0_DP
    xyra_DelRadLUwFlux = 0.0_DP
    xyra_DelRadLDwFlux = 0.0_DP


    LOOP_band_RTE : do n = 1, nbmax


      xyz_CloudREff = xyz_CloudWatREff
      call RadC2001CalcCloudOptProp(                         &
        & 'Liquid', n, xyz_CloudREff,                        & ! (in)
        & xyz_CloudExtCoef, xyz_CloudWatSSA, xyz_CloudWatAF  & ! (out)
        & )
      xyz_DelCloudWatOptDep = xyz_CloudExtCoef * xyz_DelH2OLiqMass
      !
      xyz_CloudREff = xyz_CloudIceREff
      call RadC2001CalcCloudOptProp(                         &
        & 'Ice', n, xyz_CloudREff,                           & ! (in)
        & xyz_CloudExtCoef, xyz_CloudIceSSA, xyz_CloudIceAF  & ! (out)
        & )
      xyz_DelCloudIceOptDep = xyz_CloudExtCoef * xyz_DelH2OSolMass

      call RadC2001ReduceCloudOptDep(      &
        & xyz_CloudWatSSA, xyz_CloudWatAF, & ! (in)
        & xyz_DelCloudWatOptDep            & ! (inout)
        & )
      call RadC2001ReduceCloudOptDep(      &
        & xyz_CloudIceSSA, xyz_CloudIceAF, & ! (in)
        & xyz_DelCloudIceOptDep            & ! (inout)
        & )
      !
      call CloudUtilsLocalizeCloud(  &
        & xyz_CloudCover,            & ! (in   )
        & xyz_DelCloudWatOptDep      & ! (inout)
        & )
      call CloudUtilsLocalizeCloud(  &
        & xyz_CloudCover,            & ! (in   )
        & xyz_DelCloudIceOptDep      & ! (inout)
        & )
      !
      xyz_TransCloudOneLayer = &
        & exp( - ( xyz_DelCloudWatOptDep + xyz_DelCloudIceOptDep ) * DiffFactor )
      !
      call CloudUtilsCalcOverlapCloudTrans(       &
        & xyz_TransCloudOneLayer, xyz_CloudCover, & ! (in)
        & xyrr_TransCloud                         & ! (out)
        & )


      ! Now, nothing is done when n = nbmax.
      if ( n == nbmax ) cycle


      xyrr_Trans = xyrra_TransSaved(:,:,:,:,n)

      xyrr_Trans = xyrr_Trans * xyrr_TransCloud



      call CalcIntegratedPFWithTable2D( &
        & n, xy_SurfTemp,                         &
        & xy_SurfIntPF,                           &
        & 1, jmax                                 &
        & )
      call CalcIntegratedPFWithTable3D( &
        & n, kmax, xyz_Temp,                    &
        & xyz_IntPF,                            &
        & 1, jmax                               &
        & )

      call CalcIntegratedPFWithTable2D( &
        & n, xy_SurfTemp,                         &
        & xy_IntDPFDT0,                           &
        & 1, jmax,                                &
        & .true.                                  &
        & )
      call CalcIntegratedPFWithTable2D( &
        & n, xyz_Temp(:,:,1),                     &
        & xy_IntDPFDT1,                           &
        & 1, jmax,                                &
        & .true.                                  &
        & )

      xy_SurfIntPF = xy_SurfEmis * PI * xy_SurfIntPF
      xyz_IntPF    =               PI * xyz_IntPF
      xy_IntDPFDT0 = xy_SurfEmis * PI * xy_IntDPFDT0
      xy_IntDPFDT1 =               PI * xy_IntDPFDT1



      ! Lines below are under testing.
      !
!!$      xy_SurfIntPF2 = xy_SurfIntPF
!!$      xyz_IntPF2    = xyz_IntPF
!!$      xy_IntDPFDT02 = xy_IntDPFDT0
!!$      xy_IntDPFDT12 = xy_IntDPFDT1
!!$
!!$
!!$      call RadC2001CalcIntegratedPF2D(       &
!!$        & n, xy_SurfTemp,                    &
!!$        & xy_SurfIntPF                       &
!!$        & )
!!$      call RadC2001CalcIntegratedPF3D(       &
!!$        & n, kmax, xyz_Temp,                 &
!!$        & xyz_IntPF                          &
!!$        & )
!!$
!!$      call RadC2001CalcIntegratedPF2D(       &
!!$        & n, xy_SurfTemp,                    &
!!$        & xy_IntDPFDT0,                      &
!!$        & .true.                             &
!!$        & )
!!$      call RadC2001CalcIntegratedPF2D(       &
!!$        & n, xyz_Temp(:,:,1),                &
!!$        & xy_IntDPFDT1,                      &
!!$        & .true.                             &
!!$        & )
!!$
!!$      xy_SurfIntPF = xy_SurfEmis * xy_SurfIntPF
!!$      xyz_IntPF    =               xyz_IntPF
!!$      xy_IntDPFDT0 = xy_SurfEmis * xy_IntDPFDT0
!!$      xy_IntDPFDT1 =               xy_IntDPFDT1
!!$
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          write( 20+n, * ) xy_SurfTemp(i,j), xy_SurfIntPF2(i,j), xy_SurfIntPF(i,j), &
!!$            & xy_IntDPFDT02(i,j), xy_IntDPFDT02(i,j)
!!$        end do
!!$      end do
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          write( 40+n, * ) xyz_Temp(i,j,1), xy_IntDPFDT12(i,j), xy_IntDPFDT12(i,j)
!!$        end do
!!$      end do
!!$      do k = 1, kmax
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            write( 60+n, * ) xyz_Temp(i,j,k), xyz_IntPF2(i,j,k), xyz_IntPF(i,j,k)
!!$          end do
!!$        end do
!!$      end do
!!$      call flush( 20+n )
!!$      call flush( 40+n )
!!$      call flush( 60+n )



!!$      call OLD_RadRTENonScat(                                      &
!!$        & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0, & ! (in)
!!$        & xyrr_Trans,                                          & ! (in)
!!$        & xyr_RadFlux, xyra_DelRadFlux                         & ! (out)
!!$        & )

      call RadRTENonScat(                                                &
!!$      call RadRTENonScatAnotherForm( &
        & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0, & ! (in)
        & xyrr_Trans,                                          & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux,                      & ! (out)
        & xyra_DelRadUwFlux, xyra_DelRadDwFlux               & ! (out)
        & )



!!$      i = 0
!!$      j = jmax/2+1
!!$      do k = 0, kmax
!!$        write( 6, * ) k, xyr_RadFlux(i,j,k), &
!!$          & xyr_RadFlux(i,j,k) - ( xyr_RadUwFlux(i,j,k) - xyr_RadDwFlux(i,j,k) )
!!$      end do
!!$      do k = 0, kmax
!!$        write( 6, * ) k, &
!!$          & xyra_DelRadFlux(i,j,k,0), &
!!$          & xyra_DelRadFlux(i,j,k,0) - ( xyra_DelRadUwFlux(i,j,k,0) - xyra_DelRadDwFlux(i,j,k,0) ), &
!!$          & xyra_DelRadFlux(i,j,k,1), &
!!$          & xyra_DelRadFlux(i,j,k,1) - ( xyra_DelRadUwFlux(i,j,k,1) - xyra_DelRadDwFlux(i,j,k,1) )
!!$      end do



!!$      if ( ( n == 1 ) .or. ( n == 2 ) .or. ( n == 9 ) ) then
!!$        !
!!$        ! For bands 0-340, 340-540, 1380-1900
!!$        ! merge with flux calculated with a method for middle atmosphere
!!$        !
!!$
!!$        xyrr_Trans = xyrr_TransMASaved
!!$
!!$        call RadELWV22IntegRTE(                               &
!!$          & n,                                                & ! (in )
!!$          & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans,   & ! (in )
!!$          & xyr_RadFluxMA, xyra_DelRadFluxMA                  & ! (out)
!!$          & )
!!$
!!$        call RadDcpamELWV22CutMergeFlux(          &
!!$          & xyr_Press,                            & ! (in)
!!$          & xyr_RadFlux, xyra_DelRadFlux,         & ! (inout)
!!$          & xyr_RadFluxMA, xyra_DelRadFluxMA      & ! (in) optional
!!$          & )
!!$
!!$      else if ( ( n == 4 ) .or. ( n == 6 ) .or. ( n == 8 ) ) then
!!$        !
!!$        ! For bands 800-980, 1100-1380, 1900-3000
!!$        ! flux above a pressure level is modified to be constant
!!$        !
!!$
!!$        call RadDcpamELWV22CutMergeFlux(          &
!!$          & xyr_Press,                            & ! (in)
!!$          & xyr_RadFlux, xyra_DelRadFlux          & ! (inout)
!!$          & )
!!$
!!$      end if


      xyr_RadLUwFlux     = xyr_RadLUwFlux     + xyr_RadUwFlux
      xyr_RadLDwFlux     = xyr_RadLDwFlux     + xyr_RadDwFlux
      xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadUwFlux
      xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadDwFlux

    end do LOOP_band_RTE





!!$    i = 0
!!$    j = jmax / 2 + 1
!!$    write( 73, * ) xy_SurfTemp(i,j), 0.0d0, 0.0d0, xyr_Press(i,j,0)
!!$    do k = 1, kmax
!!$      write( 73, * ) xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), xyz_QO3(i,j,k), &
!!$        & xyz_Press(i,j,k)
!!$    end do
!!$    call flush( 73 )
!!$
!!$    i = 0
!!$    j = jmax / 2 + 1
!!$    do k = 1, kmax
!!$      write( 83, * ) &
!!$        & + (     xyr_RadLFlux(i,j,k-1) - xyr_RadLFlux(i,j,k) )  &
!!$        &     / ( xyr_Press(i,j,k-1)    - xyr_Press(i,j,k) )     &
!!$        &     / 1004.6 * 9.8, &
!!$        & xyz_Press(i,j,k)
!!$    end do
!!$    call flush( 83 )
!!$
!!$    i = 0
!!$    j = jmax / 2 + 1
!!$    do k = 0, kmax
!!$      write( 93, * ) xyr_RadLFlux(i,j,k), xyr_Press(i,j,k)
!!$    end do
!!$    call flush( 93 )
!!$    stop


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine RadEarthLWV24Flux

  !--------------------------------------------------------------------------------------
!!$
!!$  subroutine RadiationDcpamELWV23CutMergeFlux( &
!!$    & xyr_Press,                              & ! (in)
!!$    & xyr_RadLFlux, xyra_DelRadLFlux,         & ! (inout)
!!$    & xyr_RadLFluxMA, xyra_DelRadLFluxMA      & ! (in) optional
!!$    & )
!!$    !
!!$    ! Radiative flux above a pressure level is modified to be constant or is merged with
!!$    ! that in middle atmosphere
!!$    !
!!$
!!$    ! USE statements
!!$    !
!!$
!!$    !
!!$    ! Grid points settings
!!$    !
!!$    use gridset, only: imax, & ! 
!!$                               ! Number of grid points in longitude
!!$      &                jmax, & ! 
!!$                               ! Number of grid points in latitude
!!$      &                kmax    ! 
!!$                               ! Number of vertical level
!!$
!!$    real(DP), intent(in   )           :: xyr_Press         (0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP), intent(inout)           :: xyr_RadLFlux      (0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP), intent(inout)           :: xyra_DelRadLFlux  (0:imax-1, 1:jmax, 0:kmax, 0:1)
!!$    real(DP), intent(in   ), optional :: xyr_RadLFluxMA    (0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP), intent(in   ), optional :: xyra_DelRadLFluxMA(0:imax-1, 1:jmax, 0:kmax, 0:1)
!!$
!!$    !
!!$    ! Work variables
!!$    !
!!$    real(DP), parameter :: BoundaryPress = 30.0d2
!!$    integer             :: xy_kcut            (0:imax-1, 1:jmax)
!!$    real(DP)            :: xy_BoundaryFlux    (0:imax-1, 1:jmax)
!!$    real(DP)            :: xya_BoundaryDelFlux(0:imax-1, 1:jmax, 0:1)
!!$    integer             :: i
!!$    integer             :: j
!!$    integer             :: k
!!$
!!$
!!$    do k = 0, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xyr_Press(i,j,k) >= BoundaryPress ) then
!!$            xy_kcut(i,j) = k
!!$          end if
!!$        end do
!!$      end do
!!$    end do
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_kcut(i,j) == kmax ) then
!!$          xy_kcut(i,j) = kmax - 1
!!$        end if
!!$      end do
!!$    end do
!!$
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        k = xy_kcut(i,j)
!!$        xy_BoundaryFlux(i,j) =                                              &
!!$          &   ( xyr_RadLFlux(i,j,k+1) - xyr_RadLFlux(i,j,k) )               &
!!$          & / log( ( xyr_Press(i,j,k+1) + 1.0d-100 ) / ( xyr_Press(i,j,k) + 1.0d-100 ) )&
!!$          & * log( BoundaryPress                     / ( xyr_Press(i,j,k) + 1.0d-100 ) )&
!!$          & + xyr_RadLFlux(i,j,k)
!!$        xya_BoundaryDelFlux(i,j,0) =                                        &
!!$          &   ( xyra_DelRadLFlux(i,j,k+1,0) - xyra_DelRadLFlux(i,j,k,0) )   &
!!$          & / log( ( xyr_Press(i,j,k+1) + 1.0d-100 ) / ( xyr_Press(i,j,k) + 1.0d-100 ) )&
!!$          & * log( BoundaryPress                     / ( xyr_Press(i,j,k) + 1.0d-100 ) )&
!!$          & + xyra_DelRadLFlux(i,j,k,0)
!!$        xya_BoundaryDelFlux(i,j,1) =                                        &
!!$          &   ( xyra_DelRadLFlux(i,j,k+1,1) - xyra_DelRadLFlux(i,j,k,1) )   &
!!$          & / log( ( xyr_Press(i,j,k+1) + 1.0d-100 ) / ( xyr_Press(i,j,k) + 1.0d-100 ) )&
!!$          & * log( BoundaryPress                     / ( xyr_Press(i,j,k) + 1.0d-100 ) )&
!!$          & + xyra_DelRadLFlux(i,j,k,1)
!!$      end do
!!$    end do
!!$
!!$    do k = 0, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xyr_Press(i,j,k) < BoundaryPress ) then
!!$            xyr_RadLFlux(i,j,k) = xy_BoundaryFlux(i,j)
!!$          end if
!!$        end do
!!$      end do
!!$    end do
!!$    do k = 0, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if ( xyr_Press(i,j,k) < BoundaryPress ) then
!!$            xyra_DelRadLFlux(i,j,k,0) = xya_BoundaryDelFlux(i,j,0)
!!$            xyra_DelRadLFlux(i,j,k,1) = xya_BoundaryDelFlux(i,j,1)
!!$          end if
!!$        end do
!!$      end do
!!$    end do
!!$
!!$
!!$    if ( present( xyr_RadLFluxMA ) ) then
!!$
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          k = xy_kcut(i,j)
!!$          xy_BoundaryFlux(i,j) =                                              &
!!$            &   ( xyr_RadLFluxMA(i,j,k+1) - xyr_RadLFluxMA(i,j,k) )           &
!!$            & / log( ( xyr_Press(i,j,k+1) + 1.0d-100 ) / ( xyr_Press(i,j,k) + 1.0d-100 ) )   &
!!$            & * log( BoundaryPress                     / ( xyr_Press(i,j,k) + 1.0d-100 ) )   &
!!$            & + xyr_RadLFluxMA(i,j,k)
!!$        end do
!!$      end do
!!$
!!$      do k = 0, kmax
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            if ( xyr_Press(i,j,k) < BoundaryPress ) then
!!$              xyr_RadLFlux(i,j,k) = xyr_RadLFlux(i,j,k) &
!!$                & + xyr_RadLFluxMA(i,j,k) - xy_BoundaryFlux(i,j)
!!$            end if
!!$          end do
!!$        end do
!!$      end do
!!$
!!$    end if
!!$
!!$    if ( present( xyra_DelRadLFluxMA ) ) then
!!$
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          k = xy_kcut(i,j)
!!$          xya_BoundaryDelFlux(i,j,0) =                                            &
!!$            &   ( xyra_DelRadLFluxMA(i,j,k+1,0) - xyra_DelRadLFluxMA(i,j,k,0) )   &
!!$            & / log( ( xyr_Press(i,j,k+1) + 1.0d-100 ) / ( xyr_Press(i,j,k) + 1.0d-100 ) )       &
!!$            & * log( BoundaryPress                     / ( xyr_Press(i,j,k) + 1.0d-100 ) )       &
!!$            & + xyra_DelRadLFluxMA(i,j,k,0)
!!$          xya_BoundaryDelFlux(i,j,1) =                                            &
!!$            &   ( xyra_DelRadLFluxMA(i,j,k+1,1) - xyra_DelRadLFluxMA(i,j,k,1) )   &
!!$            & / log( ( xyr_Press(i,j,k+1) + 1.0d-100 ) / ( xyr_Press(i,j,k) + 1.0d-100 ) )       &
!!$            & * log( BoundaryPress                     / ( xyr_Press(i,j,k) + 1.0d-100 ) )       &
!!$            & + xyra_DelRadLFluxMA(i,j,k,1)
!!$        end do
!!$      end do
!!$
!!$      do k = 0, kmax
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            if ( xyr_Press(i,j,k) < BoundaryPress ) then
!!$              xyra_DelRadLFlux(i,j,k,0) = xyra_DelRadLFlux(i,j,k,0) &
!!$                & + xyra_DelRadLFluxMA(i,j,k,0) - xya_BoundaryDelFlux(i,j,0)
!!$              xyra_DelRadLFlux(i,j,k,1) = xyra_DelRadLFlux(i,j,k,1) &
!!$                & + xyra_DelRadLFluxMA(i,j,k,1) - xya_BoundaryDelFlux(i,j,1)
!!$            end if
!!$          end do
!!$        end do
!!$      end do
!!$
!!$    end if
!!$
!!$
!!$  end subroutine RadiationDcpamELWV23CutMergeFlux

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

  subroutine RadEarthLWV24Init( &
    & FlagSnow                  &
    & )

    ! USE statements
    !

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

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 暦と日時の取り扱い
    ! Calendar and Date handler
    !
    use dc_calendar, only: DCCalConvertByUnit

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! Chou and Kouvaris (1991) による長波放射モデル
    ! Long radiation model described by Chou and Kouvaris (1991)
    !
    use rad_CK1991, only : RadCK1991Init

    ! Chou et al. (2001) による長波放射モデル
    ! Long radiation model described by Chou et al. (2001)
    !
    use rad_C2001, only : RadC2001Init

    ! 散乱を無視した放射伝達方程式
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScatInit

    ! 雲関系ルーチン
    ! Cloud-related routines
    !
    use cloud_utils, only : CloudUtilsInit


    logical, intent(in) :: FlagSnow


    real(DP)          :: DelTimeCalcTransValue
    character(STRING) :: DelTimeCalcTransUnit

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    integer :: n


    namelist /rad_Earth_LW_V2_4_nml/ &
      & FlagHighAlt,              &
      & DelTimeCalcTransValue,    &
      & DelTimeCalcTransUnit,     &
      & flag_save_time


    if ( rad_Earth_LW_V2_4_inited ) return


    FlagHighAlt           = .false.

    DelTimeCalcTransValue = 3.0
    DelTimeCalcTransUnit  = 'hrs'
    flag_save_time        = .false.


    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml,                          & ! (in)
        & nml = rad_Earth_LW_V2_4_nml,         & ! (out)
        & iostat = iostat_nml )                  ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if

    ! Handle interval time
    !
    IntTimeSave = DCCalConvertByUnit( DelTimeCalcTransValue, DelTimeCalcTransUnit, 'sec' ) ! (in)




    do n = 1, nbmax
      ! unit conversion from (cm-1) to (m-1)
      aa_BandParam(1,n) = aa_BandParam(1,n) * 1.0d2
      aa_BandParam(2,n) = aa_BandParam(2,n) * 1.0d2
    end do


    ! allocate a variable for saving transmittance
    !
    allocate( xyrra_TransSaved (0:imax-1,1:jmax,0:kmax,0:kmax,1:nbmax) )
    allocate( xyrr_TransMASaved(0:imax-1,1:jmax,0:kmax,0:kmax)         )


    call RadEarthLWV24PrepPFTable


    ! Initialization of modules used in this module
    !

    ! Chou and Kouvaris (1991) による長波放射モデル
    ! Long radiation model described by Chou and Kouvaris (1991)
    !
    call RadCK1991Init

    ! Chou et al. (2001) による長波放射モデル
    ! Long radiation model described by Chou et al. (2001)
    !
    call RadC2001Init

    ! 散乱を無視した放射伝達方程式
    ! Radiative transfer equation without considering scattering
    !
    call RadRTENonScatInit

    ! 雲関系ルーチン
    ! Cloud-related routines
    !
    call CloudUtilsInit( &
      & FlagSnow         &
      & )


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  FlagHighAlt       = %b', &
      & l = (/ FlagHighAlt /) )
    call MessageNotify( 'M', module_name, '  DelTimeCalcTrans  = %f [%c]', &
      & d = (/ DelTimeCalcTransValue /), c1 = trim( DelTimeCalcTransUnit ) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    rad_Earth_LW_V2_4_inited = .true.

  end subroutine RadEarthLWV24Init

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

  subroutine RadEarthLWV24PrepPFTable

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

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

    ! プランク関数の計算
    ! Calculate Planck function
    !
    use planck_func, only : PF, DPFDT, Integ_PF_GQ_Array2D, Integ_DPFDT_GQ_Array2D

    integer , parameter :: NGaussQuad = 5
    logical             :: FlagCheckLoopExit
    real(DP)            :: xy_TempTMP   (0:imax-1, 1:jmax)
    real(DP)            :: xy_PF        (0:imax-1, 1:jmax)
    real(DP)            :: xy_DPFDT     (0:imax-1, 1:jmax)
    real(DP)            :: xy_PFTable   (0:imax-1, 1:jmax)
    real(DP)            :: xy_DPFDTTable(0:imax-1, 1:jmax)
    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 :: GQP(:)
    real(DP)      , allocatable :: GQW(:)


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


    ! Preparation of tables for calculation of Plank function
    !
    TableTempMin       =  50.0d0
    TableTempMax       = 600.0d0
    TableTempIncrement =   0.1d0
    ntmax              = ( TableTempMax - TableTempMin ) / TableTempIncrement + 1
    allocate( a_TableTemp   (1:ntmax) )
    allocate( aa_TableIPF   (1:ntmax, 1:nbmax) )
    allocate( aa_TableIDPFDT(1:ntmax, 1:nbmax) )

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


    aa_TableIPF   (:,:) = 0.0d0
    aa_TableIDPFDT(:,:) = 0.0d0

    allocate( GQP(1:NGaussQuad) )
    allocate( GQW(1:NGaussQuad) )
    do n = 1, nbmax
      call GauLeg( &
        & aa_BandParam(1,n), aa_BandParam(2,n), NGaussQuad, & ! (in )
        & GQP, GQW                                          & ! (out)
        & )
      do m = 1, ntmax
        do l = 1, NGaussQuad
          aa_TableIPF   (m,n) = &
            & aa_TableIPF   (m,n) + PF   ( GQP(l), a_TableTemp(m) ) * GQW(l)
          aa_TableIDPFDT(m,n) = &
            & aa_TableIDPFDT(m,n) + DPFDT( GQP(l), a_TableTemp(m) ) * GQW(l)
        end do
      end do
    end do
    deallocate( GQP )
    deallocate( GQW )


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

    !      This routine is called once here, to initialize a pre-calculated table.
    n = 1
    xy_TempTMP = 300.0d0
    call CalcIntegratedPFWithTable2D( &
      & n, xy_TempTMP,                &
      & xy_PFTable,                   &
      & 1, jmax,                      &
      & .false.                       &
      & )

    do n = 1, nbmax

      FlagCheckLoopExit = .false.
      l = 1
      do

        do j = 1, jmax
          do i = 0, imax-1
            xy_TempTMP(i,j) = &
              &   a_TableTemp(1)                                                     &
              & + ( a_TableTemp(2) - a_TableTemp(1) ) * 0.5d0                        &
              & + ( a_TableTemp(2) - a_TableTemp(1) ) &
              &     * ( imax * jmax * ( l - 1 ) + imax * ( j - 1 ) + i )
          end do
        end do

        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_TempTMP(i,j) > a_TableTemp(ntmax) ) then
              xy_TempTMP(i,j) = a_TableTemp(ntmax)
              FlagCheckLoopExit = .true.
            end if
          end do
        end do


        call Integ_PF_GQ_Array2D(                             &
          & aa_BandParam(1,n), aa_BandParam(2,n), NGaussQuad, &
          & 0, imax-1, 1, jmax, xy_TempTMP,                   &
          & xy_PF                                             &
          & )
        call Integ_DPFDT_GQ_Array2D(                          &
          & aa_BandParam(1,n), aa_BandParam(2,n), NGaussQuad, & ! (in )
          & 0, imax-1, 1, jmax, xy_TempTMP,                   & ! (in )
          & xy_DPFDT                                          & ! (out)
          & )

        call CalcIntegratedPFWithTable2D( &
          & n, xy_TempTMP,                &
          & xy_PFTable,                   &
          & 1, jmax,                      &
          & .false.                       &
          & )
        call CalcIntegratedPFWithTable2D( &
          & n, xy_TempTMP,                &
          & xy_DPFDTTable,                &
          & 1, jmax,                      &
          & .true.                        &
          & )

        do j = 1, jmax
          do i = 0, imax-1
            ErrorPFInteg = &
              & abs( xy_PF   (i,j) - xy_PFTable   (i,j) ) / xy_PF   (i,j)
            if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
              call MessageNotify( 'E', module_name, 'Error of integrated PF, %f, is greater than threshold, %f, in band %d.', &
                & d = (/ ErrorPFInteg, ThresholdErrorPFInteg /), i = (/n/) )
            end if
            ErrorPFInteg = &
              & abs( xy_DPFDT(i,j) - xy_DPFDTTable(i,j) ) / xy_DPFDT(i,j)
            if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
              call MessageNotify( 'E', module_name, 'Error of integrated DPFDT, %f, is greater than threshold, %f, in band %d', &
                & d = (/ ErrorPFInteg, ThresholdErrorPFInteg /), i = (/n/) )
            end if
          end do
        end do

        if ( FlagCheckLoopExit ) exit
        l = l + 1
      end do

    end do

  end subroutine RadEarthLWV24PrepPFTable

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

  subroutine CalcIntegratedPFWithTable2D( &
    & iband, xy_Temp,                     &
    & xy_IntegPF,                         &
    & js, je,                             &
    & flag_DPFDT                          &
    & )

    ! USE statements
    !

    integer , intent(in )           :: iband
    real(DP), intent(in )           :: xy_temp   (0:imax-1, 1:jmax)
    real(DP), intent(out)           :: xy_IntegPF(0:imax-1, 1:jmax)
    integer , intent(in )           :: js
    integer , intent(in )           :: je
    logical , intent(in ), optional :: flag_DPFDT

    !
    ! local variables
    !
    real(DP) :: xyz_Temp   (0:imax-1, 1:jmax, 1)
    real(DP) :: xyz_IntegPF(0:imax-1, 1:jmax, 1)
    integer  :: j


    do j = js, je
      xyz_Temp(:,j,1) = xy_Temp(:,j)
    end do

    call CalcIntegratedPFWithTable3D( &
      & iband, 1, xyz_temp,                 &
      & xyz_IntegPF,                        &
      & js, je,                             &
      & flag_DPFDT                          &
      & )

    do j = js, je
      xy_IntegPF(:,j) = xyz_IntegPF(:,j,1)
    end do


  end subroutine CalcIntegratedPFWithTable2D

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

  subroutine CalcIntegratedPFWithTable3D( &
    & iband, km, xyz_temp,                &
    & xyz_IntegPF,                        &
    & js, je,                             &
    & flag_DPFDT                          &
    & )

    ! USE statements
    !

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

    integer , intent(in ) :: iband
    integer , intent(in ) :: km
    real(DP), intent(in ) :: xyz_temp   (0:imax-1, 1:jmax, 1:km)
    real(DP), intent(out) :: xyz_IntegPF(0:imax-1, 1:jmax, 1:km)
    logical , intent(in ), optional :: flag_DPFDT
    integer , intent(in )           :: js
    integer , intent(in )           :: je

    !
    ! local variables
    !
    logical                     :: local_flag_DPFDT

    integer                     :: xyz_TempIndex(0:imax-1, 1:jmax, 1:km)
    integer                     :: i
    integer                     :: j
    integer                     :: k
    integer                     :: m


    do k = 1, km
      do j = js, je
        do i = 0, imax-1

          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 = 1, km
        do j = js, je
          do i = 0, imax-1
            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) = &
              &   aa_TableIPF(m-1,iband)                           &
              &   * ( 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 ) ) )  &
              & + aa_TableIPF(m  ,iband)                           &
              &   * ( 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 ) ) )  &
              & + aa_TableIPF(m+1,iband)                           &
              &   * ( 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 = 1, km
        do j = js, je
          do i = 0, imax-1
            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_TableIDPFDT( m-1, iband )

            xyz_IntegPF(i,j,k) = &
              &   aa_TableIDPFDT(m-1,iband)                        &
              &   * ( 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 ) ) )  &
              & + aa_TableIDPFDT(m  ,iband)                        &
              &   * ( 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 ) ) )  &
              & + aa_TableIDPFDT(m+1,iband)                        &
              &   * ( 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 rad_Earth_LW_V2_4