!= Roewe and Liou (1978) による長波放射モデル
!
!= Long radiation model described by Roewe and Liou (1978)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: radiation_RL78.f90,v 1.6 2010-04-13 13:14:45 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module radiation_RL78
  !
  != Roewe and Liou (1978) による長波放射モデル
  !
  != Long radiation model described by Roewe and Liou (1978)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 長波放射モデル.
  !
  ! This is a model of long wave radiation. 
  !
  !== References
  !
  !  Roewe, D., and K.-N. Liou, Influence of cirrus clouds on the infrared cooling 
  !    rate in the troposphere and lower stratosphere, J. Appl. Met., 17, 92-106, 1978.
  !
  !== 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.

  ! 日付および時刻の取り扱い
  ! Date and time handler
  !
  use dc_date_types, only: DC_DIFFTIME
                              ! 日時の差を表現するデータ型.
                              ! Data type for difference about date and time

  ! Declaration statements
  !
  implicit none
  private

  ! 
  ! Public procedure
  !
  public:: RadiationRL78Flux


  character(*), parameter:: module_name = 'radiation_RL78'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: radiation_RL78.f90,v 1.6 2010-04-13 13:14:45 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

  real(DP), parameter:: AvogNum    = 6.0221415d23
  real(DP), parameter:: DiffFactor = 1.66d0


  integer , parameter:: nbmax = 32
  real(DP), save     :: a_BandParam(1:4,1:nbmax)

  real(DP), save     :: a_CO2BandParam(1:4)

!!$  integer , parameter:: nbH2Omax = 21
!!$  real(DP), save     :: a_H2OBandParam(1:4,1:nbH2Omax)
!!$  real(DP), save     :: a_H2OContParam(1:3,1:nbH2Omax)
!!$  integer , parameter:: nbO3max = 11
!!$  real(DP), save     :: a_O3BandParam(1:4,1:nbO3max)


!!$  data a_BandParam &
!!$    & /    0.0d0,   40.0d0,    0.0d0   , 0.0d0  , & !  1:nothing
!!$    !---
!!$    &     40.0d0,  160.0d0, 7240.30d0  , 0.182d0, & !  2:H2O band
!!$    &    160.0d0,  280.0d0, 6024.80d0  , 0.094d0, & !  3:H2O band
!!$    &    280.0d0,  380.0d0, 1614.10d0  , 0.081d0, & !  4:H2O band
!!$    &    380.0d0,  500.0d0,  139.03d0  , 0.080d0, & !  5:H2O band
!!$    &    500.0d0,  600.0d0,   21.64d0  , 0.068d0, & !  6:H2O band & CO2 band
!!$    &    600.0d0,  720.0d0,    2.919d0 , 0.060d0, & !  7:H2O band & CO2 band
!!$    &    720.0d0,  800.0d0,    0.386d0 , 0.059d0, & !  8:H2O band & CO2 band
!!$    !---
!!$    &    800.0d0,  900.0d0,    0.0715d0, 0.067d0, & !  9:H2O band & continuum
!!$    !---
!!$    &    900.0d0, 1000.0d0,    0.0d0   , 0.0d0  , & ! 10:H2O continuum
!!$    !---
!!$    &   1000.0d0, 1006.5d0,    3.21d-3 , 5.0d0,   & ! 11:O3 & H2O continuum
!!$    &   1006.5d0, 1013.0d0,    6.42d-3 , 5.0d0,   & ! 12:O3 & H2O continuum
!!$    &   1013.0d0, 1019.5d0,    1.28d-2 , 5.0d0,   & ! 13:O3 & H2O continuum
!!$    &   1019.5d0, 1026.0d0,    2.14d-2 , 5.5d0,   & ! 14:O3 & H2O continuum
!!$    &   1026.0d0, 1032.5d0,    2.35d-2 , 5.8d0,   & ! 15:O3 & H2O continuum
!!$    &   1032.5d0, 1039.0d0,    1.71d-2 , 8.0d0,   & ! 16:O3 & H2O continuum
!!$    &   1039.0d0, 1045.5d0,    1.18d-2 , 6.1d0,   & ! 17:O3 & H2O continuum
!!$    &   1045.5d0, 1052.0d0,    2.78d-2 , 8.4d0,   & ! 18:O3 & H2O continuum
!!$    &   1052.0d0, 1058.5d0,    3.53d-2 , 8.3d0,   & ! 19:O3 & H2O continuum
!!$    &   1058.5d0, 1065.0d0,    1.28d-2 , 6.7d0,   & ! 20:O3 & H2O continuum
!!$    !---
!!$    &   1065.0d0, 1200.0d0,    0.0d0   , 0.0d0,   & ! 21:H2O continuum
!!$    !---
!!$    &   1200.0d0, 1350.0d0,   12.65d0  , 0.089d0, & ! 22:H2O band
!!$    &   1350.0d0, 1450.0d0,  134.4d0   , 0.230d0, & ! 23:H2O band
!!$    &   1450.0d0, 1550.0d0,  632.0d0   , 0.320d0, & ! 24:H2O band
!!$    &   1550.0d0, 1650.0d0,  331.2d0   , 0.296d0, & ! 25:H2O band
!!$    &   1650.0d0, 1750.0d0,  434.1d0   , 0.452d0, & ! 26:H2O band
!!$    &   1750.0d0, 1850.0d0,  136.0d0   , 0.359d0, & ! 27:H2O band
!!$    &   1850.0d0, 1950.0d0,   35.7d0   , 0.165d0, & ! 28:H2O band
!!$    &   1950.0d0, 2050.0d0,    9.015d0 , 0.104d0, & ! 29:H2O band
!!$    &   2050.0d0, 2200.0d0,    1.529d0 , 0.116d0, & ! 30:H2O band
!!$    !---
!!$    &   2200.0d0, 9000.0d0,    0.0d0   , 0.0d0    / ! 31:nothing


  ! MEMO:
  ! Bands range from 40 to 2500 cm-1.
  !
  ! Integrated flux of Planck function below 50 cm-1 is about 0.005%. 
  ! Integrated flux of Planck function beyond 2500 cm-1 is about 0.15%. 
  ! (Sekiguchi, D-thesis, p.41)
  data a_BandParam &
    & /   40.0d0,  160.0d0, 7240.30d0  , 0.182d0, & !  1:H2O band
    &    160.0d0,  280.0d0, 6024.80d0  , 0.094d0, & !  2:H2O band
    &    280.0d0,  380.0d0, 1614.10d0  , 0.081d0, & !  3:H2O band
    &    380.0d0,  500.0d0,  139.03d0  , 0.080d0, & !  4:H2O band
    !---
    &    500.0d0,  582.0d0,   21.64d0  , 0.068d0, & !  5:H2O band
    &    582.0d0,  600.0d0,   21.64d0  , 0.068d0, & !  6:H2O band & CO2 band
    &    600.0d0,  720.0d0,    2.919d0 , 0.060d0, & !  7:H2O band & CO2 band
    &    720.0d0,  752.0d0,    0.386d0 , 0.059d0, & !  8:H2O band & CO2 band
    &    752.0d0,  800.0d0,    0.386d0 , 0.059d0, & !  9:H2O band
    !---
    &    800.0d0,  900.0d0,    0.0715d0, 0.067d0, & ! 10:H2O band & continuum
    !---
    &    900.0d0, 1000.0d0,    0.0d0   , 0.0d0  , & ! 11:H2O continuum
    !---
    &   1000.0d0, 1006.5d0,    3.21d-3 , 5.0d0,   & ! 12:O3 & H2O continuum
    &   1006.5d0, 1013.0d0,    6.42d-3 , 5.0d0,   & ! 13:O3 & H2O continuum
    &   1013.0d0, 1019.5d0,    1.28d-2 , 5.0d0,   & ! 14:O3 & H2O continuum
    &   1019.5d0, 1026.0d0,    2.14d-2 , 5.5d0,   & ! 15:O3 & H2O continuum
    &   1026.0d0, 1032.5d0,    2.35d-2 , 5.8d0,   & ! 16:O3 & H2O continuum
    &   1032.5d0, 1039.0d0,    1.71d-2 , 8.0d0,   & ! 17:O3 & H2O continuum
    &   1039.0d0, 1045.5d0,    1.18d-2 , 6.1d0,   & ! 18:O3 & H2O continuum
    &   1045.5d0, 1052.0d0,    2.78d-2 , 8.4d0,   & ! 19:O3 & H2O continuum
    &   1052.0d0, 1058.5d0,    3.53d-2 , 8.3d0,   & ! 20:O3 & H2O continuum
    &   1058.5d0, 1065.0d0,    1.28d-2 , 6.7d0,   & ! 21:O3 & H2O continuum
    !---
    &   1065.0d0, 1200.0d0,    0.0d0   , 0.0d0,   & ! 22:H2O continuum
    !---
    &   1200.0d0, 1350.0d0,   12.65d0  , 0.089d0, & ! 23:H2O band
    &   1350.0d0, 1450.0d0,  134.4d0   , 0.230d0, & ! 24:H2O band
    &   1450.0d0, 1550.0d0,  632.0d0   , 0.320d0, & ! 25:H2O band
    &   1550.0d0, 1650.0d0,  331.2d0   , 0.296d0, & ! 26:H2O band
    &   1650.0d0, 1750.0d0,  434.1d0   , 0.452d0, & ! 27:H2O band
    &   1750.0d0, 1850.0d0,  136.0d0   , 0.359d0, & ! 28:H2O band
    &   1850.0d0, 1950.0d0,   35.7d0   , 0.165d0, & ! 29:H2O band
    &   1950.0d0, 2050.0d0,    9.015d0 , 0.104d0, & ! 30:H2O band
    &   2050.0d0, 2200.0d0,    1.529d0 , 0.116d0, & ! 31:H2O band
    !---
    &   2200.0d0, 2600.0d0,    0.0d0   , 0.0d0    / ! 32:nothing
