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

    !暗黙の型宣言禁止
    implicit none

    !入力変数
    real(8)                  :: HeightHeatUp   = 0.0d0  !放射強制を与える鉛直領域の上限
    real(8)                  :: HeightHeatDown = 0.0d0  !放射強制を与える鉛直領域の上限
    real(8)                  :: HeightCoolUp   = 0.0d0  !放射強制を与える鉛直領域の上限
    real(8)                  :: HeightCoolDown = 0.0d0  !放射強制を与える鉛直領域の上限
    integer                  :: k                    !ループ変数
    integer                  :: unit

    ! NAMELIST から情報を取得
    NAMELIST /radiation_heatbalance_nml/ RadCoolRate, HeightHeatUp, HeightHeatDown, HeightCoolUp, HeightCoolDown

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

    do k = kmin, kmax-1
      if ( z_Z(k) <= HeightHeatUp .AND. HeightHeatUp < z_Z(k+1) ) then 
        IdxHeatUp = k
      elseif ( z_Z(k) <= HeightHeatDown .AND. HeightHeatDown < z_Z(k+1) ) then 
        IdxHeatDown = k
      elseif ( z_Z(k) <= HeightCoolUp .AND. HeightCoolUp < z_Z(k+1) ) then 
        IdxCoolUp = k
      elseif ( z_Z(k) <= HeightCoolDown .AND. HeightCoolDown < z_Z(k+1) ) then 
        IdxCoolDown = k
      end if
    end do

    allocate( xyz_RadHeightHeat(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_RadHeightCool(imin:imax, jmin:jmax, kmin:kmax) )

    xyz_RadHeightHeat = 1.0d0
    xyz_RadHeightHeat(:,:,IdxHeatDown:IdxHeatUp) = 1.0d0
    xyz_RadHeightCool = 0.0d0
    xyz_RadHeightCool(:,:,IdxCoolDown:IdxCoolUp) = 1.0d0

    ! Output
    !
    if (myrank == 0) then 
      call MessageNotify( "M", "Radiation", "RadCoolRate = %f", d=(/RadCoolRate/))
      call MessageNotify( "M", "Radiation", "HeightHeatUp = %f", d=(/HeightHeatUp/))
      call MessageNotify( "M", "Radiation", "HeightHeatDown= %f", d=(/HeightHeatDown/))
      call MessageNotify( "M", "Radiation", "HeightCoolUp = %f", d=(/HeightCoolUp/))
      call MessageNotify( "M", "Radiation", "HeightCoolDown= %f", d=(/HeightCoolDown/))
    end if

    call HistoryAutoAddVariable( varname='PTempRad', dims=(/'x','y','z','t'/), longname='Radiation term of potential temperature', units='K.s-1"', xtype='double')

    call HistoryAutoAddVariable( varname='ExnerRad', dims=(/'x','y','z','t'/), longname='Radiation term of exner function', units='s-1"', xtype='double')

  end subroutine Radiation_heatbalance_init

  subroutine radiation_heatbalance_forcing(xyz_Exner, xyz_PTemp, xyz_DPTempDt, xyz_DExnerDt)
    !
    ! 温位の放射強制項. 
    ! 地表面から RadHeight で指定された高度までの間で空間的に一様な加熱, 
    ! RadHeight から RadHeight2 までの間で空間的な一様な冷却を与える. 
    ! 放射強制全体として加熱と冷却がバランスするように時々刻々一様加熱率
    ! および一様冷却率を変化させる. 
    ! 冷却の振幅を与え, 加熱の振幅を調節する. 

    !暗黙の型宣言を禁止
    implicit none

    !変数定義
    real(DP), intent(in)  :: xyz_Exner(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(in)  :: xyz_PTemp(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(inout) :: xyz_DPTempDt(imin:imax, jmin:jmax, kmin:kmax)
    real(DP), intent(inout) :: xyz_DExnerDt(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_DPTempDt0(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_DExnerDt0(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_Rad(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_RadPI(imin:imax, jmin:jmax, kmin:kmax)
    real(DP)              :: xyz_DensSum(imin:imax, jmin:jmax, kmin:kmax) 
                                        !密度(基本場成分と擾乱成分の和)
    real(DP)              :: HeatSum    !ある時刻での単位時間当たりの加熱量
    real(DP)              :: CoolSum    !ある時刻での単位時間当たりの冷却量
    real(DP)              :: RadHeatRate
    integer               :: i, j, k

    ! 初期化
    !
    HeatSum = 0.0d0
    CoolSum = 0.0d0
    xyz_DPTempDt0 = xyz_DPTempDt
    xyz_DExnerDt0 = xyz_DExnerDt

    ! 密度の算出
    xyz_DensSum = PressSfc * (xyz_ExnerBZ + xyz_Exner)**( CvDry /GasRDry ) / (GasRDry * (xyz_PTempBZ + xyz_PTemp ) )   

    ! ある時刻での (加熱量/加熱率) を計算
    do k = IdxHeatDown, IdxHeatUp
      do j = 1, ny
        do i = 1, nx
          HeatSum = HeatSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k)
        end do
      end do
    end do
    
    ! ある時刻での (冷却量/冷却率) を計算
    do k = IdxCoolDown, IdxCoolUp
      do j = 1, ny
        do i = 1, nx
          CoolSum = CoolSum + z_dz(k) * CpDry * xyz_DensSum(i,j,k)
        end do
      end do
    end do

    ! 加熱率の算出
    ! RadCoolRate が負値であることに注意. 
    RadHeatRate = - RadCoolRate * CoolSum / HeatSum
    
    xyz_Rad = xyz_RadHeightHeat * RadHeatRate / ( xyz_ExnerBZ  * DayTime ) + xyz_RadHeightCool * RadCoolRate / ( xyz_ExnerBZ  * DayTime )

    xyz_DPTempDt = xyz_DPTempDt0 + xyz_Rad

    ! 圧力変化
    !
    xyz_RadPI = (xyz_VelSoundBZ ** 2.0d0) * xyz_Rad / (CpDry * xyz_VPTempBZ ** 2.0d0)
    
    xyz_DExnerDt = xyz_DExnerDt0 + xyz_RadPI

    ! output
    !
    call HistoryAutoPut(TimeN, 'PTempRad', xyz_Rad(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerRad', xyz_RadPI(1:nx,1:ny,1:nz))
    
  end subroutine radiation_heatbalance_forcing