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

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

  !モジュール読み込み
  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, dz     !Z 座標軸(スカラー格子点)
  use constants,only: CpDry       ! 1 日の長さ [s]
  use basicset, only: xyz_ExnerBZ, & !エクスナー関数の基本場
    &                 xyz_DensBZ     !基本場の密度
  use DExnerDt, only: xyz_DExnerDt_xyz
  use namelist_util, only: namelist_filename
  
  !暗黙の型宣言禁止
  implicit none

  !private 属性
  private

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

  character(*), parameter:: module_name = 'radiation_bs1998'
                              ! モジュールの名称.
                              ! Module name
  public Radiation_bs1998_init

contains

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

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer  :: k             !ループ変数

    allocate( xyz_PTempRadBS1998(imin:imax, jmin:jmax, kmin:kmax) )
    allocate( xyz_ExnerRadBS1998(imin:imax, jmin:jmax, kmin:kmax) )
    
    ! 温位の放射強制項.
    ! BS1998 で用いられている放射プロファイルを用いる
    !
    do k = kmin, kmax
      xyz_PTempRadBS1998(:,:,k) = cal_Qsub( z_Z(k) ) / CpDry / xyz_DensBZ(:,:,k) / xyz_ExnerBZ(:,:,k)
    end do
    xyz_ExnerRadBS1998 = xyz_DExnerDt_xyz( xyz_PTempRadBS1998 ) 

    ! Output
    !
    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='K.s-1',    &
      & xtype='double')

  end subroutine Radiation_bs1998_init


  function cal_Qsub( z0 )
    ! 同じルーチンが surfaceflux_bs1998.f90 にも含まれているので注意すること.

    
    implicit none

    real(DP), intent(in) :: z0
    real(DP)             :: cal_Qsub
    real(DP), parameter  :: z_Upper = 6.70d4      ! 上部頂点高度
    real(DP), parameter  :: z_Lower = 2.70d4      ! 下部頂点高度
    real(DP), parameter  :: c_Upper = 2.70d-2     ! 第 1 係数
    real(DP), parameter  :: c_Lower = 3.60d-3     ! 第 2 係数
    real(DP), parameter  :: s_Upper = 7.50d3      ! 第 1 標準偏差
    real(DP), parameter  :: s_Lower = 1.30d4      ! 第 2 標準偏差

    cal_Qsub =                                                                      &
      & (                                                                           &
      &     c_Lower * exp( - ( z0 - z_Lower )**2.0d0 / ( 2.0d0 * s_Lower **2.0d0) ) &
      &   + c_Upper * exp( - ( z0 - z_Upper )**2.0d0 / ( 2.0d0 * s_Upper **2.0d0) ) &
      &  ) 
    
    
  end function cal_Qsub


end module Radiation_bs1998
