subroutine RadCK1991CalcTrans( xyz_DelAbsMass, xyz_Press, xyz_Temp, Spec, xyrr_Trans )

    ! 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)
    character(len=*), intent(in ) :: Spec
    real(DP)        , intent(out) :: xyrr_Trans    (0:imax-1, 1:jmax, 0:kmax, 0:kmax)


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


    if ( Spec == 'CO2' ) then
      call RadCK1991CalcTransCore( NTabPress, NCO2TabAbsAmt, a_TabLog10Press, a_CO2TabLog10AbsAmt, aa_CO2TabAlpha, aa_CO2TabBeta, aa_CO2TabAbs, xyz_DelAbsMass, xyz_Press, xyz_Temp, xyrr_Trans )
    else if ( Spec == 'O3' ) then
      call RadCK1991CalcTransCore( NTabPress, NO3TabAbsAmt, a_TabLog10Press, a_O3TabLog10AbsAmt, aa_O3TabAlpha, aa_O3TabBeta, aa_O3TabAbs, xyz_DelAbsMass, xyz_Press, xyz_Temp, xyrr_Trans )
    else
      call MessageNotify( 'E', module_name, 'Specified composition, %c, is inappropriate', c1 = trim(Spec) )
    end if


  end subroutine RadCK1991CalcTrans

  subroutine RadCK1991Init


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

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


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


