!= Module Radiation_3d
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: radiation_3d.f90,v 1.2 2007-10-01 08:57:59 odakker Exp $
! Tag Name::  $Name: arare4-20080627 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!モデルの放射過程を計算するためのパッケージ型モジュール
!具体的には以下の項を計算するための関数を格納する.  
!  * 一様冷却
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!
!== Future Plans
!
!

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

  !モジュール読み込み
  use dc_message, only: MessageNotify
  use dc_iounit,  only: FileOpen

  use GridSet_3d,  only: DimXMin,     & !x 方向の配列の下限
    &                 DimXMax,     & !x 方向の配列の上限
    &                 DimYMin,     & !y 方向の配列の下限
    &                 DimYMax,     & !y 方向の配列の上限
    &                 DimZMin,     & !z 方向の配列の下限
    &                 DimZMax,     & !z 方向の配列の上限
    &                 z_Z            !Z 座標軸(スカラー格子点)
  use TimeSet,  only: DayTime        ! 1 日の長さ [s]
  use basicset_3d, only: xyz_ExnerBasicZ !エクスナー関数の基本場
  use StorePotTemp_3d,  only: StorePotTempRad, StorePotTempDamp
  
  !暗黙の型宣言禁止
  implicit none

  !private 属性
  private

  !関数を public にする. 
  public Radiation_init
  public xyz_RadHeatConst
  public xyz_NewtonCool

  !変数定義
  real(8)              :: RadHeatRate  = 0.0d0 !一様放射加熱率 [K/day]
  real(8), allocatable :: xyz_RadHeight(:,:,:)    !放射加熱が存在する領域
  
  !save 属性
  save RadHeatRate, xyz_RadHeight

contains

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

    !暗黙の型宣言禁止
    implicit none

    !入力変数
    character(*), intent(in) :: cfgfile
    real(8)                  :: RadHeightUp   !放射強制を与える鉛直領域の上限
    real(8)                  :: RadHeightDown !放射強制を与える鉛直領域の下限
    integer                  :: k             !ループ変数
    integer                  :: unit            !装置番号

    ! NAMELIST から情報を取得
    NAMELIST /radiation/ RadHeatRate, RadHeightUp, RadHeightDown

    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=radiation)
    close(unit)


    allocate( xyz_RadHeight(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) )

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

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

    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 subroutine Radiation_init


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

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

    xyz_RadHeatConst = &
      & xyz_RadHeight * RadHeatRate  &
      & / ( ( xyz_ExnerBasicZ + xyz_Exner ) * DayTime )

    call StorePotTempRad( xyz_RadHeatConst )
    
  end function xyz_RadHeatConst
  
  
  
!!!------------------------------------------------------------------------!!!
  function xyz_NewtonCool(xyz_PotTemp)
    !
    ! 温位の放射強制項. 
    ! 地表面から RadHeight で指定された高度までの間で一様放射冷却を与える. 
    !
    ! = バグ
    ! 
    ! 上記の仕様通り設計されていない. 要検討. 
    
    !暗黙の型宣言を禁止
    implicit none

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

    !放射緩和時間は 5 日
    EFTime = 5.0d0 * DayTime  
    
    do k = DimZMin, DimZMax
      xyz_PotTempMean(:,:,k) = &
        & sum( xyz_PotTemp(:,:,k) ) / &
        & (real(DimXMax - DimXMin + 1) * real(DimYMax - DimYMin + 1) )
    end do
    
    xyz_NewtonCool = - (xyz_PotTemp - xyz_PotTempMean) / EFTime 
    
    call StorePotTempDamp( xyz_NewtonCool )

  end function xyz_NewtonCool
  
end module Radiation_3d
