!= Module Radiation
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: radiation.f90,v 1.14 2011-03-30 13:53:08 sugiyama Exp $
! Tag Name::  $Name: arare4-20120511 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!モデルの放射過程を計算するためのパッケージ型モジュール
!具体的には以下の項を計算するための関数を格納する.  
!  * 一様冷却
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!
!== Future Plans
!
!

module Radiation
  !
  !モデルの放射過程を計算するためのパッケージ型モジュール
  !具体的には以下の項を計算するための関数を格納する.  
  !  * 一様冷却
  !

  !モジュール読み込み
  use dc_message, only: MessageNotify
  use debugset,   only: cpurank
  use GridSet,  only: DimXMin,     & !x 方向の配列の下限
    &                 DimXMax,     & !x 方向の配列の上限
    &                 DimZMin,     & !z 方向の配列の下限
    &                 DimZMax,     & !z 方向の配列の上限
    &                 RegZMin,     & !x 方向の配列の下限
    &                 RegXMax,     & !x 方向の配列の上限
    &                 RegXMin,     & !z 方向の配列の下限
    &                 RegZMax,     & !z 方向の配列の上限
    &                 DelZ,        & !格子間隔
    &                 s_Z            !Z 座標軸(スカラー格子点)
  use TimeSet,  only: DayTime        ! 1 日の長さ [s]
!  use basicset, only: xz_ExnerBasicZ !エクスナー関数の基本場
  use basicset, only: xz_ExnerBasicZ, & !エクスナー関数の基本場
    &                 xz_PressBasicZ, & !圧力の基本場
    &                 xz_PotTempBasicZ, & !温位の基本場
    &                 PressSfc,       & !地表面圧力
    &                 CpDry,          & !定圧比熱
    &                 CvDry,          & !定積比熱
    &                 GasRDry,        & !気体定数
    &                 Grav              !重力加速度
  use StorePotTemp,  only: StorePotTempRad, StorePotTempDamp
  
  !暗黙の型宣言禁止
  implicit none

  !private 属性
  private

  !関数を public にする. 
  public Radiation_init
  public xz_RadHeatConst
  public xz_RadHeatVary
  public xz_NewtonCool
  public xz_radheatbalance

  !変数定義
  real(8), save              :: RadHeatRate = 0.0d0  !一様放射加熱率 [K/day]
  real(8), save              :: RadCoolRate = 0.0d0  !一様放射加熱率 [K/day]
  real(8), save, allocatable :: xz_RadHeight(:,:)    !放射加熱が存在する領域
  real(8), save, allocatable :: xz_RadHeight2(:,:)   !放射加熱が存在する領域
  real(8), save              :: RadHeight1 = 0.0d0
  real(8), save              :: RadHeight2 = 0.0d0
  integer, save              :: IntHeightUp = 0.0d0  !加熱域上限の鉛直座標に対応する整数値
  integer, save              :: IntHeightUp2 = 0.0d0 !冷却域上限の鉛直座標に対応する整数値

contains