!!$    namelist /rad_RL78_nml/ &
!!$      & VMRCO2,                   &
!!$      & DelTimeCalcTransValue,    &
!!$      & DelTimeCalcTransUnit,     &
!!$      & flag_save_time


    if ( rad_CK1991_inited ) return

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


    ! Convert unit of pressure from mbar to Pa
    !
    a_TabLog10Press = 10.0d0**a_TabLog10Press
    a_TabLog10Press = a_TabLog10Press * 1.0d2
    a_TabLog10Press = log10( a_TabLog10Press )

    ! Convert unit of absorber amount from (atm cm)_{STP} to kg m-2
    !   To convert from (atm cm)_{STP} to kg m-2, the value is divided by
    !   1.0d2 / 101325.0d0 * 8.31432d0 / ( 44.0d-3 ) * 273.15d0, and
    !   1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 273.15d0, 
    !   for CO2 and O3, respectively.
    !   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.
    !
    a_CO2TabLog10AbsAmt = 10.0d0**a_CO2TabLog10AbsAmt
    a_CO2TabLog10AbsAmt = a_CO2TabLog10AbsAmt / ( 1.0d2 / 101325.0d0 * 8.31432d0 / ( 44.0d-3 ) * 273.15d0 )
    a_CO2TabLog10AbsAmt = log10( a_CO2TabLog10AbsAmt )
    !
    a_O3TabLog10AbsAmt  = 10.0d0**a_O3TabLog10AbsAmt
    a_O3TabLog10AbsAmt  = a_O3TabLog10AbsAmt / ( 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 273.15d0 )
    a_O3TabLog10AbsAmt  = log10( a_O3TabLog10AbsAmt )
    !
    ! Convert values absorptance from -1e3 * log10(A) to log10(A)
    !
    aa_CO2TabAbs = 10.0d0**( aa_CO2TabAbs / ( -1.0d3 ) )
    aa_CO2TabAbs = log10( aa_CO2TabAbs )
    !
    aa_O3TabAbs  = 10.0d0**( aa_O3TabAbs / ( -1.0d3 ) )
    aa_O3TabAbs  = log10( aa_O3TabAbs )
    !
    ! Convert values of alpha absorptance from 1e4 * alpha to alpha
    !
    aa_CO2TabAlpha = aa_CO2TabAlpha / ( 1.0d4 )
    !
    aa_O3TabAlpha  = aa_O3TabAlpha  / ( 1.0d4 )
    !
    ! Convert values of beta absorptance from 1e6 * beta to beta
    !
    aa_CO2TabBeta  = aa_CO2TabBeta  / ( 1.0d6 )
    !
    aa_O3TabBeta   = aa_O3TabBeta   / ( 1.0d6 )



    ! 印字 ; 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_ck1991_inited = .true.

  end subroutine RadCK1991Init

  subroutine RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_Tab, xy_IndexPress, xy_IndexAbsAmt, xy_Log10Press, xy_Log10AbsAmt, xy_Array )

    ! USE statements
    !


    integer , intent(in ) :: NTabPress
    integer , intent(in ) :: NTabAbsAmt
    real(DP), intent(in ) :: a_TabLog10Press (NTabPress )
    real(DP), intent(in ) :: a_TabLog10AbsAmt(NTabAbsAmt)
    real(DP), intent(in ) :: aa_Tab          (NTabAbsAmt, NTabPress)
    integer , intent(in ) :: xy_IndexPress   (0:imax-1, 1:jmax)
    integer , intent(in ) :: xy_IndexAbsAmt  (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_Log10Press   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_Log10AbsAmt  (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xy_Array        (0:imax-1, 1:jmax)

    !
    ! Work variables
    !
    real(DP) :: val1
    real(DP) :: val2
    integer  :: ip1
    integer  :: ip2
    integer  :: iw1
    integer  :: iw2
    integer  :: i
    integer  :: j


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


    do j = 1, jmax
      do i = 0, imax-1

        ip1 = xy_IndexPress (i,j) - 1
        ip2 = xy_IndexPress (i,j)
        iw1 = xy_IndexAbsAmt(i,j) - 1
        iw2 = xy_IndexAbsAmt(i,j)

        val1 = ( aa_Tab          (iw2,ip1) - aa_Tab          (iw1,ip1) ) / ( a_TabLog10AbsAmt(iw2)     - a_TabLog10AbsAmt(iw1)     ) * ( xy_Log10AbsAmt  (i,j)     - a_TabLog10AbsAmt(iw1)     ) + aa_Tab            (iw1,ip1)
        val2 = ( aa_Tab          (iw2,ip2) - aa_Tab          (iw1,ip2) ) / ( a_TabLog10AbsAmt(iw2)     - a_TabLog10AbsAmt(iw1)     ) * ( xy_Log10AbsAmt  (i,j)     - a_TabLog10AbsAmt(iw1)     ) + aa_Tab            (iw1,ip2)

        xy_Array(i,j) = ( val2                 - val1                 ) / ( a_TabLog10Press(ip2) - a_TabLog10Press(ip1) ) * ( xy_Log10Press  (i,j) - a_TabLog10Press(ip1) ) + val1

      end do
    end do

  end subroutine RadCK1991Interpolate

  subroutine RadCK1991CalcTransCore( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabAlpha, aa_TabBeta, aa_TabAbs, xyz_DelAbsMass, xyz_Press, xyz_Temp, xyrr_Trans )

    ! USE statements
    !


    integer , intent(in ) :: NTabPress
    integer , intent(in ) :: NTabAbsAmt
    real(DP), intent(in ) :: a_TabLog10Press (NTabPress )
    real(DP), intent(in ) :: a_TabLog10AbsAmt(NTabAbsAmt)
    real(DP), intent(in ) :: aa_TabAlpha     (NTabAbsAmt, NTabPress)
    real(DP), intent(in ) :: aa_TabBeta      (NTabAbsAmt, NTabPress)
    real(DP), intent(in ) :: aa_TabAbs       (NTabAbsAmt, NTabPress)
    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(out) :: xyrr_Trans      (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    !
    ! Work variables
    !
    real(DP), parameter :: RefTemp = 250.0_DP

    real(DP) :: xy_EffTemp      (0:imax-1, 1:jmax)
    real(DP) :: xy_EffPress     (0:imax-1, 1:jmax)
    real(DP) :: xy_AbsAmt       (0:imax-1, 1:jmax)
    real(DP) :: xy_Log10EffPress(0:imax-1, 1:jmax)
    real(DP) :: xy_Log10AbsAmt  (0:imax-1, 1:jmax)
    integer  :: xy_IndexPress   (0:imax-1, 1:jmax)
    integer  :: xy_IndexAbsAmt  (0:imax-1, 1:jmax)

    real(DP) :: xy_Alpha        (0:imax-1, 1:jmax)
    real(DP) :: xy_Beta         (0:imax-1, 1:jmax)
    real(DP) :: xy_AbsAtRefTemp (0:imax-1, 1:jmax)

    real(DP) :: xy_Abs          (0:imax-1, 1:jmax)

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


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


    do k = 0, kmax
      do kk = k+1, kmax

        xy_EffTemp  = 0.0_DP
        xy_EffPress = 0.0_DP
        xy_AbsAmt   = 0.0_DP
        do l = k+1, kk
          xy_EffTemp  = xy_EffTemp  + xyz_Temp (:,:,l) * xyz_DelAbsMass(:,:,l)
          xy_EffPress = xy_EffPress + xyz_Press(:,:,l) * xyz_DelAbsMass(:,:,l)
          xy_AbsAmt   = xy_AbsAmt   +                    xyz_DelAbsMass(:,:,l)
        end do
        xy_EffTemp  = xy_EffTemp  / ( xy_AbsAmt + 1.0d-100 )
        xy_EffPress = xy_EffPress / ( xy_AbsAmt + 1.0d-100 )


        xy_Log10EffPress = log10( xy_EffPress + 1.0d-100 )
        xy_Log10AbsAmt   = log10( xy_AbsAmt   + 1.0d-100 )


        do j = 1, jmax
          do i = 0, imax-1

            if ( xy_Log10EffPress(i,j) < a_TabLog10Press (1) ) then
!!$              call MessageNotify( 'E', module_name, &
!!$                & 'at k = %d and kk = %d, Log10EffPress(%d,%d) = %f < %f. too small', &
!!$                & i = (/k, kk, i, j/), &
!!$                & d = (/xy_Log10EffPress(i,j), a_TabLog10Press(1)/) )

              xy_IndexPress(i,j) = 2

            else

#ifdef VECTOR
              xy_IndexPress(i,j) = 1 + 1
              Loop_Search_Press  : do l = 1+1, NTabPress-1
                if ( a_TabLog10Press (l) < xy_Log10EffPress(i,j) ) then
                  xy_IndexPress(i,j) = l + 1
                end if
              end do Loop_Search_Press
#else
              Loop_Search_Press  : do l = 1+1, NTabPress
                if ( a_TabLog10Press (l) > xy_Log10EffPress(i,j) ) then
                  exit Loop_Search_Press
                end if
              end do Loop_Search_Press
              if ( l > NTabPress ) then
                l = NTabPress

!!$                call MessageNotify( 'E', module_name, &
!!$                  & 'at k = %d and kk = %d, Log10EffPress(%d,%d) = %f > %f. too large', &
!!$                  & i = (/k, kk, i, j/), &
!!$                  & d = (/xy_Log10EffPress(i,j), a_TabLog10Press(NTabPress)/) )

              end if
              xy_IndexPress(i,j) = l
#endif

            end if


            if ( xy_Log10AbsAmt(i,j) < a_TabLog10AbsAmt(1) ) then
!!$              call MessageNotify( 'E', module_name, &
!!$                & 'at k = %d and kk = %d, Log10AbsAmt(%d,%d) = %f < %f. too small', &
!!$                & i = (/k, kk, i, j/), &
!!$                & d = (/xy_Log10AbsAmt(i,j), a_TabLog10AbsAmt(1)/) )

              xy_IndexAbsAmt(i,j) = 2
            else

#ifdef VECTOR
              xy_IndexAbsAmt(i,j) = 1 + 1
              Loop_Search_AbsAmt : do l = 1+1, NTabAbsAmt-1
                if ( a_TabLog10AbsAmt(l) < xy_Log10AbsAmt(i,j) ) then
                  xy_IndexAbsAmt(i,j) = l + 1
                end if
              end do Loop_Search_AbsAmt
#else
              Loop_Search_AbsAmt : do l = 1+1, NTabAbsAmt
                if ( a_TabLog10AbsAmt(l) > xy_Log10AbsAmt(i,j) ) then
                  exit Loop_Search_AbsAmt
                end if
              end do Loop_Search_AbsAmt
              if ( l > NTabAbsAmt ) then

                l = NTabAbsAmt
!!$                call MessageNotify( 'E', module_name, &
!!$                  & 'at k = %d and kk = %d, Log10AbsAmt(%d,%d) = %f > %f. too large', &
!!$                  & i = (/k, kk, i, j/), &
!!$                  & d = (/xy_Log10AbsAmt(i,j), a_TabLog10AbsAmt(NTabAbsAmt)/) )
              end if
              xy_IndexAbsAmt(i,j) = l
#endif

            end if

          end do
        end do


        call RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabAlpha, xy_IndexPress, xy_IndexAbsAmt, xy_Log10EffPress, xy_Log10AbsAmt, xy_Alpha )
        call RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabBeta, xy_IndexPress, xy_IndexAbsAmt, xy_Log10EffPress, xy_Log10AbsAmt, xy_Beta )
        call RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabAbs, xy_IndexPress, xy_IndexAbsAmt, xy_Log10EffPress, xy_Log10AbsAmt, xy_AbsAtRefTemp )

        xy_AbsAtRefTemp = 10.0d0**xy_AbsAtRefTemp
        xy_Abs = xy_AbsAtRefTemp * ( 1.0_DP + xy_Alpha * ( xy_EffTemp - RefTemp ) + xy_Beta  * ( xy_EffTemp - RefTemp )**2 )

        xyrr_Trans(:,:,k,kk) = 1.0_DP - xy_Abs
      end do
    end do
    !
    ! correction
    !
    do k = 0, kmax
      do kk = k+2, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyrr_Trans(i,j,k,kk) > xyrr_Trans(i,j,k,kk-1) ) then
              xyrr_Trans(i,j,k,kk) = xyrr_Trans(i,j,k,kk-1)
            end if
          end do
        end do
      end do
    end do

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


  end subroutine RadCK1991CalcTransCore