!= Chou et al. (2001) による長波放射モデル
!
!= Long radiation model described by Chou et al. (2001)
!
! Authors::   Yoshiyuki O. TAKAHASHI
! Version::   $Id: rad_C2001.f90,v 1.4 2013/01/27 11:31:14 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module rad_C2001
  !
  != Chou et al. (2001) による長波放射モデル
  !
  != Long radiation model described by Chou et al. (2001)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 長波放射モデル.
  !
  ! This is a model of long wave radiation. 
  !
  !== 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
  !
!!$  ! 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.

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

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


  integer , parameter         :: nbmax = 10

  real(DP)             , save :: aa_BandWN              (1:2, 1:nbmax)
  real(DP)             , save :: aa_BandIntegPFCoefs    (0:5, 1:nbmax)
  real(DP)             , save :: aa_BandH2OSclParams    (1:5, 1:nbmax)
  real(DP)             , save :: aa_BandH2OLineParams   (1:8, 1:nbmax)
  real(DP)             , save :: aa_BandH2OContParams   (1:1, 1:nbmax)

  integer , parameter         :: nsbmax = 3
  real(DP)             , save :: aa_SubBandWN           (1:2, 1:nsbmax)
  real(DP)             , save :: aa_SubBandH2OSclParams (1:5, 1:nsbmax)
  real(DP)             , save :: aa_SubBandH2OLineParams(1:8, 1:nsbmax)
  real(DP)             , save :: aa_SubBandH2OContParams(1:1, 1:nsbmax)
  real(DP)             , save :: aa_SubBandCO2SclParams (1:5, 1:2)
  real(DP)             , save :: aa_SubBandCO2LineParams(1:8, 1:2)

  real(DP)             , save :: aa_BandN2OSclParams (1:5, 1:nbmax)
  real(DP)             , save :: aa_BandN2OLineParams(1:8, 1:nbmax)
  real(DP)             , save :: aa_BandCH4SclParams (1:5, 1:nbmax)
  real(DP)             , save :: aa_BandCH4LineParams(1:8, 1:nbmax)

  real(DP)             , save :: aa_WeakBandCO2SclParams (1:5, 1:nbmax)
  real(DP)             , save :: aa_WeakBandCO2LineParams(1:8, 1:nbmax)

  integer , parameter         :: ncloudparam = 4
  real(DP)             , save :: aa_BandCloudIceExtParams(1:ncloudparam, 1:nbmax)
  real(DP)             , save :: aa_BandCloudIceSSAParams(1:ncloudparam, 1:nbmax)
  real(DP)             , save :: aa_BandCloudIceAFParams (1:ncloudparam, 1:nbmax)
  real(DP)             , save :: aa_BandCloudLiqExtParams(1:ncloudparam, 1:nbmax)
  real(DP)             , save :: aa_BandCloudLiqSSAParams(1:ncloudparam, 1:nbmax)
  real(DP)             , save :: aa_BandCloudLiqAFParams (1:ncloudparam, 1:nbmax)


  ! MEMO:
  ! Bands range from 0 to 3000 cm-1.
  !

  ! Table 1
  !
  data aa_BandWN &
    & / &
    &    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 + O3 + CO2
    & 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.
    & /

  ! Table 2
  !
  data aa_BandIntegPFCoefs &
    & / &
    ! c0         c1         c2         c3         c4         c5
    &  5.344e+0, -2.062e-1,  2.533e-3, -6.863e-6,  1.012e-8, -6.267e-12, & ! 1
    &  2.715e+1, -5.404e-1,  2.950e-3,  2.723e-7, -9.338e-9,  9.968e-12, & ! 2
    & -3.486e+1,  1.113e+0, -1.301e-2,  6.496e-5, -1.182e-7,  8.042e-11, & ! 3
    & -6.051e+1,  1.409e+0, -1.208e-2,  4.405e-5, -5.674e-8,  2.566e-11, & ! 4
    & -2.669e+1,  5.283e-1, -3.445e-3,  6.072e-6,  1.252e-8, -2.155e-11, & ! 5
    & -6.727e+0,  4.226e-2,  1.044e-3, -1.292e-5,  4.740e-8, -4.486e-11, & ! 6
    &  1.879e+1, -5.836e-1,  6.968e-3, -3.939e-5,  1.012e-7, -8.230e-11, & ! 7
    &  1.034e+2, -2.513e+0,  2.375e-2, -1.069e-4,  2.184e-7, -1.370e-10, & ! 8
    & -1.048e+1,  3.821e-1, -5.227e-3,  3.441e-5, -1.108e-7,  1.409e-10, & ! 9
    &  1.677e+0,  6.540e-2, -1.813e-3,  1.291e-5, -2.672e-8,  1.979e-11  & !10
    & /


  ! Table 3
  !
  data aa_BandH2OSclParams &
    & / &
    !  pr        Tr        m         a         b
    &  500.0d2 , 250.0d0 , 1.0d0   , 0.0021d0, -1.01d-5, & ! 1:H2O
    &  500.0d2 , 250.0d0 , 1.0d0   , 0.0140d0,  5.57d-5, & ! 2:H2O
    &  -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3:H2O + CO2
    &  500.0d2 , 250.0d0 , 1.0d0   , 0.0302d0,  2.96d-4, & ! 4:H2O
    &  500.0d2 , 250.0d0 , 1.0d0   , 0.0307d0,  2.86d-4, & ! 5:H2O + O3
    &  500.0d2 , 250.0d0 , 0.77d0  , 0.0195d0,  1.11d-4, & ! 6:H2O
    &  500.0d2 , 250.0d0 , 0.5d0   , 0.0152d0,  7.61d-5, & ! 7:H2O
    &  500.0d2 , 250.0d0 , 1.0d0   , 0.0008d0, -3.52d-6, & ! 8:H2O
    &  500.0d2 , 250.0d0 , 1.0d0   , 0.0096d0,  1.64d-5, & ! 9:H2O
    &  500.0d2 , 250.0d0 , 1.0d0   , 0.0149d0,  6.20d-5  & !10:H2O
    & /


  ! Table 4
  !
  data aa_BandH2OLineParams &
    & / &
    ! k1/mu    n   dg1       dg2       dg3       dg4       dg5       dg6
    & 2.96d1 ,  6, 0.2747d0, 0.2717d0, 0.2752d0, 0.1177d0, 0.0352d0, 0.0255d0, & ! 1
    & 4.17d-1,  6, 0.1521d0, 0.3974d0, 0.1778d0, 0.1826d0, 0.0374d0, 0.0527d0, & ! 2
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3
    & 5.25d-4,  6, 0.4654d0, 0.2991d0, 0.1343d0, 0.0646d0, 0.0226d0, 0.0140d0, & ! 4
    & 5.25d-4,  6, 0.5543d0, 0.2723d0, 0.1131d0, 0.0443d0, 0.0160d0, -1.0d100, & ! 5
    & 9.37d-3,  8, 0.5955d0, 0.2693d0, 0.0953d0, 0.0335d0, 0.0064d0, -1.0d100, & ! 6
    & 4.72d-2,  9, 0.1958d0, 0.3469d0, 0.3147d0, 0.1013d0, 0.0365d0, 0.0048d0, & ! 7
    & 1.32d0 ,  6, 0.0740d0, 0.1636d0, 0.4174d0, 0.1783d0, 0.1101d0, 0.0566d0, & ! 8
    & 5.25d-4, 16, 0.1437d0, 0.2197d0, 0.3185d0, 0.2351d0, 0.0647d0, 0.0183d0, & ! 9
    & 1.06d-1,  8, 0.3153d0, 0.4604d0, 0.1326d0, 0.0798d0, 0.0119d0, -1.0d100  & !10
    & /

  ! Table 9
  !
  data aa_BandH2OContParams &
    & / &
    ! k/mu
    & -1.0d100, & ! 1
    & -1.0d100, & ! 2
    & -1.0d100, & ! 3
    & 15.8d0  , & ! 4
    &  9.40d0 , & ! 5
    &  7.75d0 , & ! 6
    &  8.78d0 , & ! 7
    & -1.0d100, & ! 8
    & -1.0d100, & ! 9
    & -1.0d100  & !10
    & /


  ! Table 1
  !
  data aa_SubBandWN &
    & / &
    &  540.0d0,  620.0d0, & ! 1:H2O + CO2
    &  620.0d0,  720.0d0, & ! 2:H2O + CO2
    &  720.0d0,  800.0d0  & ! 3:H2O + CO2
    & /

  ! Table 3
  !
  data aa_SubBandH2OSclParams &
    & / &
    !  pr       Tr       m       a          b
    &  500.0d2, 250.0d0, 1.0d0 , 0.0167d0,  8.54d-5, & ! 1
    &  500.0d2, 250.0d0, 1.0d0 , 0.0167d0,  8.54d-5, & ! 2
    &  500.0d2, 250.0d0, 1.0d0 , 0.0167d0,  8.54d-5  & ! 3
    & /

  ! Table 10 (H2O)
  !
  data aa_SubBandH2OLineParams &
    & / &
    ! k1        n  dg1       dg2       dg3       dg4       dg5       dg6
    & 1.33d-2 , 8, 0.1782d0, 0.0593d0, 0.0215d0, 0.0068d0, 0.0022d0, -1.0d100, & ! a
    & 1.33d-2 , 8, 0.0923d0, 0.1675d0, 0.0923d0, 0.0187d0, 0.0178d0, -1.0d100, & ! b
    & 1.33d-2 , 8, 0.0000d0, 0.1083d0, 0.1581d0, 0.0455d0, 0.0274d0, 0.0041d0  & ! c
    & /

  ! Table 9
  !
  data aa_SubBandH2OContParams &
    & / &
    ! k/mu
    & 109.6d0 , & ! 1
    &  54.8d0 , & ! 2
    &  27.4d0   & ! 3
    & /

  ! Table 3
  !   Memo: Parameters to scale absorber amount are given for three subbands 
  !         in Table 3 (below). But, parameters for k-distribution methods, 
  !         i.e., absorption coefficients and width of bins, are only given 
  !         for two bands of "Wings" and "Center" in Table 10.
  !         So, k-distribution parameters in 3c is regarded as same as those 
  !         in 3a. 
  !         See brief description about this in p. 23.
