subroutine radiation_sounding_init
    !
    ! サウンディングファイルから加熱率を読み込む
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !内部変数
    integer             :: AltCol = 0       !「高度」の列番号 (サウンディングファイル内)
    integer             :: SWaveCol = 0     !「短波放射による加熱率」の列番号 (サウンディングファイル内)
    integer             :: LWaveCol = 0     !「長波放射による加熱率」の列番号 (サウンディングファイル内)
    integer             :: unit             !設定ファイル用装置番号  
    integer             :: k
    integer             :: io
    character(30)       :: SoundingFile     ! サウンディングファイル
    
    integer, parameter  :: maxch=12
    character(len=200)  :: buf, eachcol(maxch)
    integer             :: i, num, MaxCol

    real(DP)            :: r_tmpQradSW(10000)
    real(DP)            :: r_tmpQradLW(10000)
    real(DP)            :: r_tmpAlt(10000)
    real(DP)            :: r_tmpQrad(10000)
    integer             :: NumRec = 0

    real(DP), allocatable :: xyr_PTempRadSndg(:,:,:)   

    logical             :: flag

   
    !設定ファイルから読み込む出力ファイル情報
    !
    NAMELIST /radiation_sounding_nml/ SoundingFile, AltCol, SWaveCol, LWaveCol, FlagDExnerDtRad
              
    !設定ファイルから出力ファイルに記載する情報を読み込む
    !
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=radiation_sounding_nml)
    close(unit)
      
    ! 初期化
    !
    io = 0
    NumRec = 0
    MaxCol = max( AltCol, max( SWaveCol, LWaveCol ) )
!    write(*,*) "MaxCol", MaxCol
    r_tmpAlt  = 0.0d0
     r_tmpQrad = 0.0d0
     r_tmpQradSW = 0.0d0
     r_tmpQradLW = 0.0d0

    ! ファイルのオープン
    !
    open (17, file=SoundingFile, status='old')
    
    ! ファイル呼び出し
    !
    do while ( io == 0 ) 
      ! 1 行分読み出し
      !
      read (17, '(a)', IOSTAT=io) buf
      
      ! 行をカンマ区切りで分割
      !
      call devidecsv( buf, eachcol, maxch, num )
      
      ! 確認
      !
!      write(*,*) num
!      do i=1, num
!        write(*,*) i, eachcol(i)(1:len_trim(eachcol(i)))
!      end do
      
      ! num の値が小さいものはヘッダとみなす. 
      !
      if (num >= MaxCol) then 
        ! 行数の計算
        !
        NumRec = NumRec + 1        
        
        ! 値の代入
        !
        if (AltCol > 0)   read( eachcol(AltCol)(1:len_trim(eachcol(AltCol))), *)   r_tmpAlt(NumRec) 
        if (SWaveCol > 0) read( eachcol(SWaveCol)(1:len_trim(eachcol(SWaveCol))), *) r_tmpQradSW(NumRec) 
        if (LWaveCol > 0) read( eachcol(LWaveCol)(1:len_trim(eachcol(LWaveCol))), *) r_tmpQradLW(NumRec) 

        r_tmpQrad(NumRec) = r_tmpQradSW(NumRec) + r_tmpQradLW(NumRec)

      end if

!      write(*,*) r_tmpAlt(NumRec), r_tmpQradSW(NumRec), r_tmpQradLW(NumRec), r_tmpQrad(NumRec) 

    end do

    ! ファイルのクローズ
    !
    close (17)

    !初期化
    !
    allocate( xyz_PTempRadSndg(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyr_PTempRadSndg(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_ExnerRadSndg(imin:imax, jmin:jmax, kmin:kmax) )

    xyz_PTempRadSndg = 0.0d0
    xyr_PTempRadSndg = 0.0d0
    xyz_ExnerRadSndg = 0.0d0

    flag = .false.
    
    do k = kmin, kmax
      do i = 1, NumRec
        if ( r_Z(k) == r_tmpAlt(i) ) then           
          xyr_PTempRadSndg(:,:,k) = r_tmpQrad(i)         
        end if
      end do
    end do

    do k = kmin, kmax
      do i = 1, NumRec
        if ( z_Z(k) == r_tmpAlt(i) ) then 
          flag = .true. 
!          write(*,*) r_tmpAlt(i), r_tmpQrad(i), '(', r_tmpQradSW(i), '+', r_tmpQradLW(i), ')'
          xyz_PTempRadSndg(:,:,k) = r_tmpQrad(i)
        end if
      end do
    end do

    if (.NOT. flag) then 
      do k = kmin+1, kmax
        xyz_PTempRadSndg(:,:,k) = ( xyr_PTempRadSndg(:,:,k-1) + xyr_PTempRadSndg(:,:,k) ) * 5.0d-1
      end do
    end if


    ! 単位換算
    !
    xyz_PTempRadSndg = xyz_PTempRadSndg / DayTime / xyz_ExnerBZ


    if ( FlagDExnerDtRad ) then    
      xyz_ExnerRadSndg = xyz_DExnerDt_xyz(xyz_PTempRadSndg)
    else
      xyz_ExnerRadSndg = 0.0d0
    end if


    ! ヒストリデータ定義
    ! 
    call HistoryAutoAddVariable( varname='PTempRad', dims=(/'x','y','z','t'/), longname='Radiation term of potential temperature', units='K.s-1', xtype='float')

    call HistoryAutoAddVariable( varname='ExnerRad', dims=(/'x','y','z','t'/), longname='Radiation term of Exner function', units='K.s-1', xtype='float')

  end subroutine radiation_sounding_init