!!$    &   2200.0d0, 2500.0d0,    0.0d0   , 0.0d0    / ! 32:nothing


!!$  data a_H2OBandParam &
!!$    & /    0.0d0,   40.0d0,    0.0d0   , 0.0d0  , & !  1:nothing
!!$    !---
!!$    &     40.0d0,  160.0d0, 7240.30d0  , 0.182d0, & !  2:H2O band
!!$    &    160.0d0,  280.0d0, 6024.80d0  , 0.094d0, & !  3:H2O band
!!$    &    280.0d0,  380.0d0, 1614.10d0  , 0.081d0, & !  4:H2O band
!!$    &    380.0d0,  500.0d0,  139.03d0  , 0.080d0, & !  5:H2O band
!!$    &    500.0d0,  600.0d0,   21.64d0  , 0.068d0, & !  6:H2O band
!!$    &    600.0d0,  720.0d0,    2.919d0 , 0.060d0, & !  7:H2O band
!!$    &    720.0d0,  800.0d0,    0.386d0 , 0.059d0, & !  8:H2O band
!!$    &    800.0d0,  900.0d0,    0.0715d0, 0.067d0, & !  9:H2O band and continuum
!!$    !---
!!$    &    900.0d0, 1000.0d0,    0.0d0   , 0.0d0  , & ! 10:H2O continuum
!!$    &   1000.0d0, 1200.0d0,    0.0d0   , 0.0d0  , & ! 11:H2O continuum
!!$    !--- 
!!$    &   1200.0d0, 1350.0d0,   12.65d0  , 0.089d0, & ! 12:H2O band
!!$    &   1350.0d0, 1450.0d0,  134.4d0   , 0.230d0, & ! 13:H2O band
!!$    &   1450.0d0, 1550.0d0,  632.0d0   , 0.320d0, & ! 14:H2O band
!!$    &   1550.0d0, 1650.0d0,  331.2d0   , 0.296d0, & ! 15:H2O band
!!$    &   1650.0d0, 1750.0d0,  434.1d0   , 0.452d0, & ! 16:H2O band
!!$    &   1750.0d0, 1850.0d0,  136.0d0   , 0.359d0, & ! 17:H2O band
!!$    &   1850.0d0, 1950.0d0,   35.7d0   , 0.165d0, & ! 18:H2O band
!!$    &   1950.0d0, 2050.0d0,    9.015d0 , 0.104d0, & ! 19:H2O band
!!$    &   2050.0d0, 2200.0d0,    1.529d0 , 0.116d0, & ! 20:H2O band
!!$    !---
!!$    &   2200.0d0, 9000.0d0,    0.0d0   , 0.0d0    / ! 21:nothing

!!$  data a_H2OContParam &
!!$    & /    0.0d0,    0.0d0,  0.0d0 , & !  1:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & !  2:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & !  3:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & !  4:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & !  5:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & !  6:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & !  7:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & !  8:
!!$    &    800.0d0,  900.0d0, 11.14d0, & !  9:
!!$    &    900.0d0, 1000.0d0,  7.36d0, & ! 10:
!!$    &   1000.0d0, 1200.0d0,  5.15d0, & ! 11:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 12:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 13:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 14:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 15:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 16:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 17:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 18:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 19:
!!$    &      0.0d0,    0.0d0,  0.0d0 , & ! 20:
!!$    &      0.0d0,    0.0d0,  0.0d0   / ! 21:


  data a_CO2BandParam &
    & /  582.0d0,  752.0d0,  718.7d0   , 0.448d0  / ! CO2 15 micron


!!$  data a_O3BandParam &
!!$    & / 1000.0d0, 1006.5d0,    3.21d-3 , 5.0d0,   & !  1:O3
!!$    &   1006.5d0, 1013.0d0,    6.42d-3 , 5.0d0,   & !  2:O3
!!$    &   1013.0d0, 1019.5d0,    1.28d-2 , 5.0d0,   & !  3:O3
!!$    &   1019.5d0, 1026.0d0,    2.14d-2 , 5.5d0,   & !  4:O3
!!$    &   1026.0d0, 1032.5d0,    2.35d-2 , 5.8d0,   & !  5:O3
!!$    &   1032.5d0, 1039.0d0,    1.71d-2 , 8.0d0,   & !  6:O3
!!$    &   1039.0d0, 1045.5d0,    1.18d-2 , 6.1d0,   & !  7:O3
!!$    &   1045.5d0, 1052.0d0,    2.78d-2 , 8.4d0,   & !  8:O3
!!$    &   1052.0d0, 1058.5d0,    3.53d-2 , 8.3d0,   & !  9:O3
!!$    &   1058.5d0, 1065.0d0,    1.28d-2 , 6.7d0,   & ! 10:O3
!!$    &   1065.0d0, 1200.0d0,    0.0d0   , 0.0d0    / ! 11: nothing


!!$  real(DP), save:: a_TMPFORO3BandParam(1:4,1:3)
!!$
!!$  data a_TMPFORO3BandParam &
!!$    & /    0.0d0, 1000.0d0,    0.0d0   , 0.0d0  , &
!!$    !---
!!$    &   1000.0d0, 1200.0d0,    0.0d0   , 0.0d0  , & ! O3
!!$    !---
!!$    &   1200.0d0, 9000.0d0,    0.0d0   , 0.0d0    /


  real(DP), save:: H2OScaleIndex
  real(DP), save:: CO2ScaleIndex
  real(DP), save:: O3ScaleIndex

  real(DP), save:: MeanMolWeight
  real(DP), save:: H2OMolWeight
  real(DP), save:: CO2MolWeight

  logical , save:: flag_save_time

  integer , save:: NGaussQuad

  real(DP), allocatable, save:: xyrra_TransSave  (:,:,:,:,:)

  type(DC_DIFFTIME), save:: IntTimeSave
                              ! 長波フラックスを計算する時間間隔.
                              ! Interval time of long wave flux calculation
  type(DC_DIFFTIME), 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(:,:)


  integer , parameter :: nRefAtm = 31
  real(DP), save:: RefAtmPro(1:5, 1:nRefAtm)


  ! Mid-Latitude Summer Atmosphere
  ! H       P         T       H2o        O3
  ! km      mbar       K       gm-3      gm-3
  data RefAtmPro &
    & / &
    &  0,   1.013e+03,   294.0,   1.40e+01,   6.0e-05, &
    &  1,   9.020e+02,   290.0,   9.30e+00,   6.0e-05, &
    &  2,   8.020e+02,   285.0,   5.85e+00,   6.0e-05, &
    &  3,   7.100e+02,   279.0,   3.43e+00,   6.2e-05, &
    &  4,   6.280e+02,   273.0,   1.89e+00,   6.4e-05, &
    &  5,   5.540e+02,   267.1,   1.00e+00,   6.6e-05, &
    &  6,   4.870e+02,   261.0,   6.09e-01,   6.9e-05, &
    &  7,   4.260e+02,   254.7,   3.71e-01,   7.5e-05, &
    &  8,   3.720e+02,   248.2,   2.10e-01,   7.9e-05, &
    &  9,   3.240e+02,   241.7,   1.18e-01,   8.6e-05, &
    & 10,   2.810e+02,   235.2,   6.43e-02,   9.0e-05, &
    & 11,   2.430e+02,   228.8,   2.19e-02,   1.1e-04, &
    & 12,   2.090e+02,   222.3,   6.46e-03,   1.2e-04, &
    & 13,   1.790e+02,   216.9,   1.66e-03,   1.5e-04, &
    & 14,   1.530e+02,   215.8,   9.95e-04,   1.8e-04, &
    & 15,   1.300e+02,   215.8,   8.40e-04,   1.9e-04, &
    & 16,   1.100e+02,   215.8,   7.10e-04,   2.1e-04, &
    & 17,   9.500e+01,   215.8,   6.14e-04,   2.4e-04, &
    & 18,   8.120e+01,   216.0,   5.24e-04,   2.8e-04, &
    & 19,   6.950e+01,   217.0,   4.46e-04,   3.2e-04, &
    & 20,   5.950e+01,   218.2,   3.80e-04,   3.4e-04, &
    & 21,   5.100e+01,   219.4,   3.24e-04,   3.6e-04, &
    & 22,   4.370e+01,   220.6,   2.76e-04,   3.6e-04, &
    & 23,   3.760e+01,   221.8,   2.36e-04,   3.4e-04, &
    & 24,   3.220e+01,   223.0,   2.01e-04,   3.2e-04, &
    & 25,   2.770e+01,   224.2,   1.72e-04,   3.0e-04, &
    & 30,   1.320e+01,   234.2,   7.85e-05,   2.0e-04, &
    & 35,   6.520e+00,   245.3,   3.70e-05,   9.2e-05, &
    & 40,   3.330e+00,   257.5,   1.80e-05,   4.1e-05, &
    & 45,   1.760e+00,   269.7,   9.09e-06,   1.3e-05, &
    & 50,   9.510e-01,   276.2,   4.80e-06,   4.3e-06  &
    & /

  real(DP), save :: VMRCO2            ! Volume mixing ratio of CO2

  logical , save:: sw_fs

  data sw_fs /.true./


