subroutine RadC2001CalcCloudOptProp( Spec, iband, xyz_REff, xyz_ExtCoef, xyz_SSA, xyz_AF )

    ! 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 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 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 RadC2001CalcTrans( Spec, iband, xyz_DelAbsMass, xyz_Press, xyz_Temp, xyz_QMix, xyrr_Trans )

    ! 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, xyz_DelAbsMassScaled )
    else
      call RadC2001LineScaleAmt( xyz_DelAbsMass, xyz_Press, xyz_Temp, RefPress, RefTemp, PressScaleIndex, CoefA, CoefB, xyz_DelAbsMassScaled )
    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 RadC2001CalcTransBand3CO2( xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2, xyrr_Trans )

    ! 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, xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2, xyrr_TransEach )

      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, xyrr_Trans )

    ! 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, xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransCont )
      Spec = 'SubBandH2OLine'
      call RadC2001CalcTrans( Spec, iband, xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransLine )

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

  subroutine RadC2001ReduceCloudOptDep( xyz_SSA, xyz_AF, xyz_DelOptDep )

    ! 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 RadC2001ContScaleH2OAmt( xyz_DelAbsMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyz_ContDelAbsMassScaled )

    ! 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 * exp( 1800.0d0 * ( 1.0d0 / xyz_Temp - 1.0d0 / 296.0d0 ) )


  end subroutine RadC2001ContScaleH2OAmt

  subroutine RadC2001LineScaleAmt( xyz_DelAbsMass, xyz_Press, xyz_Temp, RefPress, RefTemp, PressScaleIndex, CoefA, CoefB, xyz_DelAbsMassScaled )

    ! 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