!!$  data aa_SubBandCO2SclParams &
!!$    & / &
!!$    !  pr       Tr       m      a          b
!!$    &  300.0d2 , 250.0d0 , 0.5d0   , 0.0179d0,  1.02d-4, & ! 1:H2O + CO2
!!$    &   30.0d2 , 250.0d0 , 0.85d0  , 0.0042d0,  2.00d-5, & ! 2:H2O + CO2
!!$    &  300.0d2 , 250.0d0 , 0.5d0   , 0.0184d0,  1.12d-4  & ! 3:H2O + CO2
!!$    & /
  data aa_SubBandCO2SclParams &
    & / &
    !  pr        Tr        m       a          b
    &  300.0d2 , 250.0d0 , 0.5d0 , 0.0179d0,  1.02d-4, & ! Wings
    &   30.0d2 , 250.0d0 , 0.85d0, 0.0042d0,  2.00d-5  & ! Center
    & /

  ! Table 10 (CO2)
  !
  data aa_SubBandCO2LineParams &
    & / &
    ! k1        n  dg1       dg2       dg3       dg4       dg5       dg6
    & 2.66d-5,  8, 0.1395d0, 0.1407d0, 0.1549d0, 0.1357d0, 0.0182d0, 0.0220d0, & ! Wings
    & 2.66d-3,  8, 0.0766d0, 0.1372d0, 0.1189d0, 0.0335d0, 0.0169d0, 0.0059d0  & ! Center
    & /


  ! Table 5 (N2O)
  !
  data aa_BandN2OLineParams &
    & / &
    ! k1       n   dg1       dg2       dg3       dg4       dg5       dg6
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 1
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 2
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 4
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 5
    & 6.32d-2, 21, 0.9404d0, 0.0596d0, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 6
    & 5.36d-2,  8, 0.5620d0, 0.1387d0, 0.2406d0, 0.0587d0, -1.0d100, -1.0d100, & ! 7
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 8
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 9
    & 2.52d-1, 58, 0.971d0 , 0.029d0 , -1.0d100, -1.0d100, -1.0d100, -1.0d100  & !10
    & /
  data aa_BandN2OSclParams &
    & / &
    !  pr        Tr        m       a          b
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 1
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 2
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 4
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 5
    &  500.0d2 , 250.0d0 , 0.0d0 , 1.93d-3 ,  4.38d-6, & ! 6
    &  500.0d2 , 250.0d0 , 0.48d0, 1.38d-3 ,  7.48d-6, & ! 7
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 8
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 9
    &  500.0d2 , 250.0d0 , 0.0d0 , 1.45d-3 ,  3.67d-6  & !10
    & /

  ! Table 5 (CH4)
  !
  data aa_BandCH4LineParams &
    & / &
    ! k1       n   dg1       dg2       dg3       dg4       dg5       dg6
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 1
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 2
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 4
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 5
    & 5.81d-3,  1, 1.0d0   , -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 6
    & 6.29d-2, 12, 0.6107d0, 0.2802d0, 0.1073d0, 0.0018d0, -1.0d100, -1.0d100, & ! 7
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 8
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 9
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100  & !10
    & /
  data aa_BandCH4SclParams &
    & / &
    !  pr        Tr        m       a          b
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 1
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 2
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 4
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 5
    &  500.0d2 , 250.0d0 , 0.0d0 , 1.70d-2 ,  1.58d-4, & ! 6
    &  500.0d2 , 250.0d0 , 0.65d0, 5.96d-4 , -2.29d-6, & ! 7
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 8
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 9
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100  & !10
    & /

  ! Table 6
  !
  data aa_WeakBandCO2LineParams &
    & / &
    ! k1       n   dg1       dg2       dg3       dg4       dg5       dg6
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 1
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 2
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3
    & 1.92d-7,  5, 0.1216d0, 0.2436d0, 0.2498d0, 0.2622d0, 0.0781d0, 0.0427d0, & ! 4
    & 1.92d-7,  5, 0.0687d0, 0.1480d0, 0.1951d0, 0.3344d0, 0.1720d0, 0.0818d0, & ! 5
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 6
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 7
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 8
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 9
    & 2.66d-5,  8, 0.267d0 , 0.220d0 , 0.210d0 , 0.241d0 , 0.020d0 , 0.042d0   & !10
    & /
  data aa_WeakBandCO2SclParams &
    & / &
    !  pr        Tr        m       a          b
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 1
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 2
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 3
    &  500.0d2 , 250.0d0 , 0.0d0 , 3.58d-2 ,  4.05d-4, & ! 4
    &  500.0d2 , 250.0d0 , 0.0d0 , 3.43d-2 ,  3.74d-4, & ! 5
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 6
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 7
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 8
    & -1.0d100, -1.0d100, -1.0d100, -1.0d100, -1.0d100, & ! 9
    &  300.0d2 , 250.0d0 , 0.5d0 , 1.79d-2 ,  1.02d-4  & !10
    & /



  ! Table 11a
  !   Unit of effective radius of particle is micron meter.
  !   Unit of extinction parameters would be g-1 m2 ?
  data aa_BandCloudIceExtParams &
    & / &
    & -0.44171, 0.62951, 0.06465, -1.0d100, &
    & -0.13727, 0.61291, 0.28962, -1.0d100, &
    & -0.01878, 1.67680, 0.79080, -1.0d100, &
    & -0.01896, 1.06510, 0.69493, -1.0d100, &
    & -0.04788, 0.88178, 0.54492, -1.0d100, &
    & -0.02265, 1.57390, 0.76161, -1.0d100, &
    & -0.01038, 2.15640, 0.89045, -1.0d100, &
    & -0.00450, 2.51370, 0.95989, -1.0d100, &
    & -0.00044, 3.15050, 1.03750, -1.0d100, &
    & -0.02956, 1.44680, 0.71283, -1.0d100  &
    & /

  ! Table 11b
  !   Unit of effective radius of particle is micron meter.
  !   Unit of extinction parameters would be g-1 m2 ?
  data aa_BandCloudLiqExtParams &
    & / &
    & 0.08641,  0.01769, -1.5572E-3,  3.4896E-5, &
    & 0.22027,  0.00997, -1.8719E-3,  5.3112E-5, &
    & 0.38074, -0.03027,  1.0154E-3, -1.1849E-5, &
    & 0.15587,  0.00371, -7.7705E-4,  2.0547E-5, &
    & 0.05518,  0.04544, -4.2067E-3,  1.0184E-4, &
    & 0.12724,  0.04751, -5.2037E-3,  1.3711E-4, &
    & 0.30390,  0.01656, -3.5271E-3,  1.0828E-4, &
    & 0.63617, -0.06287,  2.2350E-3, -2.3177E-5, &
    & 1.15470, -0.19282,  1.2084E-2, -2.5612E-4, &
    & 0.34021, -0.02805,  1.0654E-3, -1.5443E-5  &
    & /

  ! Table 12a
  !   Unit of effective radius of particle is micron meter.
  data aa_BandCloudIceSSAParams &
    & / &
    & 0.17201,  1.2229E-2, -1.4837E-4,  5.8020E-7, &
    & 0.81470, -2.7293E-3,  9.7816E-8,  5.7650E-8, &
    & 0.54859, -4.8273E-4,  5.4353E-6, -1.5679E-8, &
    & 0.39218,  4.1717E-3, -4.8869E-5,  1.9144E-7, &
    & 0.71773, -3.3640E-3,  1.9713E-5, -3.3189E-8, &
    & 0.77345, -5.5228E-3,  4.8379E-5, -1.5151E-7, &
    & 0.74975, -5.6604E-3,  5.6475E-5, -1.9664E-7, &
    & 0.69011, -4.5348E-3,  4.9322E-5, -1.8255E-7, &
    & 0.83963, -6.7253E-3,  6.1900E-5, -2.0862E-7, &
    & 0.64860, -2.8692E-3,  2.7656E-5, -8.9680E-8  &
    & /

  ! Table 12b
  !   Unit of effective radius of particle is micron meter.
  data aa_BandCloudLiqSSAParams &
    & / &
    &  0.17201,  1.2229E-2, -1.4837E-4, 5.8020E-7, &
    & -0.01338,  9.3134E-2, -6.0491E-3, 1.3059E-4, &
    &  0.03710,  7.3211E-2, -4.4211E-3, 9.2448E-5, &
    & -0.00376,  9.3344E-2, -5.6561E-3, 1.1387E-4, &
    &  0.40212,  7.8083E-2, -5.9583E-3, 1.2883E-4, &
    &  0.57928,  5.9094E-2, -5.4425E-3, 1.2725E-4, &
    &  0.68974,  4.2334E-2, -4.9469E-3, 1.2863E-4, &
    &  0.80122,  9.4578E-3, -2.8508E-3, 9.0078E-5, &
    &  1.02340, -2.6204E-2,  4.2552E-4, 3.2160E-6, &
    &  0.05092,  7.5409E-2, -4.7305E-3, 1.0121E-4  &
    & /

  ! Table 13a
  !   Unit of effective radius of particle is micron meter.
  data aa_BandCloudIceAFParams &
    & / &
    & 0.57867, 1.0135E-2, -1.1142E-4, 4.1537E-7, &
    & 0.72259, 3.1149E-3, -1.9927E-5, 5.6024E-8, &
    & 0.76109, 4.5449E-3, -4.6199E-5, 1.6446E-7, &
    & 0.86934, 2.7474E-3, -3.1301E-5, 1.1959E-7, &
    & 0.89103, 1.8513E-3, -1.6551E-5, 5.5193E-8, &
    & 0.86325, 2.1408E-3, -1.6846E-5, 4.9473E-8, &
    & 0.85064, 2.5028E-3, -2.0812E-5, 6.3427E-8, &
    & 0.86945, 2.4615E-3, -2.3882E-5, 8.2431E-8, &
    & 0.80122, 3.1906E-3, -2.4856E-5, 7.2411E-8, &
    & 0.73290, 4.8034E-3, -4.4425E-5, 1.4839E-7  &
    & /

  ! Table 13b
  !   Unit of effective radius of particle is micron meter.
  data aa_BandCloudLiqAFParams &
    & / &
    & -0.51930,  0.20290  , -1.1747E-2,  2.3868E-4, &
    & -0.22151,  0.19708  , -1.2462E-2,  2.6646E-4, &
    &  0.14157,  0.14705  , -9.5802E-3,  2.0819E-4, &
    &  0.41590,  0.10482  , -6.9118E-3,  1.5115E-4, &
    &  0.55338,  7.7016E-2, -5.2218E-3,  1.1587E-4, &
    &  0.61384,  6.4402E-2, -4.6241E-3,  1.0746E-4, &
    &  0.67891,  4.8698E-2, -3.7021E-3,  9.1966E-5, &
    &  0.78169,  2.0803E-2, -1.4749E-3,  3.9362E-5, &
    &  0.93218, -3.3425E-2,  2.9632E-3, -6.9362E-5, &
    &  0.01649,  0.16561  , -1.0723E-2,  2.3220E-4  &
    & /


  ! 
  ! Public procedure
  !
  public :: RadC2001CalcTransBand3CO2
  public :: RadC2001CalcTransBand3H2O
  public :: RadC2001CalcTrans

  public :: RadC2001CalcCloudOptProp
  public :: RadC2001ReduceCloudOptDep

  public :: RadC2001CalcIntegratedPF3D
  public :: RadC2001CalcIntegratedPF2D

  public :: RadC2001Init


  real(DP), save :: MeanMolWeight = 28.97d-3
  real(DP), save :: H2OMolWeight  = 18.0d-3
  real(DP), save :: CO2MolWeight  = 44.0d-3



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


  character(*), parameter:: module_name = 'rad_C2001'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_C2001.f90,v 1.4 2013/01/27 11:31:14 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version


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