contains

  subroutine RadiationRL78Flux( &
    & xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, xy_SurfTemp, & ! (in )
    & xyr_RadLFlux, xyra_DelRadLFlux                                  & ! (out)
    & )


    ! 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

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav, & ! $ g $ [m s-2].
                                 ! 
                                 ! Gravitational acceleration
                         PI      ! $ \pi $ .
                                 ! Circular constant

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻.
                              ! Time of step $ t $.
      & EndTime, &            ! 計算終了時刻.
                              ! End time of calculation
      & TimesetClockStart, TimesetClockStop


    ! 日付および時刻の取り扱い
    ! Date and time handler
    !
    use dc_date, only: operator(-), operator(>=), operator(+), &
      & operator(==), toChar

    use read_time_series, only: SetValuesFromTimeSeriesWrapper

    use set_cloud, only : SetCloudLW

    ! OpenMP
    !
    !$ use omp_lib

    real(DP), intent(in ):: xyz_Temp        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QVap        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QO3         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyr_Press       (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ):: xyz_Press       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xy_SurfTemp     (0:imax-1, 1:jmax)
    real(DP), intent(out):: xyr_RadLFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


    ! 
    ! Work variables
    !
    real(DP):: xyz_QCO2              (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_CloudDelOptDep    (0:imax-1, 1:jmax, 1:kmax)

    logical :: flag_dry_atmosphere

    logical :: a_TransFlag(1:nbmax)


    integer :: n

    integer :: nt

    integer :: js
    integer :: je

    integer :: nthreads
    integer, allocatable :: a_js(:)
    integer, allocatable :: a_je(:)



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


    if ( sw_fs ) then
      call RadiationRoeweLiouInit
    end if


    nthreads = 1
    !$ nthreads  = omp_get_max_threads()

    allocate( a_js(0:nthreads-1) )
    allocate( a_je(0:nthreads-1) )

    do n = 0, nthreads-1

      if ( n == 0 ) then
        a_js(n) = 1
      else
        a_js(n) = a_je(n-1) + 1
      end if

      a_je(n) = a_js(n  ) + jmax / nthreads - 1
      if ( n + 1 <= mod( jmax, nthreads ) ) then
        a_je(n) = a_je(n) + 1
      end if

    end do



!!$    xyz_QCO2(:,:,:) = 382.0d-6 * CO2MolWeight / MeanMolWeight
    xyz_QCO2(:,:,:) = VMRCO2 * CO2MolWeight / MeanMolWeight
!!$    xyz_QO3 (:,:,:) = 0.0d0


!!$    xyz_QVap(:,:,:) = 0.0d0
!!$    xyz_QCO2(:,:,:) = 300.0d-6 * CO2MolWeight / MeanMolWeight
!!$    xyz_QCO2(:,:,:) = 0.0d0
!!$    call SetRefAtmPro(      &
!!$      & 5, .true., xyz_Press, &
!!$      & xyz_QO3                 &
!!$      & )


!!$    if ( .true. ) then
    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.


      ! Check for dry atmosphere
      !
      flag_dry_atmosphere = .false.
      if ( flag_save_time ) then
        if ( all( xyz_QVap <= 0.0d0 ) ) then
          flag_dry_atmosphere = .true.
!!$          write( 6, * ) 'Dry atmosphere'
        else
          flag_dry_atmosphere = .false.
        end if
      end if


      call SetCloudLW(    &
        & xyz_CloudDelOptDep    &
        & )




      !$OMP PARALLEL DEFAULT(PRIVATE) &
      !$OMP SHARED(nthreads,a_js,a_je, &
      !$OMP        xyz_QCO2, &
      !$OMP        xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, xy_SurfTemp,  &
      !$OMP        xyz_CloudDelOptDep, &
      !$OMP        flag_dry_atmosphere)
      !$OMP DO

      LOOP_thread_trans : do nt = 0, nthreads-1

        js = a_js(nt)
        je = a_je(nt)

        if ( js > je ) cycle

!!$      write( 6, * ) n, js, je


        call RadiationRL78CalcTotTrans( &
          & xyz_QCO2, &
          & xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, & ! (in )
          & xyz_CloudDelOptDep, &
          & flag_dry_atmosphere, &
          & js, je  &
          & )

      end do LOOP_thread_trans

      !$OMP END DO
      !$OMP END PARALLEL


    end if



    a_TransFlag = .false.
    if ( flag_save_time ) then
      do n = 1, nbmax
        if ( all( xyrra_TransSave(:,:,0,kmax,n) == 1.0d0 ) ) then
          a_TransFlag(n) = .true.
        else
          a_TransFlag(n) = .false.
        end if
      end do
    end if


    !$OMP PARALLEL DEFAULT(PRIVATE) &
    !$OMP SHARED(nthreads,a_js,a_je, &
    !$OMP        a_TransFlag,                       &
    !$OMP        xyz_Temp, xy_SurfTemp,             &
    !$OMP        xyr_RadLFlux, xyra_DelRadLFlux)
    !$OMP DO

    LOOP_thread_RTE : do nt = 0, nthreads-1

      js = a_js(nt)
      je = a_je(nt)

      if ( js > je ) cycle


      call RadiationRL78IntegRTE( &
        & a_TransFlag,                                                    & ! (in )
        & xyz_Temp, xy_SurfTemp,                                          & ! (in )
        & xyr_RadLFlux, xyra_DelRadLFlux,                                 & ! (out)
        & js, je                                                          & ! (in )
        & )

    end do LOOP_thread_RTE

    !$OMP END DO
    !$OMP END PARALLEL


    deallocate( a_js )
    deallocate( a_je )


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


  end subroutine RadiationRL78Flux

  !**************************************************************************

  subroutine RadiationRL78IntegRTE( &
    & a_TransFlag,                                                    & ! (in )
    & xyz_Temp, xy_SurfTemp,                                          & ! (in )
    & xyr_RadLFlux, xyra_DelRadLFlux,                                 & ! (out)
    & js, je                                                          & ! (in )
    & )


    ! 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

    ! 
    ! Physical constants settings
    !
    use constants, only: PI      ! $ \pi $ .
                                 ! Circular constant

    logical , intent(in ):: a_TransFlag     (1:nbmax)
    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(out):: xyr_RadLFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    integer , intent(in ):: js
    integer , intent(in ):: je



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

    real(DP):: xyz_PlankFunc         (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xy_SurfPlankFunc      (0:imax-1, 1:jmax)
    real(DP):: xy_DPFDT0             (0:imax-1, 1:jmax)
    real(DP):: xy_DPFDT1             (0:imax-1, 1:jmax)


    integer :: j
    integer :: k, kk, k1, k2
    integer :: l
    integer :: n


    do k = 0, kmax
      do j = js, je
        xyr_RadLFlux(:,j,k)   = 0.0d0
      end do
    end do
    do l = 0, 1
      do k = 0, kmax
        do j = js, je
          xyra_DelRadLFlux(:,j,k,l) = 0.0d0
        end do
      end do
    end do


    LOOP_band: do n = 1, nbmax

      do k2 = 0, kmax
        do k1 = 0, kmax
          do j = js, je
            xyrr_Trans(:,j,k1,k2) = xyrra_TransSave(:,j,k1,k2,n)
          end do
        end do
      end do


      if ( a_TransFlag(n) ) then

        call CalcIntegratedPFWithTable2D( &
          & n, xy_SurfTemp,                         &
          & xy_SurfPlankFunc,                       &
          & js, je                                  &
          & )

        do k = 0, kmax
          do j = js, je
            xyr_RadLFlux(:,j,k) = xyr_RadLFlux(:,j,k) &
              & + PI * xy_SurfPlankFunc(:,j) ! * xyrr_Trans(:,:,0,k)
          end do
        end do

        call CalcIntegratedPFWithTable2D( &
          & n, xy_SurfTemp,                         &
          & xy_DPFDT0,                              &
          & js, je,                                 &
          & .true.                                  &
          & )

        do k = 0, kmax
          do j = js, je
            xyra_DelRadLFlux(:,j,k,0) = xyra_DelRadLFlux(:,j,k,0) &
              & + PI * xy_DPFDT0(:,j) ! * xyrr_Trans(:,:,0,k)
          end do
        end do

      else

        call CalcIntegratedPFWithTable2D( &
          & n, xy_SurfTemp,                         &
          & xy_SurfPlankFunc,                       &
          & js, je                                  &
          & )
        call CalcIntegratedPFWithTable3D( &
          & n, kmax, xyz_Temp,                    &
          & xyz_PlankFunc,                        &
          & js, je                                &
          & )

        do k = 0, kmax
          do j = js, je
            xyr_RadLFlux(:,j,k) = xyr_RadLFlux(:,j,k) &
              & + PI * xy_SurfPlankFunc(:,j) * xyrr_Trans(:,j,0,k)
          end do
          do kk = 1, kmax
            do j = js, je
              xyr_RadLFlux(:,j,k) = xyr_RadLFlux(:,j,k) &
                & - PI * xyz_PlankFunc(:,j,kk) &
                & * ( xyrr_Trans(:,j,k,kk-1) - xyrr_Trans(:,j,k,kk) )
            end do
          end do
        end do

        call CalcIntegratedPFWithTable2D( &
          & n, xy_SurfTemp,                         &
          & xy_DPFDT0,                              &
          & js, je,                                 &
          & .true.                                  &
          & )
        call CalcIntegratedPFWithTable2D( &
          & n, xyz_Temp(:,:,1),                     &
          & xy_DPFDT1,                              &
          & js, je,                                 &
          & .true.                                  &
          & )

        do k = 0, kmax
          do j = js, je
            xyra_DelRadLFlux(:,j,k,0) = xyra_DelRadLFlux(:,j,k,0) &
              & + PI * xy_DPFDT0(:,j) * xyrr_Trans(:,j,0,k)
            xyra_DelRadLFlux(:,j,k,1) = xyra_DelRadLFlux(:,j,k,1) &
              & - PI * xy_DPFDT1(:,j) &
              & * ( xyrr_Trans(:,j,k,0) - xyrr_Trans(:,j,k,1) )
          end do
        end do

      end if


    end do LOOP_band


  end subroutine RadiationRL78IntegRTE

  !**************************************************************************

  subroutine RadiationRL78CalcTotTrans( &
    & xyz_QCO2, &
    & xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, & ! (in )
    & xyz_CloudDelOptDep, &
    & flag_dry_atmosphere, &
    & js, je  &
    & )


    ! 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

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav    ! $ g $ [m s-2].
                                 ! 
                                 ! Gravitational acceleration

    real(DP), intent(in ):: xyz_QCO2          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_Temp          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QVap          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QO3           (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyr_Press         (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ):: xyz_Press         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudDelOptDep(0:imax-1, 1:jmax, 1:kmax)
    logical , intent(in ):: flag_dry_atmosphere
    integer , intent(in ):: js
    integer , intent(in ):: je


    ! 
    ! Work variables
    !
    real(DP):: RefPress

    real(DP):: xyz_H2ODelAbsAmt      (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_H2ODelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_O3DelScaledAbsAmt (0:imax-1, 1:jmax, 1:kmax)

    real(DP):: xyrr_TransCO2         (0:imax-1, 1:jmax, 0:kmax, 0: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)

    integer :: j, k, k1, k2
    integer :: n



    if ( sw_fs ) then
      call RadiationRoeweLiouInit
    end if


    do k = 1, kmax
      do j = js, je
        xyz_H2ODelAbsAmt(:,j,k) = &
          &   xyz_QVap(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do

    RefPress = 1.0d5

    do k = 1, kmax
      do j = js, je
        xyz_H2ODelScaledAbsAmt(:,j,k) = &
          &   ( xyz_Press(:,j,k) / RefPress )**H2OScaleIndex &
          & * xyz_QVap(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do
    do k = 1, kmax
      do j = js, je
        xyz_CO2DelScaledAbsAmt(:,j,k) = &
          &   ( xyz_Press(:,j,k) / RefPress )**CO2ScaleIndex &
          & * xyz_QCO2(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do
    do k = 1, kmax
      do j = js, je
        xyz_O3DelScaledAbsAmt(:,j,k) = &
          &   ( xyz_Press(:,j,k) / RefPress )**O3ScaleIndex &
          & * xyz_QO3(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do



    !-----

    ! Transmission in CO2 15 micron band
    !
    call RadiationRoeweLiouCalcCO2Trans( &
      & xyz_CO2DelScaledAbsAmt, & ! (in )
      & xyrr_TransCO2,        & ! (out)
      & js, je                  & ! (in )
      & )


    !-----

    ! Transmission by cloud

!!$      if ( all( xyz_CloudDelOptDep <= 0.0d0 ) ) then
!!$
!!$        xyrr_TransCloud(:,:,:,:) = 1.0d0
!!$
!!$      else

    do k = 1, kmax
      do j = js, je
        xyz_TransCloudOneLayer(:,j,k) = exp( - xyz_CloudDelOptDep(:,j,k) * DiffFactor )
      end do
    end do

    do k1 = 0, kmax
      k2 = k1
      do j = js, je
        xyrr_TransCloud(:,j,k1,k2) = 1.0d0
      end do
      do k2 = k1+1, kmax
        do j = js, je
          xyrr_TransCloud(:,j,k1,k2) = xyrr_TransCloud(:,j,k1,k2-1) &
            & * xyz_TransCloudOneLayer(:,j,k2)
        end do
      end do
    end do
    do k1 = 0, kmax
      do k2 = 0, k1-1
        do j = js, je
          xyrr_TransCloud(:,j,k1,k2) = xyrr_TransCloud(:,j,k2,k1)
        end do
      end do
    end do

!!$      end if

    !-----

    do n = 1, nbmax

      if ( n <= 5 ) then
        ! H2O band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( &
            & n, xyz_H2ODelScaledAbsAmt, & ! (in )
            & xyrr_Trans,            & ! (out)
            & js, je                      & ! (in )
            & )
        end if
      else if ( n <= 8 ) then
        ! H2O band & CO2 band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( &
            & n, xyz_H2ODelScaledAbsAmt, & ! (in )
            & xyrr_Trans,            & ! (out)
            & js, je                      & ! (in )
            & )
        end if
        do k2 = 0, kmax
          do k1 = 0, kmax
            do j = js, je
              xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k1,k2) * xyrr_TransCO2(:,j,k1,k2)
            end do
          end do
        end do
      else if ( n == 9 ) then
        ! H2O band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( &
            & n, xyz_H2ODelScaledAbsAmt, & ! (in )
            & xyrr_Trans,            & ! (out)
            & js, je                      & ! (in )
            & )
        end if
      else if ( n == 10 ) then
        ! H2O band & continuum
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( &
            & n, xyz_H2ODelScaledAbsAmt, & ! (in )
            & xyrr_Trans,            & ! (out)
            & js, je                      & ! (in )
            & )
          call RadiationRLCalcH2OContTrans( &
            & n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, & ! (in   )
            & xyrr_Trans                                           & ! (inout)
            & )
        end if
      else if ( n == 11 ) then
        ! H2O continuum
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
          call RadiationRLCalcH2OContTrans( &
            & n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, & ! (in   )
            & xyrr_Trans                                           & ! (inout)
            & )
        end if
      else if ( n <= 21 ) then
        ! O3 band & H2O continuum
        call RadiationRoeweLiouCalcBandTrans( &
          & n, xyz_O3DelScaledAbsAmt, & ! (in )
          & xyrr_Trans,            & ! (out)
          & js, je                      & ! (in )
          & )
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k1,k2)
              end do
            end do
          end do
        else
          call RadiationRLCalcH2OContTrans( &
            & n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, & ! (in   )
            & xyrr_Trans                                           & ! (inout)
            & )
        end if
      else if ( n == 22 ) then
        ! H2O continuum
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
          call RadiationRLCalcH2OContTrans( &
            & n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, & ! (in   )
            & xyrr_Trans                                           & ! (inout)
            & )
        end if
      else if ( n <= 31 ) then
        ! H2O band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( &
            & n, xyz_H2ODelScaledAbsAmt, & ! (in )
            & xyrr_Trans,            & ! (out)
            & js, je                      & ! (in )
            & )
        end if
      else if ( n == 32 ) then
        do k2 = 0, kmax
          do k1 = 0, kmax
            do j = js, je
              xyrr_Trans(:,j,k1,k2) = 1.0d0
            end do
          end do
        end do
      else
        write( 6, * ) 'Unexpected number of band, ', n, '.'
        stop
      end if

      do k2 = 0, kmax
        do k1 = 0, kmax
          do j = js, je
            xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k1,k2) * xyrr_TransCloud(:,j,k1,k2)

            xyrra_TransSave(:,j,k1,k2,n) = xyrr_Trans(:,j,k1,k2)
          end do
        end do
      end do

    end do


  end subroutine RadiationRL78CalcTotTrans

  !**************************************************************************

  subroutine RadiationRoeweLiouCalcBandTrans( &
    & iband, xyz_DelScaledAbsAmt, & ! (in )
    & xyrr_Trans,                 & ! (out)
    & js, je                      & ! (in )
    & )

    ! 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

    integer , intent(in ):: iband
    real(DP), intent(in ):: xyz_DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyrr_Trans         (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    integer , intent(in ):: js
    integer , intent(in ):: je

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

    integer :: j
    integer :: k1
    integer :: k2
    integer :: k3
    integer :: n


    n = iband


    do k2 = 0, kmax
      do k1 = 0, kmax
        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = 0.0d0
        end do
      end do
    end do

    do k1 = 0, kmax
      do k2 = k1, kmax

        do k3 = k1+1, k2
          do j = js, je
            xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) &
              & + xyz_DelScaledAbsAmt(:,j,k3)
          end do
        end do

        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) * DiffFactor
        end do

        do j = js, je
          xyrr_Trans(:,j,k1,k2) = &
            & - a_BandParam(3,n) * xyrr_AbsAmt(:,j,k1,k2) &
            & / sqrt( 1.0d0 + a_BandParam(3,n) / a_BandParam(4,n) &
            &               * xyrr_AbsAmt(:,j,k1,k2) )
          xyrr_Trans(:,j,k1,k2) = exp( xyrr_Trans(:,j,k1,k2) )
        end do

      end do
    end do


    do k1 = 0, kmax
      do k2 = 0, k1-1
        do j = js, je
          xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k2,k1)
        end do
      end do
    end do
    do k1 = 0, kmax
      k2 = k1
      do j = js, je
        xyrr_Trans(:,j,k1,k2) = 1.0d0
      end do
    end do


  end subroutine RadiationRoeweLiouCalcBandTrans

  !**************************************************************************
!!$
!!$  subroutine RadiationRoeweLiouCalcH2OTrans( &
!!$    & iband, xyz_DelScaledAbsAmt, & ! (in )
!!$    & xyrr_Trans            & ! (out)
!!$    & )
!!$
!!$    ! 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
!!$
!!$    integer , intent(in ):: iband
!!$    real(DP), intent(in ):: xyz_DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP), intent(out):: xyrr_Trans         (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$
!!$
!!$    ! 
!!$    ! Work variables
!!$    !
!!$    real(DP):: xyrr_AbsAmt       (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$
!!$    integer :: k1
!!$    integer :: k2
!!$    integer :: k3
!!$    integer :: n
!!$
!!$
!!$    n = iband
!!$
!!$
!!$    xyrr_AbsAmt(:,:,:,:) = 0.0d0
!!$
!!$    do k1 = 0, kmax
!!$      do k2 = k1, kmax
!!$
!!$        do k3 = k1+1, k2
!!$          xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) &
!!$            & + xyz_DelScaledAbsAmt(:,:,k3)
!!$        end do
!!$
!!$        xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) * DiffFactor
!!$
!!$        xyrr_Trans(:,:,k1,k2) = &
!!$          & - a_H2OBandParam(3,n) * xyrr_AbsAmt(:,:,k1,k2) &
!!$          & / sqrt( 1.0d0 + a_H2OBandParam(3,n) / a_H2OBandParam(4,n) &
!!$          &               * xyrr_AbsAmt(:,:,k1,k2) )
!!$        xyrr_Trans(:,:,k1,k2) = exp( xyrr_Trans(:,:,k1,k2) )
!!$
!!$      end do
!!$    end do
!!$
!!$
!!$    do k1 = 0, kmax
!!$      do k2 = 0, k1-1
!!$        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
!!$      end do
!!$    end do
!!$    do k1 = 0, kmax
!!$      k2 = k1
!!$      xyrr_Trans(:,:,k1,k2) = 1.0d0
!!$    end do
!!$
!!$
!!$  end subroutine RadiationRoeweLiouCalcH2OTrans
!!$
  !**************************************************************************

  subroutine RadiationRLCalcH2OContTrans( &
    & iband, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, & ! (in   )
    & xyrr_Trans                                              & ! (inout)
    & )

    ! 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

    integer , intent(in   ):: iband
    real(DP), intent(in   ):: xyz_Temp        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ):: xyz_QVap        (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_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xyrr_Trans      (0:imax-1, 1:jmax, 0:kmax, 0:kmax)


    ! 
    ! Work variables
    !
    real(DP):: H2OContConstant
    real(DP):: a_H2OContParam(3)
    real(DP):: RefTemp
    real(DP):: gamma
    real(DP):: xyz_PressH2O      (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DelOptDep     (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyrr_H2OContOptDep(0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    integer :: k1
    integer :: k2
    integer :: k3
    integer :: n


    if ( ( iband < 9 ) .and. ( iband > 11 ) ) then
      write( 6, * ) 'Unexpected number of band in RadiationRLCalcH2OContTrans.'
      stop
    end if


    ! These constants are obtained from Roberts et al. (1976) [Applied Optics]
    a_H2OContParam(1) = 1.25d-22 ! mol-1 cm2 atm-1
    a_H2OContParam(2) = 1.67d-19 ! mol-1 cm2 atm-1
    a_H2OContParam(3) = 7.87d-3  ! (cm)

    ! unit conversion
    a_H2OContParam(1) = a_H2OContParam(1) &
      & / ( H2OMolWeight / AvogNum ) * 1.0d-4 / 101325.0d0 ! (kg-1 m2 Pa-1)
    a_H2OContParam(2) = a_H2OContParam(2) &
      & / ( H2OMolWeight / AvogNum ) * 1.0d-4 / 101325.0d0 ! (kg-1 m2 Pa-1)
    a_H2OContParam(3) = a_H2OContParam(3) * 1.0d-2    ! (m)


    n = iband


    ! This constant is deribed from Roberts et al. (1976) [Applied Optics]
!!$    gamma = 0.005d0
    gamma = 0.0d0

    xyz_PressH2O = xyz_Press * xyz_QVap * MeanMolWeight / H2OMolWeight

    RefTemp = 1800.0d0

    H2OContConstant =                                                &
      & a_H2OContParam(1)                                           &
      & + a_H2OContParam(2)                                         &
      & * exp( - a_H2OContParam(3)                                  &
      &          * ( a_BandParam(1,n) + a_BandParam(2,n) ) * 0.5d0 )

    xyz_DelOptDep = H2OContConstant                               &
      & * exp( RefTemp * ( 1.0d0 / xyz_Temp - 1.0d0 / 296.0d0 ) ) &
      & * ( xyz_PressH2O + gamma * ( xyz_Press - xyz_PressH2O ) ) &
      & * xyz_H2ODelAbsAmt

    xyrr_H2OContOptDep(:,:,:,:) = 0.0d0
    do k1 = 0, kmax
      do k2 = k1, kmax

        do k3 = k1+1, k2
          xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) &
            & + xyz_DelOptDep(:,:,k3)
        end do

        xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) &
          & * DiffFactor

        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k1,k2) &
          & * exp( - xyrr_H2OContOptDep(:,:,k1,k2) )

      end do
    end do

    do k1 = 0, kmax
      do k2 = 0, k1-1
        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
      end do
    end do
    do k1 = 0, kmax
      k2 = k1
      xyrr_Trans(:,:,k1,k2) = 1.0d0
    end do


  end subroutine RadiationRLCalcH2OContTrans

  !**************************************************************************
!!$
!!$  subroutine RadiationRoeweLiouCalcH2OContTrans0( &
!!$    & iband, xyz_Temp, xyz_QVap, xyz_Press, xyz_DelAbsAmt, & ! (in   )
!!$    & xyrr_Trans                                           & ! (inout)
!!$    & )
!!$
!!$    ! 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
!!$
!!$    integer , intent(in   ):: iband
!!$    real(DP), intent(in   ):: xyz_Temp     (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP), intent(in   ):: xyz_QVap     (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_DelAbsAmt(0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP), intent(inout):: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$
!!$
!!$    ! 
!!$    ! Work variables
!!$    !
!!$    real(DP):: RefTemp
!!$    real(DP):: gamma
!!$    real(DP):: xyz_PressH2O      (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP):: xyz_DelOptDep     (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP):: xyrr_H2OContOptDep(0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$
!!$    integer :: k1
!!$    integer :: k2
!!$    integer :: k3
!!$    integer :: n
!!$
!!$
!!$    if ( ( iband < 9 ) .and. ( iband > 11 ) ) then
!!$      write( 6, * ) 'Unexpected number of band in RadiationRoeweLiouCalcH2OContTrans.'
!!$      stop
!!$    end if
!!$
!!$    n = iband
!!$
!!$
!!$    ! This constant is deribed from Roberts et al. (1976) [Applied Optics]
!!$    gamma = 0.005d0
!!$
!!$    xyz_PressH2O = xyz_Press * xyz_QVap * MeanMolWeight / H2OMolWeight
!!$
!!$    RefTemp = 1800.0d0
!!$    xyz_DelOptDep = a_H2OContParam(3,n)                           &
!!$      & * exp( RefTemp * ( 1.0d0 / xyz_Temp - 1.0d0 / 296.0d0 ) ) &
!!$      & * ( xyz_PressH2O + gamma * ( xyz_Press - xyz_PressH2O ) ) &
!!$      & * xyz_DelAbsAmt
!!$
!!$    xyrr_H2OContOptDep(:,:,:,:) = 0.0d0
!!$    do k1 = 0, kmax
!!$      do k2 = k1, kmax
!!$
!!$        do k3 = k1+1, k2
!!$          xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) &
!!$            & + xyz_DelOptDep(:,:,k3)
!!$        end do
!!$
!!$        xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) &
!!$          & * DiffFactor
!!$
!!$        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k1,k2) &
!!$          & * exp( - xyrr_H2OContOptDep(:,:,k1,k2) )
!!$
!!$      end do
!!$    end do
!!$
!!$    do k1 = 0, kmax
!!$      do k2 = 0, k1-1
!!$        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
!!$      end do
!!$    end do
!!$    do k1 = 0, kmax
!!$      k2 = k1
!!$      xyrr_Trans(:,:,k1,k2) = 1.0d0
!!$    end do
!!$
!!$
!!$  end subroutine RadiationRoeweLiouCalcH2OContTrans0
!!$
  !**************************************************************************

  subroutine RadiationRoeweLiouCalcCO2Trans( &
!!$    & iband, xyz_CO2DelScaledAbsAmt, & ! (in )
    & xyz_CO2DelScaledAbsAmt, & ! (in )
    & xyrr_Trans,             & ! (out)
    & js, je                  & ! (in )
    & )

    ! 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

!!$    integer , intent(in ):: iband
    real(DP), intent(in ):: xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyrr_Trans            (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    integer , intent(in ):: js
    integer , intent(in ):: je


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

    integer :: j
    integer :: k1
    integer :: k2
    integer :: k3
!!$    integer :: n


!!$    n = iband
!!$    n = 2


    do k2 = 0, kmax
      do k1 = 0, kmax
        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = 0.0d0
        end do
      end do
    end do

    do k1 = 0, kmax
      do k2 = k1, kmax

        do k3 = k1+1, k2
          do j = js, je
            xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) &
              & + xyz_CO2DelScaledAbsAmt(:,j,k3)
          end do
        end do

        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) * DiffFactor
        end do

        do j = js, je
          xyrr_Trans(:,j,k1,k2) = &
            & - a_CO2BandParam(3) * xyrr_AbsAmt(:,j,k1,k2) &
            & / sqrt( 1.0d0 + a_CO2BandParam(3) / a_CO2BandParam(4) &
            &               * xyrr_AbsAmt(:,j,k1,k2) )
          xyrr_Trans(:,j,k1,k2) = exp( xyrr_Trans(:,j,k1,k2) )
        end do

      end do
    end do


    do k1 = 0, kmax
      do k2 = 0, k1-1
        do j = js, je
          xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k2,k1)
        end do
      end do
    end do
    do k1 = 0, kmax
      k2 = k1
      do j = js, je
        xyrr_Trans(:,j,k1,k2) = 1.0d0
      end do
    end do


  end subroutine RadiationRoeweLiouCalcCO2Trans

  !**************************************************************************
