subroutine RadMars15m( Time, DelTime, xyz_Temp, xyr_Press, xyz_Press, xy_SurfTemp, xyr_DOD067, QeRat, SSA, xy_SurfEmis, xyr_Rad15mFlux, xyra_DelRad15mFlux )

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

    real(DP)    , intent(in ) :: Time
    real(DP)    , intent(in ) :: DelTime
    real(DP)    , intent(in ) :: xyz_Temp          (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_Rad15mFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(out) :: xyra_DelRad15mFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP)    , intent(in ) :: xyr_DOD067        (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(in ) :: QeRat
    real(DP)    , intent(in ) :: SSA
    real(DP)    , intent(in ) :: xy_SurfEmis       (0:imax-1, 1:jmax)


    !
    ! local variables
    !


    ! 実行文 ; Executable statement
    !

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

    call rad15m_lowatm_newscheme2006( Time, DelTime, xyz_Temp, xyr_Press, xyz_Press, xy_SurfTemp, xyr_DOD067, QeRat, SSA, xy_SurfEmis, xyr_Rad15mFlux, xyra_DelRad15mFlux )

  end subroutine RadMars15m

  subroutine RadMars15mInit

    ! モジュール引用 ; USE statements
    !

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

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

    use ckd_module, only : ckd_input, ckdp, nband


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

    character(STRING) :: rad15mkg_fn
    character(STRING) :: rad15mnf_fn

    integer           :: m

    namelist /rad_Mars_15m_nml/ rad15mkg_fn, Rad15mInt


    ! 実行文 ; Executable statement
    !

    if ( rad_Mars_15m_inited ) return

    ! デフォルト値の設定
    ! Default values settings
    !

    rad15mkg_fn = "./kg15m"
!!$    rad15mnf_fn = "./nlte15mfactor"
    Rad15mInt   = 925.0_DP

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = rad_Mars_15m_nml, iostat = iostat_nml )         ! (out)
      close( unit_nml )

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


!!$    allocate( rad_gp   ( im, jm,   km ) )
!!$    allocate( rad_gph  ( im, jm, 0:km ) )
!!$    allocate( rad_gt   ( im, jm,   km ) )
!!$    allocate( rad_gts  ( im, jm,   1  ) )
!!$    allocate( rad_gdod ( im, jm, 0:km ) )


    nras = 1
    nrps = 0

!!$    allocate( sgmh_f          ( km*nvr+1 ), &
!!$      &       sgm_f           ( km*nvr   ) )
!!$    allocate( gph_f    ( im, jm, km*nvr+1 ), &
!!$      &       gp_f     ( im, jm, km*nvr   ), &
!!$      &       gth_f    ( im, jm, km*nvr+1 ) )
!!$
!!$    allocate( gvmr_f   ( im, jm, km*nvr  , nras + nrps ) )
!!$    allocate( mmmass_f ( im, jm, km*nvr   ) )
!!$    allocate( ac_f     ( im, jm, km*nvr  , nras        ) )
!!$
!!$    allocate( gdod_f   ( im, jm, km*nvr+1 ) )
!!$
!!$
    allocate( xyra_Trans  (0:imax-1, 1:jmax, 0:kmax, 0:kmax) )
!!$    allocate( pfh_f    ( im, jm, km*nvr+1 ) )
!!$
!!$    allocate( uwflh_f  ( im, jm, km*nvr+1 ), &
!!$      &       dwflh_f  ( im, jm, km*nvr+1 ) )


    allocate( trans_i2i_toa(0:imax-1, 1:jmax, 0:kmax), trans_i2i_boa(0:imax-1, 1:jmax, 0:kmax), trans_i2i_s  (0:imax-1, 1:jmax, 0:kmax), trans_i2m_lli(0:imax-1, 1:jmax, 0:kmax, 1:kmax), trans_i2m_uli(0:imax-1, 1:jmax, 0:kmax, 1:kmax) )

    trans_i2i_toa(:,:,:)   = 1.0d100
    trans_i2i_boa(:,:,:)   = 1.0d100
    trans_i2i_s  (:,:,:)   = 1.0d100
    trans_i2m_lli(:,:,:,:) = 1.0d100
    trans_i2m_uli(:,:,:,:) = 1.0d100


    !
    ! check
    !
    if( nras .ne. 1 ) then
      write( 6, * ) 'nras is not 1.'
      write( 6, * ) nras
      stop
    end if

    call ckd_input( rad15mkg_fn )

    ! check
    if( nband /= 1 ) then
      write( 6, * ) ' nband is not 1.'
      write( 6, * ) nband
      stop
    end if

    nwnl = 0
    do m = 1, nband
      nwnl = nwnl + ckdp( m ) % ng
    end do