contains

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

  subroutine RadC2001CalcCloudOptProp(       &
    & Spec, iband, xyz_REff,                       & ! (in)
    & xyz_ExtCoef, xyz_SSA, xyz_AF                 & ! (out)
    & )

    ! USE statements
    !


    character(len=*), intent(in ) :: SPEC
    integer         , intent(in ) :: iband
    real(DP)        , intent(in ) :: xyz_REff   (0:imax-1, 1:jmax, 1:kmax)
    real(DP)        , intent(out) :: xyz_ExtCoef(0:imax-1, 1:jmax, 1:kmax)
    real(DP)        , intent(out) :: xyz_SSA    (0:imax-1, 1:jmax, 1:kmax)
    real(DP)        , intent(out) :: xyz_AF     (0:imax-1, 1:jmax, 1:kmax)


    !
    ! Work variables
    !
    real(DP) :: xyz_REffInMicron(0:imax-1, 1:jmax, 1:kmax)
    integer  :: l


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


    xyz_REffInMicron = xyz_REff * 1.0d6

    xyz_ExtCoef = 0.0_DP
    xyz_SSA     = 0.0_DP
    xyz_AF      = 0.0_DP

    if ( Spec == 'Ice' ) then
      ! Eq. (6.4a)
      xyz_ExtCoef = aa_BandCloudIceExtParams(1,iband) &
        &     + aa_BandCloudIceExtParams(2,iband) &
        &       / xyz_REffInMicron**aa_BandCloudIceExtParams(3,iband)
      !  Convert unit from g-1 m2 to kg-1 m2
      xyz_ExtCoef = xyz_ExtCoef * 1.0d3
      do l = 1, ncloudparam
        ! Eq. (6.5)
        xyz_SSA = xyz_SSA + aa_BandCloudIceSSAParams(l,iband) * xyz_REffInMicron**(l-1)
        ! Eq. (6.6)
        xyz_AF  = xyz_AF  + aa_BandCloudIceAFParams (l,iband) * xyz_REffInMicron**(l-1)
      end do
    else if ( Spec == 'Liquid' ) then
      do l = 1, ncloudparam
        ! Eq. (6.4b)
        xyz_ExtCoef = xyz_ExtCoef + aa_BandCloudLiqExtParams(l,iband) * xyz_REffInMicron**(l-1)
      end do
      !  Convert unit from g-1 m2 to kg-1 m2
      xyz_ExtCoef = xyz_ExtCoef * 1.0d3
      do l = 1, ncloudparam
        ! Eq. (6.5)
        xyz_SSA = xyz_SSA + aa_BandCloudLiqSSAParams(l,iband) * xyz_REffInMicron**(l-1)
        ! Eq. (6.6)
        xyz_AF  = xyz_AF  + aa_BandCloudLiqAFParams (l,iband) * xyz_REffInMicron**(l-1)
      end do
    else
      call MessageNotify( 'E', module_name, 'Unsupported specie, %c', c1 = trim( Spec ) )
    end if


  end subroutine RadC2001CalcCloudOptProp

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

  subroutine RadC2001ReduceCloudOptDep(      &
    & xyz_SSA, xyz_AF,                             & ! (in)
    & xyz_DelOptDep                                & ! (inout)
    & )

    ! USE statements
    !


    real(DP), intent(in   ) :: xyz_SSA      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_AF       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout) :: xyz_DelOptDep(0:imax-1, 1:jmax, 1:kmax)


    !
    ! Work variables
    !
    real(DP) :: xyz_FuncF(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: a_CoefA(1:4)
    integer  :: l

    data a_CoefA / 0.5d0, 0.3738d0, 0.0076d0, 0.1186d0 /


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


    ! Extinction coefficient is modified to account for effect of scattering.

    ! Eq. (6.12)
    !
    xyz_FuncF = 0.0_DP
    do l = 1, 4
      xyz_FuncF = xyz_FuncF + a_CoefA(l) * xyz_AF**(l-1)
    end do

    ! Eq. (6.11)
    !
    xyz_DelOptDep = ( 1.0_DP - xyz_SSA * xyz_FuncF ) * xyz_DelOptDep


  end subroutine RadC2001ReduceCloudOptDep

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

  subroutine RadC2001CalcTransBand3CO2(           &
    & xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2,    & ! (in)
    & xyrr_Trans                                        & ! (out)
    & )

    ! USE statements
    !


    real(DP)        , intent(in ) :: xyz_DelCO2Mass(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 ) :: xyz_QCO2      (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_TransEach (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    character(len=STRING) :: Spec
    integer               :: iband
    integer               :: k
    integer               :: kk


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


    xyrr_Trans = 0.0_DP

    Spec = 'SubBandCO2'
    do iband = 1, 2
      call RadC2001CalcTrans(                               &
        & Spec, iband,                                      & ! (in)
        & xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2,    & ! (in)
        & xyrr_TransEach                                    & ! (out)
        & )

      do k = 0, kmax
        do kk = 0, kmax
          if ( k == kk ) then
            xyrr_Trans(:,:,k,kk) = 1.0_DP
          else
            xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk) + xyrr_TransEach(:,:,k,kk)
          end if
        end do
      end do
    end do


  end subroutine RadC2001CalcTransBand3CO2

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

  subroutine RadC2001CalcTransBand3H2O(                    &
    & xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
    & xyrr_Trans                                           & ! (out)
    & )

    ! USE statements
    !


    real(DP), intent(in ) :: xyz_DelH2OVapMass(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 ) :: xyz_QH2OVap      (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_TransCont (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    real(DP)              :: xyrr_TransLine (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    character(len=STRING) :: Spec
    integer               :: iband
    integer               :: k
    integer               :: kk


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


    xyrr_Trans = 0.0_DP

    do iband = 1, nsbmax
      Spec = 'SubBandH2OCont'
      call RadC2001CalcTrans(                                  &
        & Spec, iband,                                         & ! (in)
        & xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
        & xyrr_TransCont                                       & ! (out)
        & )
      Spec = 'SubBandH2OLine'
      call RadC2001CalcTrans(                                  &
        & Spec, iband,                                         & ! (in)
        & xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
        & xyrr_TransLine                                       & ! (out)
        & )

      do k = 0, kmax
        do kk = 0, kmax
          if ( k == kk ) then
            xyrr_Trans(:,:,k,kk) = 1.0_DP
          else
            xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk) &
              & + xyrr_TransCont(:,:,k,kk) * xyrr_TransLine(:,:,k,kk)
          end if
        end do
      end do
    end do


  end subroutine RadC2001CalcTransBand3H2O

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

  subroutine RadC2001CalcTrans(                         &
    & Spec, iband,                                      & ! (in)
    & xyz_DelAbsMass, xyz_Press, xyz_Temp, xyz_QMix,    & ! (in)
    & xyrr_Trans                                        & ! (out)
    & )

    ! USE statements
    !


    character(len=*), intent(in ) :: Spec
    integer         , intent(in ) :: iband
    real(DP)        , intent(in ) :: xyz_DelAbsMass(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 ) :: xyz_QMix      (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) :: RefPress
    real(DP) :: RefTemp
    real(DP) :: PressScaleIndex
    real(DP) :: CoefA
    real(DP) :: CoefB
    real(DP) :: CoefC            ! width of bin (dg) ! This notation is used in Chou et al. (1993)
    integer  :: CoefN
    real(DP) :: AbsCoef

    integer  :: NCoefC
    real(DP) :: a_CoefC(1:6)


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

    real(DP) :: xyz_TransOneLayer     (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyrr_TransElem        (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    integer  :: k
    integer  :: kk
    integer  :: l


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


    if ( Spec == 'H2OCont' ) then
      RefPress        =     -1.0d100
      RefTemp         =     -1.0d100
      PressScaleIndex =     -1.0d100
      CoefA           =     -1.0d100
      CoefB           =     -1.0d100
      AbsCoef         =      aa_BandH2OContParams(1  ,iband)
      a_CoefC(1:1)    =      1.0_DP
      a_CoefC(2:6)    =      1.0d100
      CoefN           = 1
      NCoefC          = 1
    else if ( Spec == 'H2OLine' ) then
      RefPress        =      aa_BandH2OSclParams (1  ,iband)
      RefTemp         =      aa_BandH2OSclParams (2  ,iband)
      PressScaleIndex =      aa_BandH2OSclParams (3  ,iband)
      CoefA           =      aa_BandH2OSclParams (4  ,iband)
      CoefB           =      aa_BandH2OSclParams (5  ,iband)
      AbsCoef         =      aa_BandH2OLineParams(1  ,iband)
      CoefN           = int( aa_BandH2OLineParams(2  ,iband) )
      a_CoefC         =      aa_BandH2OLineParams(3:8,iband)
      NCoefC          = 6
    else if ( Spec == 'SubBandH2OCont' ) then
      RefPress        =     -1.0d100
      RefTemp         =     -1.0d100
      PressScaleIndex =     -1.0d100
      CoefA           =     -1.0d100
      CoefB           =     -1.0d100
      AbsCoef         =      aa_SubBandH2OContParams(1  ,iband)
      a_CoefC(1:1)    =      1.0_DP
      a_CoefC(2:6)    =      1.0d100
      CoefN           = 1
      NCoefC          = 1
    else if ( Spec == 'SubBandH2OLine' ) then
      RefPress        =      aa_SubBandH2OSclParams (1  ,iband)
      RefTemp         =      aa_SubBandH2OSclParams (2  ,iband)
      PressScaleIndex =      aa_SubBandH2OSclParams (3  ,iband)
      CoefA           =      aa_SubBandH2OSclParams (4  ,iband)
      CoefB           =      aa_SubBandH2OSclParams (5  ,iband)
      AbsCoef         =      aa_SubBandH2OLineParams(1  ,iband)
      CoefN           = int( aa_SubBandH2OLineParams(2  ,iband) )
      a_CoefC         =      aa_SubBandH2OLineParams(3:8,iband)
      NCoefC          = 6
    else if ( Spec == 'SubBandCO2' ) then
      RefPress        =      aa_SubBandCO2SclParams (1  ,iband)
      RefTemp         =      aa_SubBandCO2SclParams (2  ,iband)
      PressScaleIndex =      aa_SubBandCO2SclParams (3  ,iband)
      CoefA           =      aa_SubBandCO2SclParams (4  ,iband)
      CoefB           =      aa_SubBandCO2SclParams (5  ,iband)
      AbsCoef         =      aa_SubBandCO2LineParams(1  ,iband)
      CoefN           = int( aa_SubBandCO2LineParams(2  ,iband) )
      a_CoefC         =      aa_SubBandCO2LineParams(3:8,iband)
      NCoefC          = 6
    else if ( Spec == 'N2OLine' ) then
      RefPress        =      aa_BandN2OSclParams (1  ,iband)
      RefTemp         =      aa_BandN2OSclParams (2  ,iband)
      PressScaleIndex =      aa_BandN2OSclParams (3  ,iband)
      CoefA           =      aa_BandN2OSclParams (4  ,iband)
      CoefB           =      aa_BandN2OSclParams (5  ,iband)
      AbsCoef         =      aa_BandN2OLineParams(1  ,iband)
      CoefN           = int( aa_BandN2OLineParams(2  ,iband) )
      a_CoefC         =      aa_BandN2OLineParams(3:8,iband)
      NCoefC          = 6
    else if ( Spec == 'CH4Line' ) then
      RefPress        =      aa_BandCH4SclParams (1  ,iband)
      RefTemp         =      aa_BandCH4SclParams (2  ,iband)
      PressScaleIndex =      aa_BandCH4SclParams (3  ,iband)
      CoefA           =      aa_BandCH4SclParams (4  ,iband)
      CoefB           =      aa_BandCH4SclParams (5  ,iband)
      AbsCoef         =      aa_BandCH4LineParams(1  ,iband)
      CoefN           = int( aa_BandCH4LineParams(2  ,iband) )
      a_CoefC         =      aa_BandCH4LineParams(3:8,iband)
      NCoefC          = 6
    else if ( Spec == 'WeakBandCO2Line' ) then
      RefPress        =      aa_WeakBandCO2SclParams (1  ,iband)
      RefTemp         =      aa_WeakBandCO2SclParams (2  ,iband)
      PressScaleIndex =      aa_WeakBandCO2SclParams (3  ,iband)
      CoefA           =      aa_WeakBandCO2SclParams (4  ,iband)
      CoefB           =      aa_WeakBandCO2SclParams (5  ,iband)
      AbsCoef         =      aa_WeakBandCO2LineParams(1  ,iband)
      CoefN           = int( aa_WeakBandCO2LineParams(2  ,iband) )
      a_CoefC         =      aa_WeakBandCO2LineParams(3:8,iband)
      NCoefC          = 6
    else
      call MessageNotify( 'E', module_name, 'Unsupported specie, %c', c1 = trim( Spec ) )
    end if

    if ( AbsCoef < 0.0_DP ) then
      xyrr_Trans = 1.0_DP
      return
    end if


    if ( ( Spec == 'H2OCont' ) .or. ( Spec == 'SubBandH2OCont' ) ) then
      call RadC2001ContScaleH2OAmt(                      &
        & xyz_DelAbsMass, xyz_Press, xyz_Temp, xyz_QMix, & ! (in)
        & xyz_DelAbsMassScaled                           & ! (out)
        & )
    else
      call RadC2001LineScaleAmt(                            &
        & xyz_DelAbsMass,                                   & ! (in)
        & xyz_Press, xyz_Temp,                              & ! (in)
        & RefPress, RefTemp, PressScaleIndex, CoefA, CoefB, & ! (in)
        & xyz_DelAbsMassScaled                              & ! (out)
        & )
    end if

    xyz_TransOneLayer = exp( - AbsCoef * xyz_DelAbsMassScaled )

    do k = 0, kmax
      kk = k
      xyrr_TransElem(:,:,k,kk) = 1.0_DP
      do kk = k+1, kmax
        xyrr_TransElem(:,:,k,kk) = xyrr_TransElem(:,:,k,kk-1) * xyz_TransOneLayer(:,:,kk)
      end do
    end do


    ! initialization
    do k = 0, kmax
      do kk = k+1, kmax
        xyrr_Trans(:,:,k,kk) = 0.0_DP
      end do
    end do

    ! Summation of fitted exponential functions
    Loop_Sum : do l = 1, NCoefC

      CoefC = a_CoefC(l)

      if ( CoefC > 0.0_DP ) then

        do k = 0, kmax
          do kk = k+1, kmax
            xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk) &
              & + CoefC * xyrr_TransElem(:,:,k,kk)
          end do
        end do

        ! Calculate fitted exponential function used in next loop
        if ( l /= NCoefC ) then
          do k = 0, kmax
            do kk = k+1, kmax
              xyrr_TransElem(:,:,k,kk) = xyrr_TransElem(:,:,k,kk)**CoefN
            end do
          end do
        end if

      end if

    end do Loop_Sum

    do k = 0, kmax
      do kk = k, k
        xyrr_Trans(:,:,k,kk) = 1.0_DP
      end do
    end do
    do k = 0, kmax
      do kk = 0, k-1
        xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,kk,k)
      end do
    end do


  end subroutine RadC2001CalcTrans

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

  subroutine RadC2001ContScaleH2OAmt(                   &
    & xyz_DelAbsMass, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
    & xyz_ContDelAbsMassScaled                          & ! (out)
    & )

    ! Equation (4.21), (4.19)

    ! USE statements
    !


    real(DP), intent(in ) :: xyz_DelAbsMass          (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 ) :: xyz_QH2OVap             (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_ContDelAbsMassScaled(0:imax-1, 1:jmax, 1:kmax)


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


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


    xyz_PressH2O            = xyz_Press * xyz_QH2OVap * MeanMolWeight / H2OMolWeight
    xyz_ContDelAbsMassScaled = xyz_DelAbsMass                      &
      & * xyz_PressH2O / 101325d0                                  & ! Pa to atm.
      & * exp( 1800.0d0 * ( 1.0d0 / xyz_Temp - 1.0d0 / 296.0d0 ) )


  end subroutine RadC2001ContScaleH2OAmt

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

  subroutine RadC2001LineScaleAmt(                      &
    & xyz_DelAbsMass,                                   & ! (in)
    & xyz_Press, xyz_Temp,                              & ! (in)
    & RefPress, RefTemp, PressScaleIndex, CoefA, CoefB, & ! (in)
    & xyz_DelAbsMassScaled                              & ! (out)
    & )

    ! Equation (4.4)

    ! USE statements
    !

    real(DP), intent(in ) :: xyz_DelAbsMass      (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 ) :: RefPress
    real(DP), intent(in ) :: RefTemp
    real(DP), intent(in ) :: PressScaleIndex
    real(DP), intent(in ) :: CoefA
    real(DP), intent(in ) :: CoefB
    real(DP), intent(out) :: xyz_DelAbsMassScaled(0:imax-1, 1:jmax, 1:kmax)


    !
    ! Work variables
    !


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


    xyz_DelAbsMassScaled = xyz_DelAbsMass                                     &
      & * ( xyz_Press / RefPress )**PressScaleIndex                           &
      & * ( 1.0_DP + CoefA * ( xyz_Temp - RefTemp ) + CoefB * ( xyz_Temp - RefTemp )**2 )


  end subroutine RadC2001LineScaleAmt

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

  subroutine RadC2001CalcIntegratedPF3D(  &
    & iband, km, xyz_Temp,                &
    & xyz_IntegPF,                        &
    & flag_DPFDT                          &
    & )

    ! USE statements
    !


    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

    !
    ! local variables
    !
    logical                     :: local_flag_DPFDT

!!$    integer                     :: xyz_TempIndex(0:imax-1, 1:jmax, 1:km)
    integer                     :: i
    integer                     :: j
    integer                     :: k
    integer                     :: l


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


    ! Temperature check
    !
    do k = 1, km
      do j = 1, jmax
        do i = 0, imax-1
          if ( ( xyz_Temp(i,j,k) < 150.0_DP ) .or. &
            &  ( xyz_Temp(i,j,k) > 350.0_DP ) ) then
            call MessageNotify( 'M', module_name, &
              & 'Temperature is not appropriate, Temp(%d,%d,%d) = %f.', &
              & i = (/i, j, k/), d = (/xyz_Temp(i,j,k)/) )
          end if
        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
      ! Calculation of integrated Planck function

      xyz_IntegPF = aa_BandIntegPFCoefs(0,iband)
      do l = 1, 5
        xyz_IntegPF = xyz_IntegPF + aa_BandIntegPFCoefs(l,iband) * xyz_Temp**l
      end do

    else
      ! Calculation of derivative of integrated Planck function

      xyz_IntegPF = aa_BandIntegPFCoefs(1,iband)
      do l = 2, 5
        xyz_IntegPF = xyz_IntegPF + aa_BandIntegPFCoefs(l,iband) * l * xyz_Temp**(l-1)
      end do

    end if


  end subroutine RadC2001CalcIntegratedPF3D

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

  subroutine RadC2001CalcIntegratedPF2D(  &
    & iband, xy_Temp,                     &
    & xy_IntegPF,                         &
    & 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)
    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)


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


    xyz_Temp(:,:,1) = xy_Temp

    call RadC2001CalcIntegratedPF3D(        &
      & iband, 1, xyz_temp,                 &
      & xyz_IntegPF,                        &
      & flag_DPFDT                          &
      & )

    xy_IntegPF = xyz_IntegPF(:,:,1)


  end subroutine RadC2001CalcIntegratedPF2D

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

  subroutine RadC2001Init


!!$    ! ファイル入出力補助
!!$    ! File I/O support
!!$    !
!!$    use dc_iounit, only: FileOpen
!!$
!!$    ! NAMELIST ファイル入力に関するユーティリティ
!!$    ! Utilities for NAMELIST file input
!!$    !
!!$    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid


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

    integer :: l

!!$    namelist /rad_C1991_nml/ &
!!$      & flag_save_time


    if ( rad_C2001_inited ) return


!!$    ! 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_dcpam_E_LW_V2_nml, & ! (out)
!!$        & iostat = iostat_nml )                ! (out)
!!$      close( unit_nml )
!!$
!!$      call NmlutilMsg( iostat_nml, module_name ) ! (in)
!!$    end if


    ! Convert unit from cm-1 to m-1
    aa_BandWN = aa_BandWN * 1.0d2

    ! Convert unit from gm-1 cm2 to kg-1 m2
    do l = 1, nbmax
      aa_BandH2OLineParams(1,l) = aa_BandH2OLineParams(1,l) * 1.0d-4 * 1.0d3
    end do
    do l = 1, nbmax
      if ( aa_BandH2OContParams(1,l) > 0.0_DP ) then
        aa_BandH2OContParams(1,l) = aa_BandH2OContParams(1,l) * 1.0d-4 * 1.0d3
      end if
    end do


    ! Convert unit from cm-1 to m-1
    aa_SubBandWN = aa_SubBandWN * 1.0d2

    ! Convert unit from gm-1 cm2 to kg-1 m2
    do l = 1, nsbmax
      aa_SubBandH2OLineParams(1,l) = aa_SubBandH2OLineParams(1,l) * 1.0d-4 * 1.0d3
    end do
    do l = 1, nsbmax
      aa_SubBandH2OContParams(1,l) = aa_SubBandH2OContParams(1,l) * 1.0d-4 * 1.0d3
    end do

    ! Convert unit from {(cm-atm)_{STP}}^{-1} to m2 kg-1
    !   MEMO: In a calculation below, GasRUniv variable in constants module can be used. 
    !         But, I do not use it, since unit of value is just converted by 
    !         multiplying a factor. Of course, non-use of GasRUniv must not cause 
    !         significant effect on result.
    do l = 1, 2
      aa_SubBandCO2LineParams(1,l) = aa_SubBandCO2LineParams(1,l) &
        & * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 44.0d-3 ) * 273.15d0
    end do


    ! 印字 ; 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) )


    rad_C2001_inited = .true.

  end subroutine RadC2001Init

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

end module rad_C2001