!!$
!!$  subroutine RadiationRoeweLiouCalcO3Trans( &
!!$    & iband, xyz_DelScaledAbsAmt, & ! (in )
!!$    & xyrr_Trans            & ! (inout)
!!$    & )
!!$
!!$    ! 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
!!$
!!$    integer , intent(in   ):: iband
!!$    real(DP), intent(in   ):: xyz_DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP), intent(inout):: xyrr_Trans         (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$
!!$
!!$    ! 
!!$    ! Work variables
!!$    !
!!$    real(DP):: xyrr_TransTentative(0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$    real(DP):: xyrr_AbsAmt        (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$    real(DP):: xy_OpDep           (0:imax-1, 1:jmax)
!!$
!!$    integer :: k1
!!$    integer :: k2
!!$    integer :: k3
!!$    integer :: n
!!$
!!$
!!$!    stop 'Now cannot use'
!!$
!!$
!!$    xyrr_AbsAmt(:,:,:,:) = 0.0d0
!!$
!!$    do k1 = 0, kmax
!!$      do k2 = k1, kmax
!!$
!!$        do k3 = k1+1, k2
!!$          xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) &
!!$            & + xyz_DelScaledAbsAmt(:,:,k3)
!!$        end do
!!$
!!$        xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) * DiffFactor
!!$
!!$      end do
!!$    end do
!!$
!!$
!!$    xyrr_TransTentative = 0.0d0
!!$
!!$    do n = 1, nbO3max-1
!!$      do k1 = 0, kmax
!!$        do k2 = k1, kmax
!!$          xy_OpDep(:,:) = &
!!$            & a_O3BandParam(3,n) * xyrr_AbsAmt(:,:,k1,k2)             &
!!$            & / sqrt( 1.0d0 + a_O3BandParam(3,n) / a_O3BandParam(4,n) &
!!$            &         * xyrr_AbsAmt(:,:,k1,k2) )
!!$
!!$          xyrr_TransTentative(:,:,k1,k2) = xyrr_TransTentative(:,:,k1,k2) &
!!$            & + exp( - xy_OpDep(:,:) )                                    &
!!$            & * ( a_O3BandParam(2,n      ) - a_O3BandParam(1,n) )         &
!!$            & / ( a_O3BandParam(2,nbO3max) - a_O3BandParam(1,1) )
!!$        end do
!!$      end do
!!$    end do
!!$
!!$    n = nbO3max
!!$    do k1 = 0, kmax
!!$      do k2 = k1, kmax
!!$        xyrr_TransTentative(:,:,k1,k2) = xyrr_TransTentative(:,:,k1,k2) &
!!$          & + 1.0d0                                                     &
!!$          & * ( a_O3BandParam(2,n      ) - a_O3BandParam(1,n) )         &
!!$          & / ( a_O3BandParam(2,nbO3max) - a_O3BandParam(1,1) )
!!$      end do
!!$    end do
!!$
!!$    do k1 = 0, kmax
!!$      do k2 = k1, kmax
!!$        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k1,k2) &
!!$          & * xyrr_TransTentative(:,:,k1,k2)
!!$      end do
!!$    end do
!!$
!!$    do k1 = 0, kmax
!!$      do k2 = 0, k1-1
!!$        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
!!$      end do
!!$    end do
!!$    do k1 = 0, kmax
!!$      k2 = k1
!!$      xyrr_Trans(:,:,k1,k2) = 1.0d0
!!$    end do
!!$
!!$
!!$  end subroutine RadiationRoeweLiouCalcO3Trans
!!$
  !**************************************************************************

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

    ! USE statements
    !

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

    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
    !

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

    ! メッセージ出力
    ! 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
    integer                     :: n


    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

  !**************************************************************************
