!= 邁｡蜊俶叛蟆�Δ繧ｸ繝･繝ｼ繝ｫ: 縺ｨ縺ゅｋ鬮伜ｺｦ鬆伜沺繧剃ｸ讒伜�蜊ｴ繝ｻ蜉�辭ｱ縺吶ｋ
!
! Authors::   YAMASHITA Tatsuya, SUGIYAMA Ko-ichiro
! Version::   $Id: radiation_heatbalance.f90,v 1.13 2014/07/08 01:04:18 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]

module Radiation_HeatBalance
  !
  ! 邁｡蜊俶叛蟆�Δ繧ｸ繝･繝ｼ繝ｫ: 縺ｨ縺ゅｋ鬮伜ｺｦ鬆伜沺繧剃ｸ讒伜�蜊ｴ繝ｻ蜉�辭ｱ縺吶ｋ
  !

  !繝｢繧ｸ繝･繝ｼ繝ｫ隱ｭ縺ｿ霎ｼ縺ｿ
  use dc_types, only: DP, STRING

  !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
  implicit none

  !private 螻樊ｧ
  private

  !螟画焚螳夂ｾｩ
  real(DP),save   :: RadCoolRate = 0.0d0    !荳讒俶叛蟆�刈辭ｱ邇� [K/day]
  integer, save   :: IdxHeatUp   = 0        !蜉�辭ｱ蝓滉ｸ企剞縺ｮ驩帷峩蠎ｧ讓吶↓蟇ｾ蠢懊☆繧区紛謨ｰ蛟､
  integer, save   :: IdxHeatDown = 0        !蜀ｷ蜊ｴ蝓滉ｸ企剞縺ｮ驩帷峩蠎ｧ讓吶↓蟇ｾ蠢懊☆繧区紛謨ｰ蛟､
  integer, save   :: IdxCoolUp   = 0        !蜉�辭ｱ蝓滉ｸ企剞縺ｮ驩帷峩蠎ｧ讓吶↓蟇ｾ蠢懊☆繧区紛謨ｰ蛟､
  integer, save   :: IdxCoolDown = 0        !蜀ｷ蜊ｴ蝓滉ｸ企剞縺ｮ驩帷峩蠎ｧ讓吶↓蟇ｾ蠢懊☆繧区紛謨ｰ蛟､
  real(DP),save, allocatable :: xyz_RadHeightHeat(:,:,:)  !謾ｾ蟆�刈辭ｱ縺悟ｭ伜惠縺吶ｋ鬆伜沺
  real(DP),save, allocatable :: xyz_RadHeightCool(:,:,:)  !謾ｾ蟆�刈辭ｱ縺悟ｭ伜惠縺吶ｋ鬆伜沺

  real(DP), save  :: FactorDExnerDtRad = 1.0d0

  character(STRING), parameter:: module_name = 'Radiation_Heatbalance'

  !髢｢謨ｰ繧� public 縺ｫ縺吶ｋ.
  public Radiation_heatbalance_init
  public Radiation_heatbalance_forcing

contains

!!!------------------------------------------------------------------------!!!
  subroutine Radiation_heatbalance_init
    !
    !NAMELIST 縺九ｉ蜀ｷ蜊ｴ邇�, 蜀ｷ蜊ｴ縺吶ｋ鬮伜ｺｦ鬆伜沺, 蜉�辭ｱ縺吶ｋ鬮伜ｺｦ鬆伜沺, 繧定ｨｭ螳壹☆繧�.
    !蜉�辭ｱ邇��豈弱せ繝�ャ繝苓ｨ育ｮ励☆繧九�縺ｧ, 蛻晄悄蛹悶Ν繝ｼ繝√Φ縺ｮ荳ｭ縺ｧ縺ｯ豎ｺ繧√↑縺�.
    !

    use dc_types,          only : DP
    use dc_iounit,         only : FileOpen
    use dc_message,        only : MessageNotify
    use gtool_historyauto, only : HistoryAutoAddVariable
    use namelist_util,     only : namelist_filename
    use gridset,           only : imin,     & !x 譁ｹ蜷代�驟榊�縺ｮ荳矩剞
      &                           imax,     & !x 譁ｹ蜷代�驟榊�縺ｮ荳企剞
      &                           jmin,     & !y 譁ｹ蜷代�驟榊�縺ｮ荳矩剞
      &                           jmax,     & !y 譁ｹ蜷代�驟榊�縺ｮ荳企剞
      &                           kmin,     & !z 譁ｹ蜷代�驟榊�縺ｮ荳矩剞
      &                           kmax        !z 譁ｹ蜷代�驟榊�縺ｮ荳企剞
    use axesset,           only : z_Z         !Z 蠎ｧ讓呵ｻｸ(繧ｹ繧ｫ繝ｩ繝ｼ譬ｼ蟄千せ)

    !證鈴ｻ吶�蝙句ｮ｣險遖∵ｭ｢
    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, &
      & FactorDExnerDtRad

    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
    !
    call MessageNotify( "M", &
      & module_name, "RadCoolRate = %f", d=(/RadCoolRate/))
    call MessageNotify( "M", &
      & module_name, "HeightHeatUp = %f", d=(/HeightHeatUp/))
    call MessageNotify( "M", &
      & module_name, "HeightHeatDown= %f", d=(/HeightHeatDown/))
    call MessageNotify( "M", &
      & module_name, "HeightCoolUp = %f", d=(/HeightCoolUp/))
    call MessageNotify( "M", &
      & module_name, "HeightCoolDown= %f", d=(/HeightCoolDown/))
    call MessageNotify( "M", &
      & module_name, "FactorDExnerDtRad= %f", d=(/ FactorDExnerDtRad /))

    ! 繝偵せ繝医Μ繝��繧ｿ螳夂ｾｩ
    !
    call HistoryAutoAddVariable(  &
      & varname='DPTempDtRad', &
      & 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='s-1"',    &
      & xtype='float')

  end subroutine Radiation_heatbalance_init