!!$    call increase_vreso_boundary( km, nvr, sgmh, sgmh_f, "log" )
!!$    do k = 1, km * nvr
!!$      sgm_f( k ) = sqrt( sgmh_f( k ) * sgmh_f( k+1 ) )
!!$    end do


!!$    call rad15m_readnlte15mfac( rad15mnf_fn )


    !
    ! This routine must be called after rad15m_readkgtbl.
    !
!!$      call rad15m_rv_read( time )
!!$    call rad15m_rv_read_newscheme2006( time )


    rad_Mars_15m_inited = .true.


  end subroutine RadMars15mInit

  subroutine calc_lnp( xyz_Press, xyz_lnPress )

    real(DP), intent(in ) :: xyz_Press  (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !


    xyz_lnPress = log( xyz_Press + 1.0d-20 )


  end subroutine calc_lnp

  subroutine calc_rteq_use_meantrans_arr3d( dlambda, xy_SurfEmis, trans_i2i_toa, trans_i2i_boa, trans_i2i_s, trans_i2m_lli, trans_i2m_uli, xyr_PF, xy_SurfPF, netflh )

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

    real(DP), intent(in ) :: dlambda
    real(DP), intent(in ) :: xy_SurfEmis(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: trans_i2i_toa(0:imax-1, 1:jmax, 0:kmax)    ! f_{1/2}    T_{k+1/2,1/2}
    real(DP), intent(in ) :: trans_i2i_boa(0:imax-1, 1:jmax, 0:kmax)    ! f_{km+1/2} T_{k+1/2,km+1/2}
    real(DP), intent(in ) :: trans_i2i_s  (0:imax-1, 1:jmax, 0:kmax)    ! f_{s}      T_{k+1/2,km+1/2}
    real(DP), intent(in ) :: trans_i2m_lli(0:imax-1, 1:jmax, 0:kmax, 1:kmax) ! upper layer interface
    real(DP), intent(in ) :: trans_i2m_uli(0:imax-1, 1:jmax, 0:kmax, 1:kmax) ! lower layer interface
    real(DP), intent(in ) :: xyr_PF   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfPF(0:imax-1, 1:jmax)
    real(DP), intent(out) :: netflh(0:imax-1, 1:jmax, 0:kmax)



    !
    ! local variables
    !
    integer :: i, j, k, k2


    do k = 0, kmax
      do j = 1, jmax
        do i = 0, imax-1
          netflh(i,j,k) = 0.0d0
        end do
      end do
    end do

    do k = 0, kmax

      do j = 1, jmax
        do i = 0, imax-1
          netflh(i,j,k) = netflh(i,j,k) + PI * xy_SurfEmis(i,j) * xy_SurfPF(i,j) * dlambda * trans_i2i_s  (i,j,k) - PI * xyr_PF(i,j,0   ) * dlambda * trans_i2i_boa(i,j,k) + PI * xyr_PF(i,j,kmax) * dlambda * trans_i2i_toa(i,j,k)
        end do
      end do

      do k2 = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            netflh(i,j,k) = netflh(i,j,k) - PI * xyr_PF(i,j,k2  ) * dlambda * trans_i2m_uli(i,j,k,k2) + PI * xyr_PF(i,j,k2-1) * dlambda * trans_i2m_lli(i,j,k,k2)
          end do
        end do
      end do

    end do


  end subroutine calc_rteq_use_meantrans_arr3d

  subroutine calc_trans_mp_arr3d( nras, nrps, xyr_Press, xyza_VMR, xyz_MMMass, xyza_AC, xyr_DOD, xyra_Trans )

    use constants , only : Grav

    integer , intent(in ) :: nras
    integer , intent(in ) :: nrps
    real(DP), intent(in ) :: xyr_Press(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xyza_VMR (0:imax-1, 1:jmax, 1:kmax, 1:nras+nrps)
    real(DP), intent(in ) :: xyz_MMMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyza_AC(0:imax-1, 1:jmax, 1:kmax, 1:nras)
    real(DP), intent(in ) :: xyr_DOD(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyra_Trans  (0:imax-1, 1:jmax, 0:kmax, 0:kmax)


    !
    ! local variables
    !
    real(DP)     :: xyz_DelOpDep(0:imax-1, 1:jmax, 1:kmax)
    real(DP)     :: xyz_DelTrans(0:imax-1, 1:jmax, 1:kmax)
    real(DP)     :: xy_Trans(0:imax-1, 1:jmax)
    real(DP), parameter :: DifFac = 1.66_DP


    integer :: k, k2, n
    integer :: ks, ke



    xyra_Trans = 1.0d100


    xyz_DelOpDep = 0.0_DP
    do n = 1, nras
      do k = 1, kmax
        xyz_DelOpDep(:,:,k) = xyz_DelOpDep(:,:,k) + xyza_AC(:,:,k,n) * xyza_VMR(:,:,k,n) / xyz_MMMass(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
      end do
    end do

    !
    ! add dust optical depth
    !
    do k = 1, kmax
      xyz_DelOpDep(:,:,k) = xyz_DelOpDep(:,:,k) + xyr_DOD(:,:,k-1) - xyr_DOD(:,:,k)
    end do

    xyz_DelTrans = exp( - xyz_DelOpDep * DifFac )


    !
    ! transmission for "zero thickness" layer ( = 1.0 )
    !
    do ks = 0, kmax
      ke = ks
      xyra_Trans(:,:,ks,ke) = 1.0_DP
    end do

    do ks = 0, kmax
      xy_Trans = 1.0_DP
      do ke = ks+1, kmax
        xy_Trans = xy_Trans * xyz_DelTrans(:,:,ke)
        xyra_Trans(:,:,ks,ke) = xy_Trans
      end do
    end do

    do ks = 0, kmax
      do ke = 0, ks-1
        xyra_Trans(:,:,ks,ke) = xyra_Trans(:,:,ke,ks)
      end do
    end do


  end subroutine calc_trans_mp_arr3d

  subroutine findindices( ks, ke, xyz_Temp, xyz_lnPress, iband, xyz_jj, xyz_kk )

    use ckd_module, only : ckdp

    integer , intent(in ) :: ks
    integer , intent(in ) :: ke
    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, ks:ke)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: iband
    integer , intent(out) :: xyz_jj(0:imax-1, 1:jmax, ks:ke)
    integer , intent(out) :: xyz_kk(0:imax-1, 1:jmax, ks:ke)


    !
    ! local variables
    !
    integer :: i, j, k, l


    do k = ks, ke
      do j = 1, jmax
        do i = 0, imax-1
          xyz_kk(i,j,k) = 1
        end do
      end do
    end do

    do l = 1+1, ckdp( iband ) % nt - 1
      do k = ks, ke
        do j = 1, jmax
          do i = 0, imax-1
            if( ckdp( iband ) % t( l ) .le. xyz_Temp(i,j,k) ) xyz_kk(i,j,k) = l
          end do
        end do
      end do
    end do

    do k = ks, ke
      do j = 1, jmax
        do i = 0, imax-1
          xyz_jj(i,j,k) = 1
        end do
      end do
    end do
    do l = 1+1, ckdp( iband ) % nlnp - 1
      do k = ks, ke
        do j = 1, jmax
          do i = 0, imax-1
            if( ckdp( iband ) % lnp( l ) <= xyz_lnPress(i,j,k) ) xyz_jj(i,j,k) = l
          end do
        end do
      end do
    end do


  end subroutine findindices

  subroutine findindices2D( xy_Temp, xy_lnPress, iband, xy_jj, xy_kk )

    real(DP), intent(in ) :: xy_Temp   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_lnPress(0:imax-1, 1:jmax)
    integer , intent(in ) :: iband
    integer , intent(out) :: xy_jj(0:imax-1, 1:jmax)
    integer , intent(out) :: xy_kk(0:imax-1, 1:jmax)


    !
    ! local variables
    !
    real(DP) :: xyz_Temp   (0:imax-1, 1:jmax, 1:1)
    real(DP) :: xyz_lnPress(0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_jj(0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_kk(0:imax-1, 1:jmax, 1:1)


    xyz_Temp   (:,:,1) = xy_Temp
    xyz_lnPress(:,:,1) = xy_lnPress

    call findindices( 1, 1, xyz_Temp, xyz_lnPress, iband, xyz_jj, xyz_kk )

    xy_jj = xyz_jj(:,:,1)
    xy_kk = xyz_kk(:,:,1)


  end subroutine findindices2D

  subroutine findindices3D( xyz_Temp, xyz_lnPress, iband, xyz_jj, xyz_kk )

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: iband
    integer , intent(out) :: xyz_jj(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(out) :: xyz_kk(0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !


    call findindices( 1, kmax, xyz_Temp, xyz_lnPress, iband, xyz_jj, xyz_kk )


  end subroutine findindices3D

  subroutine getlnac_givenindices( xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, xyz_AC )

    use ckd_module, only : ckdp

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_jj  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_kk  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: ig
    integer , intent(in ) :: iband
    real(DP), intent(out) :: xyz_AC  (0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !
    real(DP) :: lnac1, lnac2
    integer  :: i, j, k


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

          lnac1 = ( ckdp(iband)%lnac( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)+1 ) - ckdp(iband)%lnac( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   ) ) / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 ) - ckdp(iband)%t( xyz_kk(i,j,k)   ) ) * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) ) + ckdp(iband)%lnac( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   )
          lnac2 = ( ckdp(iband)%lnac( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)+1 ) - ckdp(iband)%lnac( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   ) ) / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 ) - ckdp(iband)%t( xyz_kk(i,j,k)   ) ) * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) ) + ckdp(iband)%lnac( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   )

          xyz_AC(i,j,k) = ( lnac2 - lnac1 ) / ( ckdp( iband ) % lnp( xyz_jj(i,j,k)+1 ) - ckdp( iband ) % lnp( xyz_jj(i,j,k)   ) ) * ( xyz_lnPress(i,j,k) - ckdp( iband ) % lnp( xyz_jj(i,j,k) ) ) + lnac1
        end do
      end do
    end do


  end subroutine getlnac_givenindices

  subroutine getpf_arr3d_norat( xyz_Temp, xy_SurfTemp, iband, xyr_PF, xy_SurfPF )

    use planck_func, only : Integ_PF_GQ_Array3D, Integ_PF_GQ_Array2D

    use ckd_module, only : ckdp

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfTemp(0:imax-1, 1:jmax)
    integer , intent(in ) :: iband
    real(DP), intent(out) :: xyr_PF   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xy_SurfPF(0:imax-1, 1:jmax)


    !
    ! local variables
    !
    integer :: ncp_pfint
    integer :: i, j, k


    ncp_pfint = 5

    call Integ_PF_GQ_Array3D( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), ncp_pfint, 0, imax-1, 1, jmax, 0, kmax, xyz_Temp, xyr_PF )
    call Integ_PF_GQ_Array2D( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), ncp_pfint, 0, imax-1, 1, jmax, xy_SurfTemp, xy_SurfPF )

    do k = 0, kmax
      do j = 1, jmax
        do i = 0, imax-1
          xyr_PF(i,j,k) = xyr_PF(i,j,k) / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )
        end do
      end do
    end do
    do j = 1, jmax
      do i = 0, imax-1
        xy_SurfPF(i,j) = xy_SurfPF(i,j) / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )
      end do
    end do


  end subroutine getpf_arr3d_norat

  subroutine getpfr_givenindices( ks, ke, xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, xyz_PFRat )

    use ckd_module, only: ckdp

    integer , intent(in ) :: ks
    integer , intent(in ) :: ke
    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, ks:ke)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: xyz_jj  (0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: xyz_kk  (0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: ig, iband
    real(DP), intent(out) :: xyz_PFRat(0:imax-1, 1:jmax, ks:ke)


    !
    ! local variables
    !
    real(DP) :: pfr1, pfr2
    integer  :: i, j, k, l


    do k = ks, ke
      do j = 1, jmax
        do i = 0, imax-1

          pfr1 = ( ckdp(iband)%pfr( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)+1 ) - ckdp(iband)%pfr( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   ) ) / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 ) - ckdp(iband)%t( xyz_kk(i,j,k)   ) ) * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) ) + ckdp(iband)%pfr( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   )
          pfr2 = ( ckdp(iband)%pfr( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)+1 ) - ckdp(iband)%pfr( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   ) ) / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 ) - ckdp(iband)%t( xyz_kk(i,j,k)   ) ) * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) ) + ckdp(iband)%pfr( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   )


          xyz_PFRat(i,j,k) = ( pfr2 - pfr1 ) / ( ckdp( iband ) % lnp( xyz_jj(i,j,k)+1 ) - ckdp( iband ) % lnp( xyz_jj(i,j,k)   ) ) * ( xyz_lnPress(i,j,k) - ckdp( iband ) % lnp( xyz_jj(i,j,k) ) ) + pfr1
        end do
      end do
    end do


  end subroutine getpfr_givenindices

  subroutine getpfr_givenindices2D( xy_Temp, xy_lnPress, xy_jj, xy_kk, ig, iband, xy_PFRat )

    real(DP), intent(in ) :: xy_Temp   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_lnPress(0:imax-1, 1:jmax)
    integer , intent(in ) :: xy_jj  (0:imax-1, 1:jmax)
    integer , intent(in ) :: xy_kk  (0:imax-1, 1:jmax)
    integer , intent(in ) :: ig, iband
    real(DP), intent(out) :: xy_PFRat (0:imax-1, 1:jmax)


    !
    ! local variables
    !
    real(DP) :: xyz_Temp   (0:imax-1, 1:jmax, 1:1)
    real(DP) :: xyz_lnPress(0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_jj  (0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_kk  (0:imax-1, 1:jmax, 1:1)
    real(DP) :: xyz_PFRat(0:imax-1, 1:jmax, 1:1)


    xyz_Temp   (:,:,1) = xy_Temp
    xyz_lnPress(:,:,1) = xy_lnPress
    xyz_jj  (:,:,1) = xy_jj
    xyz_kk  (:,:,1) = xy_kk

    call getpfr_givenindices( 1, 1, xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, xyz_PFRat )

    xy_PFRat(:,:) = xyz_PFRat(:,:,1)


  end subroutine getpfr_givenindices2D

  subroutine getpfr_givenindices3D( xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, xyz_PFRat )

    use ckd_module, only: ckdp

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_jj  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_kk  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: ig, iband
    real(DP), intent(out) :: xyz_PFRat(0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !

    call getpfr_givenindices( 1, kmax, xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, xyz_PFRat )


  end subroutine getpfr_givenindices3D

  subroutine m2ckdpindices( m, ig, iband )

    use ckd_module, only : ckdp, nband

    integer, intent(in ) :: m
    integer, intent(out) :: ig
    integer, intent(out) :: iband


    !
    ! local variables
    !
    integer :: num


    ! The comments below will be removed.


    num = 0
    do iband = 1, nband
      if( num + ckdp( iband ) % ng .ge. m ) exit
      num = num + ckdp( iband ) % ng
    end do
    if( iband > nband ) then
      write( 6, * ) 'Unexpected m'
      write( 6, * ) m
      stop
    end if
    ig = m - num
    if( ig > ckdp( iband ) % ng ) then
      write( 6, * ) 'Unexpected ig'
      write( 6, * ) iband, ig
      stop
    end if


  end subroutine m2ckdpindices

  subroutine rad15m_lowatm_newscheme2006( Time, DelTime, xyz_Temp, xyr_Press, xyz_Press, xy_SurfTemp, xyr_DOD067, QeRat, SSA, xy_SurfEmis, xyr_Rad15mFlux, xyra_DelRad15mFlux )


    use constants , only : Grav, CpDry

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

    use ckd_module, only : ckdp


    real(DP)    , intent(in ) :: Time
    real(DP)    , intent(in ) :: DelTime
    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_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in ) :: xy_SurfTemp(0:imax-1, 1:jmax)
    real(DP)    , intent(in ) :: xyr_DOD067(0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(in ) :: QeRat
    real(DP)    , intent(in ) :: SSA
    real(DP)    , intent(in ) :: xy_SurfEmis(0:imax-1, 1:jmax)

    real(DP)    , intent(out) :: xyr_Rad15mFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(out) :: xyra_DelRad15mFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)

    !
    ! local variables
    !
    real(DP) :: xyr_Temp  (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyz_MMMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyza_VMR(0:imax-1, 1:jmax, 1:kmax, 1:nras+nrps )
    real(DP) :: xyza_AC (0:imax-1, 1:jmax, 1:kmax, 1:nras      )
    real(DP) :: xyr_PF   (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xy_SurfPF(0:imax-1, 1:jmax)

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

    real(DP) :: weight_integral
    integer  :: ig, iband

    integer  :: i, j, k, l, m, n
    integer  :: k2

    !
    ! dod      : dust optical depth
    !
    real(DP) :: xyr_DOD(0:imax-1, 1:jmax, 0:kmax)

    !
    ! local variables for pfint
    !

    real(DP) :: MinPress
    real(DP) :: MaxPress

    integer  :: iband_reserve
    real(DP) :: xy_lnPs    (0:imax-1, 1:jmax)
    real(DP) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax )
    integer  :: xyz_jj(0:imax-1, 1:jmax, 1:kmax)
    integer  :: xyz_kk(0:imax-1, 1:jmax, 1:kmax)
    integer  :: xy_jj (0:imax-1, 1:jmax)
    integer  :: xy_kk (0:imax-1, 1:jmax)


    ! Surface temperature for calculation of gradient of radiative flux
    real(DP) :: xy_SurfTemp_for_gradcalc(0:imax-1, 1:jmax)
    ! Indices for calculation of gradient of radiative flux
    integer  :: jjs_for_gradcalc(0:imax-1, 1:jmax), kks_for_gradcalc(0:imax-1, 1:jmax)

    real(DP)     :: xyr_PFRat   (0:imax-1, 1:jmax, 0:kmax)
    real(DP)     :: xyz_PFRat   (0:imax-1, 1:jmax, 1:kmax)
    real(DP)     :: xy_SurfPFRat(0:imax-1, 1:jmax)

    logical, save :: FlagCalcTrans

    data FlagCalcTrans / .false. /


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


    k = 0
    do j = 1, jmax
      do i = 0, imax-1
!!$        gth(i,j,k) = gt(i,j,k+1)
        xyr_Temp(i,j,k) = ( xyz_Temp(i,j,2) - xyz_Temp(i,j,1) ) / log( xyz_Press(i,j,2) / xyz_Press(i,j,1) ) * log( xyr_Press(i,j,k) / xyz_Press(i,j,1) ) + xyz_Temp(i,j,1)
      end do
    end do
    do k = 1, kmax-1
      do j = 1, jmax
        do i = 0, imax-1
          xyr_Temp(i,j,k) = ( xyz_Temp(i,j,k+1) - xyz_Temp(i,j,k) ) / log( xyz_Press(i,j,k+1) / xyz_Press(i,j,k) ) * log( xyr_Press(i,j,k  ) / xyz_Press(i,j,k) ) + xyz_Temp(i,j,k)
        end do
      end do
    end do
    k = kmax
    do j = 1, jmax
      do i = 0, imax-1
        xyr_Temp(i,j,k) = xyz_Temp(i,j,k)
      end do
    end do


!!$    do k = 1, km*nvr+1
!!$      do ij = ijs, ije
!!$        gph_f( ij, 1, k ) = sgmh_f( k ) * gph( ij, 1, km+1 )
!!$      end do
!!$    end do
!!$    call calc_lnp( im, jm, km*nvr+1, gph_f , glnph_f , ijs, ije )
!!$
!!$    call increase_vreso_boundary_arr3d( im, jm, km, nvr, gth, gth_f, &
!!$      & "linear", ijs, ije )



    if (  .not. FlagCalcTrans ) then
      if ( Time - dble( int( Time / Rad15mInt ) ) * Rad15mInt < DelTime ) then
        call MessageNotify( 'M', module_name, 'Transmittance is not saved, but criterion for transmittance calculation is met.' )
      else
        call MessageNotify( 'M', module_name, 'Transmittance is not saved, and criterion for transmittance calculation ' // 'is not met. However, transmittance will be calculated.' )
      end if
    end if


    !
    ! Calculation of transmission
    !
    if( ( .not. FlagCalcTrans ) .or. ( Time - dble( int( Time / Rad15mInt )  ) * Rad15mInt ) < DelTime ) then

      FlagCalcTrans = .true.

!!$      call MessageNotify( 'M', module_name, 'Transmission is calculated.' )

      !
      ! Calculation of "absorption" dust optical depth
      ! This formulation is obtained from Forget et al. [1999].
      !
      do k = 0, kmax
        do j = 1, jmax
          do i = 0, imax-1
            xyr_DOD(i,j,k) = ( 1.0d0 - SSA ) * xyr_DOD067(i,j,k) * QeRat
          end do
        end do
      end do

!!$      call increase_vreso_boundary_arr3d( im, jm, km, nvr, gdod, gdod_f, &
!!$        & "log", ijs, ije )


      !
      ! check pressure
      !
      MinPress = 1.0d100
      MaxPress = 0.0d0
      do j = 1, jmax
        do i = 0, imax-1
          MinPress = min( MinPress, xyz_Press(i,j,kmax) )
          MaxPress = max( MaxPress, xyz_Press(i,j,1   ) )
        end do
      end do
      if( ckdp(1)%lnp(1) > log(MinPress) ) then
        write( 6, * ) 'MARS: pressure is too small.'
        write( 6, * ) MinPress, exp(ckdp(1)%lnp(1))
        stop
      end if
      if( ckdp(1)%lnp(ckdp(1)%nlnp) .lt. log(MaxPress) ) then
        write( 6, * ) 'MARS: pressure is too large.'
        write( 6, * ) MaxPress, exp(ckdp(1)%lnp(ckdp(1)%nlnp))
        stop
      end if


      xyz_MMMass = 43.5d0 * AMU
      do n = 1, nras + nrps
        do k = 1, kmax
          do j = 1, jmax
            do i = 0, imax-1
              xyza_VMR(i,j,k,n) = VMRCO2
            end do
          end do
        end do
      end do


!!$      do n = 1, nras + nrps
!!$        call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$          & gvmrh(:,:,:,n), gvmr_f(:,:,:,n), "log", &
!!$          & ijs, ije )
!!$      end do
!!$      call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$        & mmmassh, mmmass_f, "linear", &
!!$        & ijs, ije )

!!$      call calc_lnp( im, jm, km+1    , gph   , glnph   , ijs, ije )
      call calc_lnp( xyz_Press, xyz_lnPress )
      xy_lnPs(:,:) = log( xyr_Press(:,:,0) )


      !
      ! initialization
      !
      do k = 0, kmax
        do j = 1, jmax
          do i = 0, imax-1
            trans_i2i_toa(i,j,k) = 0.0d0         ! f_{1/2}    T_{k+1/2,1/2}
            trans_i2i_boa(i,j,k) = 0.0d0         ! f_{km+1/2} T_{k+1/2,km+1/2}
            trans_i2i_s  (i,j,k) = 0.0d0         ! f_{s}      T_{k+1/2,km+1/2}
          end do
        end do
      end do
      do k2 = 1, kmax
        do k = 0, kmax
          do j = 1, jmax
            do i = 0, imax-1
              trans_i2m_uli(i,j,k,k2) = 0.0d0
              trans_i2m_lli(i,j,k,k2) = 0.0d0
            end do
          end do
        end do
      end do


      !
      ! loop for wavenumber
      !

      iband_reserve = 0

      do m = 1, nwnl

        call m2ckdpindices( m, ig, iband )


        if( iband .ne. iband_reserve ) then
          call findindices3D( xyz_Temp, xyz_lnPress, iband, xyz_jj, xyz_kk )
          call findindices2D( xy_SurfTemp, xy_lnPs, iband, xy_jj, xy_kk )

          iband_reserve = iband
        end if


        ! IMPORTANT!
        ! This loop for n is confusing.
        ! We have to reconsider about it.
        ! Maybe, the component of ckdp structure has to be reconsidered. 
        ! Now, it cannot include multiple radiatively active species. 
        ! (yot, 2010/09/12)
        !
        do n = 1, nras
          call getlnac_givenindices( xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, xyza_AC(:,:,:,n) )
        end do
        do n = 1, nras
          xyza_AC(:,:,:,n) = exp( xyza_AC(:,:,:,n) )
        end do

!!$        do n = 1, nras
!!$          call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$            & ach(:,:,:,n), ac_f(:,:,:,n), "log", &
!!$            & ijs, ije )
!!$        end do


        call calc_trans_mp_arr3d( nras, nrps, xyr_Press, xyza_VMR, xyz_MMMass, xyza_AC, xyr_DOD, xyra_Trans )


        call getpfr_givenindices3D( xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, xyz_PFRat )
        xyr_PFRat(:,:,0) = xyz_PFRat(:,:,1)
        do k = 1, kmax-1
          xyr_PFRat(:,:,k) = ( xyz_PFRat(:,:,k) + xyz_PFRat(:,:,k+1) ) * 0.5_DP
        end do
        xyr_PFRat(:,:,kmax) = xyz_PFRat(:,:,kmax)
        call getpfr_givenindices2D( xy_SurfTemp, xy_lnPs, xy_jj, xy_kk, ig, iband, xy_SurfPFRat )

        do k = 0, kmax
          do j = 1, jmax
            do i = 0, imax-1
              trans_i2i_toa(i,j,k) = trans_i2i_toa(i,j,k) + xyra_Trans(i,j,k,kmax) * xyr_PFRat(i,j,kmax) * ckdp(iband)%weight(ig)
              trans_i2i_boa(i,j,k) = trans_i2i_boa(i,j,k) + xyra_Trans(i,j,k,0) * xyr_PFRat(i,j,0) * ckdp(iband)%weight(ig)
              trans_i2i_s  (i,j,k) = trans_i2i_s  (i,j,k) + xyra_Trans(i,j,k,0) * xy_SurfPFRat(i,j) * ckdp(iband)%weight(ig)
            end do
          end do
        end do

        do k2 = 1, kmax
          do k = 0, kmax
            do j = 1, jmax
              do i = 0, imax-1
                trans_i2m_uli(i,j,k,k2) = trans_i2m_uli(i,j,k,k2) + ( xyra_Trans(i,j,k,k2-1) + xyra_Trans(i,j,k,k2) ) * 0.5d0 * xyr_PFRat(i,j,k2  ) * ckdp(iband)%weight(ig)
                trans_i2m_lli(i,j,k,k2) = trans_i2m_lli(i,j,k,k2) + ( xyra_Trans(i,j,k,k2-1) + xyra_Trans(i,j,k,k2) ) * 0.5d0 * xyr_PFRat(i,j,k2-1) * ckdp(iband)%weight(ig)
              end do
            end do
          end do
        end do


      end do


!!$      call rad15m_rv_put_newscheme2006( time, gt, gph, gp, gts, gdod, ijs, ije )

    else

      if ( trans_i2i_toa(0,1,1) > 1.0d99 ) then
        write( 6, * ) 'transmission function would not be calculated.'
        stop
      end if

    end if


    ! Is this OK?
    iband = 1

    call getpf_arr3d_norat( xyr_Temp, xy_SurfTemp, iband, xyr_PF, xy_SurfPF )


    call calc_rteq_use_meantrans_arr3d( ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) ), xy_SurfEmis, trans_i2i_toa, trans_i2i_boa, trans_i2i_s, trans_i2m_lli, trans_i2m_uli, xyr_PF, xy_SurfPF, xyr_Rad15mFlux )


    do l = 0, 1
      do k = 0, kmax
        do j = 1, jmax
          do i = 0, imax-1
            xyra_DelRad15mFlux(i,j,k,l) = 0.0_DP
          end do
        end do
      end do
    end do