!!$
!!$  subroutine CalcIntegratedPF( iband, NGaussQuad, GQP, GQW, Temp, IPF )
!!$
!!$
!!$    integer , intent(in ) :: iband
!!$    integer , intent(in ) :: NGaussQuad
!!$    real(DP), intent(in ) :: GQP(1:NGaussQuad)
!!$    real(DP), intent(in ) :: GQW(1:NGaussQuad)
!!$    real(DP), intent(in ) :: Temp
!!$    real(DP), intent(out) :: IPF
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer :: l
!!$
!!$
!!$    IPF = 0.0d0
!!$
!!$    do l = 1, NGaussQuad
!!$      IPF = IPF + pf( GQP(l), Temp ) * GQW( l )
!!$    end do
!!$
!!$
!!$  end subroutine CalcIntegratedPF
!!$
  !**************************************************************************

  subroutine CalcIntegratedPF( wn1, wn2, num, im, jm, km, temp, pfinted )


    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: im, jm, km
    real(DP), intent(in ) :: temp   ( im, jm, km )
    real(DP), intent(out) :: pfinted( im, jm, km )


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


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

    pfinted(:,:,:) = 0.0d0

    do l = 1, num
      do k = 1, km
        do j = 1, jm
          do i = 1, im
            pfinted(i,j,k) = pfinted(i,j,k) &
              & + pf( x( l ), temp(i,j,k) ) * w( l )
          end do
        end do
      end do
    end do


  end subroutine CalcIntegratedPF

  !**************************************************************************

  subroutine CalcIntegratedPF2D( wn1, wn2, num, im, jm, temp, pfinted )


    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: im, jm
    real(DP), intent(in ) :: temp   ( im, jm )
    real(DP), intent(out) :: pfinted( im, jm )


    !
    ! local variables
    !
    real(dp)     :: temp3d( im, jm, 1 ), pfinted3d( im, jm, 1 )


    temp3d(:,:,1) = temp(:,:)
    call CalcIntegratedPF( wn1, wn2, num, im, jm, 1, temp3d, pfinted3d )
    pfinted(:,:) = pfinted3d(:,:,1)


  end subroutine CalcIntegratedPF2D

  !**************************************************************************
  ! function pf
  !
  ! Calculate Planc function at wanumber x and temperature t
  !
  ! x         : Wavenumber (m^{-1})
  ! t         : Temperature (K)
  !**************************************************************************

  real(DP) function pf(x,t)

    real(DP), intent(in) :: x, t


    real(dp), parameter :: &
      & c      = 2.99792458d8           , &
      & planc  = 6.6260755d-34          , &
      & boltz  = 1.380658d-23


    pf=2.0d0*planc*c*c*x*x*x &
      /(exp(planc*c*(x+1.0d-10)/(boltz*t))-1.0d0)

    return
  end function pf

  !**************************************************************************

  real(DP) function DPFDT( &
    & WN,    & ! (in )
    & Temp   & ! (in )
    & )

    ! USE statements
    !

    real(DP), intent(in ) :: WN
    real(DP), intent(in ) :: Temp


    real(DP), parameter :: &
      & c      = 2.99792458d8           , &
      & planc  = 6.6260755d-34          , &
      & boltz  = 1.380658d-23
    real(DP) :: ExpTerm
    real(DP) :: PF


    ExpTerm = exp( planc * c * ( WN + 1.0d-10 ) / ( boltz * Temp ) )

    PF = 2.0d0 * planc * c * c * WN * WN * WN &
      / ( ExpTerm - 1.0d0 )

    DPFDT = &
      & 1.0d0 / ( 2.0d0 * c * WN * WN * boltz ) &
      & * ( pf / Temp )**2 &
      & * ExpTerm


  end function DPFDT

  !**************************************************************************

  subroutine Integ_DPFDT_GQ_Array2D( &
    & WN1, WN2, Num, xy_Temp, & ! (in )
    & xy_DPFDTInted           & ! (out)
    & )

    ! 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 ) :: WN1
    real(DP), intent(in ) :: WN2
    integer , intent(in ) :: Num
    real(DP), intent(in ) :: xy_Temp      (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xy_DPFDTInted(0:imax-1, 1:jmax)


    !
    ! local variables
    !
    real(DP):: GP( Num )
    real(DP):: GW( Num )
    real(DP):: xy_DPFDT(0:imax-1, 1:jmax)
    integer :: l


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

    xy_DPFDTInted(:,:) = 0.0d0

    do l = 1, num
      call DPFDT_Array2D( &
        & GP(l),       & ! (in )
        & xy_Temp,  & ! (in )
        & xy_DPFDT  & ! (out)
        & )

      xy_DPFDTInted(:,:) = xy_DPFDTInted(:,:) + xy_DPFDT(:,:) * GW(l)
    end do


  end subroutine Integ_DPFDT_GQ_Array2D

  !**************************************************************************

  subroutine DPFDT_Array2D( &
    & WN,       & ! (in )
    & xy_Temp,  & ! (in )
    & xy_DPFDT  & ! (out)
    & )

    ! 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 ) :: WN
    real(DP), intent(in ) :: xy_Temp (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xy_DPFDT(0:imax-1, 1:jmax)


    real(DP), parameter :: &
      & c      = 2.99792458d8           , &
      & planc  = 6.6260755d-34          , &
      & boltz  = 1.380658d-23
    real(DP) :: xy_ExpTerm(0:imax-1, 1:jmax)
    real(DP) :: xy_PF     (0:imax-1, 1:jmax)


    xy_ExpTerm(:,:) = exp( planc * c * ( WN + 1.0d-10 ) / ( boltz * xy_Temp ) )

    xy_PF(:,:) = 2.0d0 * planc * c * c * WN * WN * WN &
      / ( xy_ExpTerm(:,:) - 1.0d0 )

    xy_DPFDT(:,:) = &
      & 1.0d0 / ( 2.0d0 * c * WN * WN * boltz ) &
      & * ( xy_pf(:,:) / xy_Temp(:,:) )**2 &
      & * xy_ExpTerm(:,:)


  end subroutine DPFDT_ARRAY2D

  !**************************************************************************
  ! subroutine gauleg
  !**************************************************************************
  ! This code is the same as sh_gauleg subrouting in sh_module.f90.
  ! This code was obtained from Press et al. [1993], p141.
  !**************************************************************************

  subroutine gauleg( x1, x2, n, x, w )

    ! 
    ! Physical constants settings
    !
    use constants, only: PI   ! $ \pi $ .
                              ! Circular constant

    real(DP), intent(in ) :: x1,x2
    integer , intent(in ) :: n
    real(DP), intent(out) :: x(n),w(n)


    real(DP):: eps
    real(DP):: z1, z, xm, xl, pp, p3, p2, p1
    integer :: i, j, m


    !
    ! changed at 2005/09/14
    !
!!$      eps = 1.0e-11
    eps = 1.0e-15

    m  = ( n+1 ) / 2
    xm = ( x2+x1 ) / 2.0d0
    xl = ( x2-x1 ) / 2.0d0

    do i = 1, m
      z = cos( pi*dble(i-0.25)/(dble(n+0.5)) )

100   continue
      p1 = 1.0d0
      p2 = 0.0d0
!!$         do j = 1, n
      do j = 1, n
        p3 = p2
        p2 = p1
        p1 = ( (2.0d0*dble(j)-1.0d0)*z*p2-(dble(j)-1.0d0)*p3 ) / dble(j)
      enddo

      pp = dble(n) * (z*p1-p2) / (z*z-1.0d0)
      z1 = z
      z  = z1 - p1 / pp

      !
      ! changed at 2005/09/14
      !
!!$         if( abs( z-z1 ) .gt. eps ) go to 100
!!$         if( abs( z-z1 ) / z1 .gt. eps ) go to 100
      !
      ! changed before 2008/08/03
      !
      if( abs( z-z1 ) / abs( z1 + 1.0d-200 ) .gt. eps ) go to 100

      x( i     ) = xm-xl*z
      x( n+1-i ) = xm+xl*z
      w( i     ) = 2.0d0*xl/((1.0d0-z*z)*pp*pp)
      w( n+1-i ) = w(i)
    end do


  end subroutine gauleg

  !**************************************************************************

  subroutine RadiationRoeweLiouInit


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

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

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

    ! 日付および時刻の取り扱い
    ! Date and time handler
    !
    use dc_date, only: DCDiffTimeCreate

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

    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

    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


    namelist /radiation_RL78_nml/ &
      & VMRCO2,                   &
      & DelTimeCalcTransValue,    &
      & DelTimeCalcTransUnit,     &
      & flag_save_time

    VMRCO2                = 382.0d-6

    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 = radiation_RL78_nml,       & ! (out)
        & iostat = iostat_nml )             ! (out)
      close( unit_nml )

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

    ! Handle interval time
    !
    call DCDiffTimeCreate( IntTimeSave,               & ! (out)
      & DelTimeCalcTransValue, DelTimeCalcTransUnit )   ! (in)


