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

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

  !モジュール読み込み
  use dc_types, only: DP, STRING
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify
  use gtool_historyauto, only: HistoryAutoAddVariable

  use mpi_wrapper,only: myrank
  use gridSet,  only: imin,     & !x 方向の配列の下限
    &                 imax,     & !x 方向の配列の上限
    &                 jmin,     & !y 方向の配列の下限
    &                 jmax,     & !y 方向の配列の上限
    &                 kmin,     & !z 方向の配列の下限
    &                 kmax        !z 方向の配列の上限
  use axesset,  only: z_Z         !Z 座標軸(スカラー格子点)
  use constants,only: DayTime     ! 1 日の長さ [s]
  use basicset, only: xyz_ExnerBZ !エクスナー関数の基本場
  use namelist_util, only: namelist_filename
  
  !暗黙の型宣言禁止
  implicit none

  !private 属性
  private

  !変数定義
  real(DP), save, allocatable, public :: xyz_RadHeatVary(:,:,:)   !放射加熱が存在する領域
  real(DP), save, allocatable, public :: xyz_RadHeatConst(:,:,:)  !放射加熱が存在する領域

  public Radiation_Simple_init

contains

!!!------------------------------------------------------------------------!!!
  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_RadHeatVary(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_RadHeatConst(imin:imax, jmin:jmax, kmin:kmax) )

    
    ! 温位の放射強制項. 
    ! 地表面から RadHeight で指定された高度までの間で一様放射冷却を与える. 
    do k = kmin, kmax
      if ( z_Z(k) <= HeightDown  ) then
        xyz_RadHeatConst(:,:,k) = 0.0d0 
      elseif( z_Z(k) >= HeightUp ) then
        xyz_RadHeatConst(:,:,k) = 0.0d0 
      else
        xyz_RadHeatConst(:,:,k) = RadHeatRate / DayTime / xyz_ExnerBZ(:,:,k)
      end if
    end do
    
    ! 温位の放射強制項.
    ! 地表面から 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_RadHeatVary(:,:,k) = RadHeatRate / DayTime / xyz_ExnerBZ(:,:,k)
        
      elseif ( z_Z(k) > HeightDown .AND. z_Z(k) <= HeightUP) then
        xyz_RadHeatVary(:,:,k) =  &
          & RadHeatRate / DayTime / xyz_ExnerBZ(:,:,k) &
          & * (HeightUP - z_Z(k)) / (HeightUP - HeightDown) 
        
      else if (z_Z(k) > HeightUP) then
        xyz_RadHeatVary(:,:,k) = 0.0d0
      end if
    end do
     

    ! 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')

  end subroutine Radiation_Simple_init

end module Radiation_Simple