!!$    do k = kmax, 0, -1
!!$      write( 6, * ) gph(0,1,k), gr15mnetflh(0,1,k)
!!$    end do
!!$    stop


!!$      ij = ( ije - ijs + 1 ) / 2
!!$      k  = km + 1
!!$!      write( 6, * ) 'RAD15M : ', gr15mnetflh(ij,1,k), rad_gdr15mnetfldtsh0(ij,1,k-1) * ( gts(ij,1) - rad_gtsbase(ij,1,1) ), emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1)
!!$!      write( 61, * ) gr15mnetflh(ij,1,k-1), rad_gdr15mnetfldtsh0(ij,1,k-1) * ( gts(ij,1) - rad_gtsbase(ij,1,1) ), emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1), gt(ij,1,km), gt(ij,1,km-1), gt(ij,1,km-2), gt(ij,1,km-3), gt(ij,1,km-4), gt(ij,1,km-5)
!!$      write( 6, * ) 'RAD15M : ', gr15mnetflh(ij,1,k), &
!!$        & emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1)
!!$      write( 61, * ) gr15mnetflh(ij,1,k-1), &
!!$        & emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1), &
!!$        & gt(ij,1,km), gt(ij,1,km-1), gt(ij,1,km-2), gt(ij,1,km-3), &
!!$        & gt(ij,1,km-4), gt(ij,1,km-5)
!!$      call flush( 61 )

    !
    ! output variables
    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        goru(i,j) = uwflh_sum(i,j,kmax)
!!$        gord(i,j) = 0.0d0
!!$        gsru(i,j) = uwflh_sum(i,j,0)
!!$        gsrd(i,j) = dwflh_sum(i,j,0)
!!$        gor (i,j) = goru(i,j) - gord(i,j)
!!$        gsr (i,j) = gsru(i,j) - gsrd(i,j)
!!$      end do
!!$    end do



  end subroutine rad15m_lowatm_newscheme2006