!!$    do n = 1, nbmax
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
!!$      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
!!$    end do



    do n = 1, 12-1
      ! unit conversion from (cm-1) to (m-1)
      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
      ! unit conversion from (g-1 cm2) to (kg-1 m2)
      a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
    end do
    ! O3 unit conversion
    do n = 12, 21
      ! unit conversion from (cm-1) to (m-1)
      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
      ! unit conversion from "wrong (g-1 cm2)" to (kg-1 m2)
      !   The values written by Roewe and Liou (1978) seems to be wrong.
      !   In order to convert cm-1 atm-1 to g-1 cm2, the value has to be MULTIPLIED
      !   by ~400.
      !   The factor ~400 is
      !   1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0 * 1.0d4 * 1.0d-3.
      !   However, it seems that the value was DIVIDED by ~400 to convert the value
      !   in cm-1 atm-1 to g-1 cm2 by Roewe and Liou (1978).
      !   So, first, the value written in Roewe and Liou (1978) is multiplied by ~400
      !   to obtain the correct value in cm-1 atm-1.
      !   Then the value is multiplied by ~40 to convert unit from cm-1 atm-1 to m2 kg-1.
      !   The factor ~40 is
      !   1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0. 
      a_BandParam(3,n) = a_BandParam(3,n) &
        & * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0 &
        & * 1.0d4 * 1.0d-3 &
        & * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0
    end do
    do n = 21+1, nbmax
      ! unit conversion from (cm-1) to (m-1)
      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
      ! unit conversion from (g-1 cm2) to (kg-1 m2)
      a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
    end do


