!= Chou and Lee (1996) による短波放射モデル
!
!= Short wave radiation model described by Chou and Lee (1996)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: rad_CL1996.f90,v 1.4 2011/06/19 11:12:42 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module rad_CL1996
  !
  != Chou and Lee (1996) による短波放射モデル
  !
  != Short wave radiation model described by Chou and Lee (1996)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 短波放射モデル.
  !
  ! This is a model of short wave radiation. 
  !
  !== References
  !
  !  Chou, M.-D., and K.-T. Lee, 
  !    Parameterizations for the absorption of solar radiation by water vapor and ozone, 
  !    J. Atmos. Sci., 53, 1203-1208, 1996.
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
!!$  ! NAMELIST#radiation_DennouAGCM_nml
  !

  ! USE statements
  !

  ! 
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 
  ! Grid points settings
  !
  use gridset, only: imax, & ! 
                             ! Number of grid points in longitude
    &                jmax, & ! 
                             ! Number of grid points in latitude
    &                kmax    ! 
                             ! Number of vertical level

  ! Declaration statements
  !
  implicit none
  private

  !
  ! Public procedure
  !
  public :: RadCL1996NumBands
  public :: RadCL1996UVVISParams
  public :: RadCL1996IRH2ONumKDFBin
  public :: RadCL1996IRH2OKDFParams
  public :: RadCL1996ScaleH2OVapMass
  public :: RadCL1996Init



  integer , parameter :: nband1 = 8  ! * 14500 to 57143 cm-1 (0.175 to 0.70 micron)
  integer , parameter :: nband2 = 3  ! *  2600 to 14500 cm-1 (0.70-10 micron)
  integer , parameter :: nkdf     = 10

  real(DP), save      :: a_UVVIFracSolarFlux(1:nband1)
  real(DP), save      :: a_UVVISO3AbsCoef   (1:nband1)
  real(DP), save      :: a_UVVISRayScatCoef (1:nband1)

  real(DP), save      :: a_IRH2Okdfk   (1:nkdf)                        ! k
  real(DP), save      :: aa_IRH2Okdfdgi(1:nkdf,nband1+1:nband1+nband2) ! k dist. func.



  logical , save     :: rad_cl1996_inited

  data rad_cl1996_inited /.false./




  ! Table 1
  ! Unit of k is g-1 cm2.
  data a_IRH2Okdfk &
    & / &
    & 0.0010 , 0.0133 , 0.0422 , 0.1334 , 0.4217 , 1.334  , 5.623  , 31.62  , 177.8  , 1000.0    &
    & /

  ! Table 1
  data aa_IRH2Okdfdgi &
    & / &
    & 0.20673, 0.03497, 0.03011, 0.02260, 0.01336, 0.00696, 0.00441, 0.00115, 0.00026, 0.00000,  &  ! 0.7 -1.22
    & 0.08236, 0.01157, 0.01133, 0.01143, 0.01240, 0.01258, 0.01381, 0.00650, 0.00244, 0.00094,  &  ! 1.22-2.27
    & 0.01074, 0.00360, 0.00411, 0.00421, 0.00389, 0.00326, 0.00499, 0.00465, 0.00245, 0.00145   &  ! 2.27-10
    & /


!!$  data aa_KDFParams &
!!$    & / &
!!$    !  (0.7-10)  
!!$    &  0.29983, 0.05014, 0.04555, 0.03824, 0.02965, 0.02280, 0.02321, 0.01230, 0.00515, 0.00239   &
!!$    & /