!!!------------------------------------------------------------------------!!!
  subroutine radiation_heatbalance_forcing(xyz_Exner, xyz_PTemp, xyz_DPTempDt, xyz_DExnerDt)
    !
    ! 貂ｩ菴阪�謾ｾ蟆�ｼｷ蛻ｶ鬆�.
    ! 蝨ｰ陦ｨ髱｢縺九ｉ RadHeight 縺ｧ謖�ｮ壹＆繧後◆鬮伜ｺｦ縺ｾ縺ｧ縺ｮ髢薙〒遨ｺ髢鍋噪縺ｫ荳讒倥↑蜉�辭ｱ,
    ! RadHeight 縺九ｉ RadHeight2 縺ｾ縺ｧ縺ｮ髢薙〒遨ｺ髢鍋噪縺ｪ荳讒倥↑蜀ｷ蜊ｴ繧剃ｸ弱∴繧�.
    ! 謾ｾ蟆�ｼｷ蛻ｶ蜈ｨ菴薙→縺励※蜉�辭ｱ縺ｨ蜀ｷ蜊ｴ縺後ヰ繝ｩ繝ｳ繧ｹ縺吶ｋ繧医≧縺ｫ譎ゅ�綾縲�ｸ讒伜刈辭ｱ邇�
    ! 縺翫ｈ縺ｳ荳讒伜�蜊ｴ邇�ｒ螟牙喧縺輔○繧�.
    ! 蜀ｷ蜊ｴ縺ｮ謖ｯ蟷�ｒ荳弱∴, 蜉�辭ｱ縺ｮ謖ｯ蟷�ｒ隱ｿ遽縺吶ｋ.

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蜻ｼ縺ｳ蜃ｺ縺�
    !
    use gtool_historyauto, only : HistoryAutoPut
    use timeset,           only : TimeN
    use gridset,           only : imin,     & !x 譁ｹ蜷代�驟榊�縺ｮ荳矩剞
      &                           imax,     & !x 譁ｹ蜷代�驟榊�縺ｮ荳企剞
      &                           jmin,     & !y 譁ｹ蜷代�驟榊�縺ｮ荳矩剞
      &                           jmax,     & !y 譁ｹ蜷代�驟榊�縺ｮ荳企剞
      &                           kmin,     & !z 譁ｹ蜷代�驟榊�縺ｮ荳矩剞
      &                           kmax,     & !z 譁ｹ蜷代�驟榊�縺ｮ荳企剞
      &                           nx, ny, nz
    use DExnerDt,          only : xyz_DExnerDt_xyz
    use axesset,           only : z_dz

    use constants,         only : DayTime,  & ! 1 譌･縺ｮ髟ｷ縺� [s]
      &                           PressBasis, & !貂ｩ菴阪�蝓ｺ貅門悸蜉�
      &                           CpDry,      & !螳壼悸豈皮�
      &                           CvDry,      & !螳夂ｩ肴ｯ皮�
      &                           GasRDry       !豌嶺ｽ灘ｮ壽焚
    use basicset,          only : xyz_ExnerBZ,  &!繧ｨ繧ｯ繧ｹ繝翫�髢｢謨ｰ縺ｮ蝓ｺ譛ｬ蝣ｴ
      &                           xyz_PTempBZ    !貂ｩ菴阪�蝓ｺ譛ｬ蝣ｴ

    !證鈴ｻ吶�蝙句ｮ｣險繧堤ｦ∵ｭ｢
    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 = PressBasis                            &
      & * (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_DExnerDt_xyz( xyz_Rad ) * FactorDExnerDtRad

    ! output
    !
    call HistoryAutoPut(TimeN, 'DPTempDtRad', 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

end module Radiation_HeatBalance