!!$    do n = 1, nbH2Omax
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_H2OBandParam(1,n) = a_H2OBandParam(1,n) * 1.0d2
!!$      a_H2OBandParam(2,n) = a_H2OBandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_H2OBandParam(3,n) = a_H2OBandParam(3,n) * 1.0d-4 * 1.0d3
!!$
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_H2OContParam(1,n) = a_H2OContParam(1,n) * 1.0d2
!!$      a_H2OContParam(2,n) = a_H2OContParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2 atm-1) to (kg-1 m2 Pa-1)
!!$      a_H2OContParam(3,n) = a_H2OContParam(3,n) * 1.0d-4 * 1.0d3 * 1.0d-5
!!$    end do

    ! unit conversion from (cm-1) to (m-1)
    a_CO2BandParam(1) = a_CO2BandParam(1) * 1.0d2
    a_CO2BandParam(2) = a_CO2BandParam(2) * 1.0d2
    ! unit conversion from (g-1 cm2) to (kg-1 m2)
    a_CO2BandParam(3) = a_CO2BandParam(3) * 1.0d-4 * 1.0d3

!!$    do n = 1, nbO3max
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_O3BandParam(1,n) = a_O3BandParam(1,n) * 1.0d2
!!$      a_O3BandParam(2,n) = a_O3BandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_O3BandParam(3,n) = a_O3BandParam(3,n) * 1.0d-4 * 1.0d3
!!$    end do