!!!------------------------------------------------------------------------!!!
  subroutine Radiation_init(cfgfile)
    !
    !NAMELIST から放射強制の設定を取得
    !

    !暗黙の型宣言禁止
    implicit none

    !入力変数
    character(*), intent(in) :: cfgfile
    real(8)                  :: RadHeightUp   = 0.0d0  !放射強制を与える鉛直領域の上限
    real(8)                  :: RadHeightDown = 0.0d0  !放射強制を与える鉛直領域の下限
    real(8)                  :: RadHeightUp2  = 0.0d0  !放射強制を与える鉛直領域の上限
    real(8)                  :: RadHeightDown2= 0.0d0  !放射強制を与える鉛直領域の下限
    integer                  :: k                      !ループ変数

    ! NAMELIST から情報を取得
    NAMELIST /radiation/ RadHeatRate, RadHeightUp, RadHeightDown, &
      & RadCoolRate, RadHeight1, RadHeight2, RadHeightUp2, RadHeightDown2

    open (10, FILE=cfgfile)
    read(10, NML=radiation)
    close(10)

    ! IntHeightUp の特定
    do k = DimZMin, DimZMax
      if ( (RadHeightUp - DelZ <= s_Z(k)) .and. (s_Z(k) < RadHeightUp) ) then
        IntHeightUp = k
      else if ( (RadHeightUp2 - DelZ <= s_Z(k)) .and. (s_Z(k) < RadHeightUp2) ) then
        IntHeightUp2 = k
      end if
    end do

    allocate( xz_RadHeight(DimXMin:DimXMax, DimZMin:DimZMax) )
    allocate( xz_RadHeight2(DimXMin:DimXMax, DimZMin:DimZMax) )

    ! 放射強制が存在する領域の設定
    ! RadHeightDown < s_Z < RadHeightUp の範囲では値が 1 となる
    ! 係数を用意する.
    !
    xz_RadHeight  = 1.0d0 
    xz_RadHeight2 = 1.0d0 

    do k = DimZMin, DimZMax
      if ( s_Z(k) <= RadHeightDown  ) then
        xz_RadHeight(:,k) = 0.0d0 
      elseif( s_Z(k) >= RadHeightUp ) then
        xz_RadHeight(:,k) = 0.0d0 
      end if
    end do

    do k = DimZMin, DimZMax
      if ( s_Z(k) <= RadHeightDown2  ) then
        xz_RadHeight2(:,k) = 0.0d0 
      elseif( s_Z(k) >= RadHeightUp2 ) then
        xz_RadHeight2(:,k) = 0.0d0 
      end if
    end do

    if (cpurank == 0) then 
      call MessageNotify( "M", &
        & "Radiation", "RadHeatRate = %f", d=(/RadHeatRate/))
      call MessageNotify( "M", &
        & "Radiation", "RadHeightUp = %f", d=(/RadHeightUP/))
      call MessageNotify( "M", &
        & "Radiation", "RadHeightDown= %f", d=(/RadHeightDown/))
    end if
  end subroutine Radiation_init


!!!------------------------------------------------------------------------!!!
  function xz_RadHeatConst(xz_Exner)
    !
    ! 温位の放射強制項. 
    ! 地表面から RadHeight で指定された高度までの間で一様放射冷却を与える. 
    !
    
    !暗黙の型宣言を禁止
    implicit none

    !変数定義
    real(8), intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                       !エクスナー関数の擾乱場
    real(8)               :: xz_RadHeatConst(DimXMin:DimXMax, DimZMin:DimZMax)
                                       !放射強制

    xz_RadHeatConst = &
      & xz_RadHeight * RadHeatRate  &
      & / ( ( xz_ExnerBasicZ + xz_Exner ) * DayTime )

    call StorePotTempRad( xz_RadHeatConst )
    
  end function xz_RadHeatConst
  
  
!!!------------------------------------------------------------------------!!!
  function xz_RadHeatVary(xz_Exner)
  !
  ! 温位の放射強制項.
  ! 地表面から RadHeight1 で指定された高度までの間は2[K/day](RadHeatRate)の
  ! 一様な放射冷却を与える.
  ! RadHeight1 から RadHeight2 までは (RadHeight2 - s_Z(k)) / (RadHeight2 - RadHeight1)
  ! の割合で放射冷却の効果を減少させる (RadHeight2でちょうどゼロにする)
  ! Nakajima and Matsuno(1988),中島(1944)を参考にした
  !

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

  !変数定義
  real(8), intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                     !エクスナー関数の擾乱場
  real(8)               :: xz_RadHeatVary(DimXMin:DimXMax, DimZMin:DimZMax)
                                     !放射強制
  integer               :: k

  do k = RegZMin+1, DimZMax
   if (s_Z(k).le.RadHeight1) then
    xz_RadHeatVary(:,k) = &
      & RadHeatRate  &
      & / ( ( xz_ExnerBasicZ(:,k) + xz_Exner(:,k) ) * DayTime )

   else if (s_Z(k).gt.RadHeight1.and.s_Z(k).le.RadHeight2) then
    xz_RadHeatVary(:,k) = &
      & RadHeatRate  * (RadHeight2 - s_Z(k)) &
      & / ( (RadHeight2 - RadHeight1) * ( xz_ExnerBasicZ(:,k) + xz_Exner(:,k) ) * DayTime )

   else if (s_Z(k).gt.RadHeight2) then
    xz_RadHeatVary(:,k) = 0.0d0

   end if
  end do

  call StorePotTempRad( xz_RadHeatVary )

end function xz_RadHeatVary

