subroutine Radiation_Simple_init
    !
    !NAMELIST から放射強制の設定を取得
    !

    !暗黙の型宣言禁止
    implicit none

    !入力変数
    real(DP)                 :: HeightUp   = 0.0d0  !放射強制を与える鉛直領域の上限
    real(DP)                 :: HeightDown = 0.0d0  !放射強制を与える鉛直領域の下限
    real(DP)                 :: RadHeatRate = 0.0d0   !一様放射加熱率 [K/day]
    integer                  :: k                      !ループ変数
    integer                  :: unit

    ! NAMELIST から情報を取得
    NAMELIST /radiation_simple_nml/ RadHeatRate, HeightUp, HeightDown

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=radiation_simple_nml)
    close(unit)

    allocate( xyz_PTempRadVary(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_PTempRadConst(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_ExnerRadVary(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_ExnerRadConst(imin:imax, jmin:jmax, kmin:kmax) )

    
    ! 温位の放射強制項. 
    ! 地表面から RadHeight で指定された高度までの間で一様放射冷却を与える. 
    do k = kmin, kmax
      if ( z_Z(k) <= HeightDown  ) then
        xyz_PTempRadConst(:,:,k) = 0.0d0 
      elseif( z_Z(k) >= HeightUp ) then
        xyz_PTempRadConst(:,:,k) = 0.0d0 
      else
        xyz_PTempRadConst(:,:,k) = RadHeatRate / DayTime / xyz_ExnerBZ(:,:,k)
      end if
    end do
    
    xyz_ExnerRadConst = xyz_DExnerDt_xyz(xyz_PTempRadConst)
    
    ! 温位の放射強制項.
    ! 地表面から Height1 で指定された高度までの間は2[K/day](RadHeatRate)の
    ! 一様な放射冷却を与える.
    ! Height1 から Height2 までは (Height2 - z_Z(k)) / (Height2 - Height1)
    ! の割合で放射冷却の効果を減少させる (Height2でちょうどゼロにする)
    ! Nakajima and Matsuno(1988),中島(1944)を参考にした
    !
    do k = kmin, kmax
      if ( z_Z(k) <= HeightDown ) then
        xyz_PTempRadVary(:,:,k) = RadHeatRate / DayTime / xyz_ExnerBZ(:,:,k)
        
      elseif ( z_Z(k) > HeightDown .AND. z_Z(k) <= HeightUP) then
        xyz_PTempRadVary(:,:,k) = RadHeatRate / DayTime / xyz_ExnerBZ(:,:,k) * (HeightUP - z_Z(k)) / (HeightUP - HeightDown) 
        
      else if (z_Z(k) > HeightUP) then
        xyz_PTempRadVary(:,:,k) = 0.0d0
      end if
    end do
 
    xyz_ExnerRadVary = xyz_DExnerDt_xyz(xyz_PTempRadVary) 


    ! Output
    !
    if (myrank == 0) then 
      call MessageNotify( "M", "Radiation", "RadHeatRate = %f", d=(/RadHeatRate/))
      call MessageNotify( "M", "Radiation", "HeightUp = %f", d=(/HeightUP/))
      call MessageNotify( "M", "Radiation", "HeightDown= %f", d=(/HeightDown/))
    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_Simple_init