!!$    do n = 1, 3
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_TMPFORO3BandParam(1,n) = a_TMPFORO3BandParam(1,n) * 1.0d2
!!$      a_TMPFORO3BandParam(2,n) = a_TMPFORO3BandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_TMPFORO3BandParam(3,n) = a_TMPFORO3BandParam(3,n) * 1.0d-4 * 1.0d3
!!$    end do


    ! Values by Miyoshi and Morita (19??)
    !
    H2OScaleIndex = 1.0d0
    CO2ScaleIndex = 0.9d0
    O3ScaleIndex  = 0.4d0


    NGaussQuad = 5


    MeanMolWeight = 28.0d-3
    H2OMolWeight  = 18.0d-3
    CO2MolWeight  = 44.0d-3


    allocate( xyrra_TransSave  (0:imax-1, 1:jmax, 0:kmax, 0:kmax, 1:nbmax) )


    ! Preparation of tables for calculation of Plank function
    !
    TableTempMin       =  50.0d0
    TableTempMax       = 400.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( &
        & a_BandParam(1,n), a_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 CalcIntegratedPF2D( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, &
          imax, jmax, xy_TempTMP, xy_PF )
        call Integ_DPFDT_GQ_Array2D( &
          & a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, 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

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


    do n = 1, nRefAtm
      RefAtmPro(2,n) = RefAtmPro(2,n) * 1.0d2
!!$      RefAtmPro(5,n) = RefAtmPro(5,n) * 1.0d-3 / ( 48.0d0 * 1.66d-27 ) &
!!$        & / ( RefAtmPro(2,n) / ( 1.38d-23 * RefAtmPro(3,n) ) )
      RefAtmPro(5,n) = RefAtmPro(5,n) * 1.0d-3 &
        & / ( RefAtmPro(2,n) / ( 8.314d0 / 48.0d-3 * RefAtmPro(3,n) ) )
    end do



    sw_fs = .false.


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


  end subroutine RadiationRoeweLiouInit



  !****************************************************************************

  subroutine SetRefAtmPro(      &
    & index, sw_log, xyz_Press, &
    & xyz_Array                 &
    & )

    use dc_types
    use gridset, only: imax, &
      &                jmax, &
      &                kmax

    integer , intent(in ) :: index
    logical , intent(in ) :: sw_log
    real(DP), intent(in ) :: xyz_Press(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_Array(0:imax-1,1:jmax,1:kmax)

    !
    ! local variables
    !
    real(DP):: z_RefAtmProArray(1:nRefAtm)
    integer :: k, kk


    if ( sw_log ) then
      do k = 1, nRefAtm
        z_RefAtmProArray(k) = log( RefAtmPro(index,k) + 1.0d-100 )
      end do
    else
      do k = 1, nRefAtm
        z_RefAtmProArray(k) =      RefAtmPro(index,k)
      end do
    end if


    do k = 1, kmax
      if( xyz_Press(0,1,k) <= RefAtmPro(2,nRefAtm) ) then
        xyz_Array(0,1,k) = z_RefAtmProArray(k)
      else
        search_loop : do kk = 2, nRefAtm
          if( RefAtmPro(2,kk) < xyz_Press(0,1,k) ) exit search_loop
        end do search_loop
        if( kk > nRefAtm ) &
          stop 'Unexpected error in setting temperature profile'
        xyz_Array(0,1,k) = ( z_RefAtmProArray( kk ) - z_RefAtmProArray( kk-1 ) ) &
          / ( log( RefAtmPro(2,kk)  / RefAtmPro(2,kk-1) ) ) &
          * ( log( xyz_Press(0,1,k) / RefAtmPro(2,kk-1) ) ) &
          + z_RefAtmProArray( kk-1 )
      end if
    end do


    if ( sw_log ) then
      do k = 1, kmax
        xyz_Array(0,1,k) = exp( xyz_Array(0,1,k) )
      end do
    else
      do k = 1, kmax
        xyz_Array(0,1,k) =      xyz_Array(0,1,k)
      end do
    end if

    do k = 1, kmax
      xyz_Array(:,:,k) = xyz_Array(0,1,k)
    end do


  end subroutine SetRefAtmPro



end module radiation_RL78