!!$  data aa_UVVISBandParams &
!!$    & / &
!!$    &   , & ! 0.175-0.225              (UV-C)
!!$    &   , & ! 0.225-0.245, 0.260-0.280 (UV-C)
!!$    &   , & ! 0.245-0.260              (UV-C)
!!$    &   , & ! 0.280-0.295              (UV-B)
!!$    &   , & ! 0.295-0.310              (UV-B)
!!$    &   , & ! 0.310-0.320              (UV-B)
!!$    &   , & ! 0.320-0.400              (UV-A)
!!$    &     & ! 0.400-0.700              (PAR)
!!$    & /



  ! Table 3.
  data a_UVVIFracSolarFlux &
    & /    0.00057,   0.00367,   0.00083,   0.00417,   0.00600,   0.00556,   0.05913,    0.39081 /

  ! Table 3.
  ! Unit is (cm-atm)_{STP}^{-1}
  data a_UVVISO3AbsCoef &
    & /   30.47   , 187.24   , 301.92   ,  42.83   ,   7.09   ,   1.25   ,   0.0345 ,    0.0539  /

  ! Table 3.
  ! Unit is (mb)^{-1}
  data a_UVVISRayScatCoef &
    & /    0.00604,   0.00170,   0.00222,   0.00132,   0.00107,   0.00091,   0.00055,    0.00012 /



  character(*), parameter:: module_name = 'rad_CL1996'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_CL1996.f90,v 1.4 2011/06/19 11:12:42 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  !--------------------------------------------------------------------------------------

  subroutine RadCL1996NumBands( &
    & nbands1, nbands2 & ! (out)
    & )

    integer, intent(out) :: nbands1
    integer, intent(out) :: nbands2


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_cl1996_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    nbands1 = nband1
    nbands2 = nband2

  end subroutine RadCL1996NumBands

  !--------------------------------------------------------------------------------------

  subroutine RadCL1996IRH2ONumKDFBin( &
    & nbin & ! (out)
    & )

    integer, intent(out) :: nbin


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_cl1996_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    nbin = nkdf

  end subroutine RadCL1996IRH2ONumKDFBin

  !--------------------------------------------------------------------------------------

  subroutine RadCL1996ScaleH2OVapMass(  &
    & xyz_Temp, xyz_DelH2OVapMass, xyz_Press, & ! (in )
    & xyz_DelH2OVapMassScaled                 & ! (out)
    & )


    ! USE statements
    !

    real(DP), intent(in ):: xyz_Temp               (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelH2OVapMass      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_Press              (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyz_DelH2OVapMassScaled(0:imax-1, 1:jmax, 1:kmax)


    !
    ! Work variables
    !
    real(DP), parameter :: H2OScaleIndex = 0.8_DP
    real(DP), parameter :: RefPress      = 300.0d2
    real(DP), parameter :: RefTemp       = 240.0d0


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_cl1996_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    xyz_DelH2OVapMassScaled =                                         &
      &   ( xyz_Press / RefPress )**H2OScaleIndex                     &
      & * exp( 0.00135_DP * ( xyz_Temp - RefTemp ) )                  &
      & * xyz_DelH2OVapMass


  end subroutine RadCL1996ScaleH2OVapMass

  !--------------------------------------------------------------------------------------

  subroutine RadCL1996IRH2OKDFParams(        &
    & iband, ikdfbin,                              & ! (in )
    & KDFAbsCoef, KDFWeight                        & ! (out)
    & )


    ! USE statements
    !

    integer , intent(in ):: iband
    integer , intent(in ):: ikdfbin
    real(DP), intent(out):: KDFAbsCoef
    real(DP), intent(out):: KDFWeight


    !
    ! Work variables
    !
    integer :: l
    integer :: m


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_cl1996_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    l = iband
    m = ikdfbin

    KDFAbsCoef = a_IRH2Okdfk   (m)

    KDFWeight  = aa_IRH2Okdfdgi(m,l)


  end subroutine RadCL1996IRH2OKDFParams

  !--------------------------------------------------------------------------------------

  subroutine RadCL1996UVVISParams(                    &
    & iband,                                                & ! (in )
    & UVVISFracSolarFlux, UVVISO3AbsCoef, UVVISRayScatCoef  & ! (out)
    & )


    ! USE statements
    !

    integer , intent(in ):: iband
    real(DP), intent(out):: UVVISFracSolarFlux
    real(DP), intent(out):: UVVISO3AbsCoef
    real(DP), intent(out):: UVVISRayScatCoef


    !
    ! Work variables
    !
    integer :: l


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_cl1996_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    l = iband

    UVVISFracSolarFlux = a_UVVIFracSolarFlux(l)
    UVVISO3AbsCoef     = a_UVVISO3AbsCoef   (l)
    UVVISRayScatCoef   = a_UVVISRayScatCoef (l)


  end subroutine RadCL1996UVVISParams

  !--------------------------------------------------------------------------------------

  subroutine RadCL1996Init


!!$    ! NAMELIST ファイル入力に関するユーティリティ
!!$    ! Utilities for NAMELIST file input
!!$    !
!!$    use namelist_util, only: NmlutilMsg

!!$    ! ファイル入出力補助
!!$    ! File I/O support
!!$    !
!!$    use dc_iounit, only: FileOpen


!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    use gtool_historyauto, only: HistoryAutoAddVariable


!!$    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
!!$                              ! Unit number for NAMELIST file open
!!$    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
!!$                              ! IOSTAT of NAMELIST read

!!$    ! NAMELIST 変数群
!!$    ! NAMELIST group name
!!$    !
!!$    namelist /rad_CL1996_nml/ !&
!!$      & ShortAtmosAlbedo
!!$          !
!!$          ! デフォルト値については初期化手続 "rad_CL1996#RadCL1996Init"
!!$          ! のソースコードを参照のこと.
!!$          !
!!$          ! Refer to source codes in the initialization procedure
!!$          ! "rad_LH74#RadLH74Init" for the default values.
!!$          !



    if ( rad_cl1996_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !

!!$    ! NAMELIST の読み込み
!!$    ! NAMELIST is input
!!$    !
!!$    if ( trim(namelist_filename) /= '' ) then
!!$      call FileOpen( unit_nml, &          ! (out)
!!$        & namelist_filename, mode = 'r' ) ! (in)
!!$
!!$      rewind( unit_nml )
!!$      read( unit_nml,                     & ! (in)
!!$        & nml = rad_CL1996_nml,     & ! (out)
!!$        & iostat = iostat_nml )             ! (out)
!!$      close( unit_nml )
!!$
!!$      call NmlutilMsg( iostat_nml, module_name ) ! (in)
!!$    end if



    ! Unit is changed of k from g-1 cm2 to kg-1 m2.
    !
    a_IRH2Okdfk = a_IRH2Okdfk * 1.0d3 * 1.0d-4

    ! Convert unit of O3 absorption coefficient from (cm-atm)^{-1} to (kg m-2)^{-1}
    !   In order to convert unit from cm-1 atm-1 to m2 kg-1, multiply by
    !     1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 273.15d0.
    !
    a_UVVISO3AbsCoef = a_UVVISO3AbsCoef &
        & * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 273.15d0

    ! Convert unit of Rayleigh scattering coefficient from (mb)^{-1} to (kg m-2)^{-1}
    !   memo:
    !     1    mbar = 1e2 Pa
    !     1e-2 mbar = 1   Pa
    !
    ! In order to convert unit from (mb)^{-1} to (kg m-2)^{-1}, scattering 
    ! coefficient is multiplied by 1.0d-1 * g, where g is the gravitational 
    ! acceleration. In the old version of the code, a variable Grav in constants 
    ! module is used for gravitational acceleration, as is shown in a line below. 
    ! However, this variable is replaced with a constant value 9.8_DP. 
    ! This is because the value of Grav may be changed to a value which is 
    ! different from the Earth's value. (Of course, the scattering coefficient 
    ! in unit of (kg m-2)^{-1} is independent on gravitational acceleration.)
    !
!!$    a_UVVISRayScatCoef = a_UVVISRayScatCoef * 1.0d-2 * Grav
    a_UVVISRayScatCoef = a_UVVISRayScatCoef * 1.0d-2 * 9.8_DP



    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
!!$    call MessageNotify( 'M', module_name, 'ShortAtmosAlbedo = %f', d = (/ ShortAtmosAlbedo /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    rad_cl1996_inited = .true.

  end subroutine RadCL1996Init

  !--------------------------------------------------------------------------------------

end module rad_CL1996