!!!------------------------------------------------------------------------!!!
  function xz_NewtonCool(xz_PotTemp)
    !
    ! 温位の放射強制項. 
    ! 地表面から RadHeight で指定された高度までの間で一様放射冷却を与える. 
    !
    
    !暗黙の型宣言を禁止
    implicit none

    !変数定義
    real(8), intent(in)   :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                       !温位の擾乱場
    real(8)               :: xz_NewtonCool(DimXMin:DimXMax, DimZMin:DimZMax)
                                       !放射強制
    real(8)               :: xz_PotTempMean(DimXMin:DimXMax, DimZMin:DimZMax)
                                       !温位擾乱の東西平均
    real(8)               :: EFTime    !放射緩和時間
    integer               :: k

    !放射緩和時間は 5 日
    EFTime = 5.0d0 * DayTime  
    
    do k = DimZMin, DimZMax
      xz_PotTempMean(:,k) = sum( xz_PotTemp(:,k) ) / real(DimXMax - DimXMin + 1)
    end do
    
    xz_NewtonCool = - xz_PotTempMean / EFTime 
    
    call StorePotTempDamp( xz_NewtonCool )

  end function xz_NewtonCool


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

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

    !変数定義
    real(8), intent(in)   :: xz_ExnerBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in)   :: xz_PotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(in)   :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)               :: xz_RadHeatBalance(DimXMin:DimXMax, DimZMin:DimZMax)
                                        !放射強制
    real(8)               :: xz_DensSum(DimXMin:DimXMax, DimZMin:DimZMax) 
                                        !密度(基本場成分と擾乱成分の和)
    real(8)               :: HeatSum    !ある時刻での単位時間当たりの加熱量
    real(8)               :: CoolSum    !ある時刻での単位時間当たりの冷却量
    integer               :: i
    integer               :: k

    ! 密度の算出
    xz_DensSum = PressSfc &
      & * (xz_ExnerBasicZ + xz_Exner)**( CvDry /GasRDry ) &
      & / (GasRDry * (xz_PotTempBasicZ + xz_PotTemp ) )

    HeatSum = 0.0d0
    CoolSum = 0.0d0

    !加熱するのは, 高度 0 - 1000 m => 格子点は 1 (RegZMin + 1) から IntHeightUp
    !冷却するのは, 高度 1000 - 15000 m => 格子点は IntHeightUp + 1 から IntHeightUp2
    !*** 抜粋 *** 
    ! IntHeightUp の特定
    !do k = DimZMin, DimZMax
    !  if ( (RadHeightUp - DelZ <= s_Z(k)) .and. (s_Z(k) < RadHeightUp) ) then
    !    IntHeightUp = k
    !  else if ( (RadHeightUp2 - DelZ <= s_Z(k)) .and. (s_Z(k) < RadHeightUp2) ) then
    !    IntHeightUp2 = k
    !  end if
    !end do
    !
    ! 下端                           1000 
    !   |--*--|--*--|--*--|--*--|--*--|--*--|--*--|
    !      1                       5  
    !
    !  1000                                      Max
    !   |--*--|--*--|--*--|--*--|--*--|--*--|--*--|
    !      6                                   I

    ! ある時刻での (加熱量/加熱率) を計算
    do k = RegZmin + 1, IntHeightUp 
       do i = RegXmin + 1, RegXMax
          HeatSum = HeatSum + DelZ * CpDry * xz_DensSum(i,k)
       end do
    end do

    ! ある時刻での (冷却量/冷却率) を計算
    do k = IntHeightUp + 1, IntHeightUp2 
       do i = RegXmin + 1, RegXMax
          CoolSum = CoolSum + DelZ * CpDry * xz_DensSum(i,k)
       end do
    end do

    ! 加熱率の算出
    ! RadCoolRate が負値であることに注意. 
    RadHeatRate = - RadCoolRate * CoolSum / HeatSum

    xz_RadHeatBalance = &
      & xz_RadHeight * RadHeatRate  &
      & / ( xz_ExnerBasicZ  * DayTime ) &
      & + xz_RadHeight2 * RadCoolRate  &
      & / ( xz_ExnerBasicZ  * DayTime )

    call StorePotTempRad( xz_RadHeatBalance )
    
  end function xz_RadHeatBalance
  
end module Radiation
