!= Cho and Arking (1981) による短波放射モデル
!
!= Short wave radiation model described by Lacis and Hansen (1981)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: radiation_CA81.f90,v 1.1 2010-01-11 01:28:10 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module radiation_CA81
  !
  != Cho and Arking (1981) による短波放射モデル
  !
  != Short wave radiation model described by Lacis and Hansen (1974)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 短波放射モデル.
  !
  ! This is a model of short wave radiation. 
  !
  !== References
  !
  !  Cho, M.-D., and A. Arking, An efficient method for computing the absorption of 
  !    solar radiation by water vapor, J. Atmos. Sci., 38, 798-807, 1981.
  !
  !== 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.

  ! Declaration statements
  !
  implicit none
  private

  ! 
  ! Public procedure
  !
  public :: RadiationCA81Flux
  public :: RadiationCA81NumKDFBin
  public :: RadiationCA81H2ODelOptDep


  real(DP), save     :: H2OScaleIndex

  real(DP), save     :: DiffFactor

  real(DP), save     :: ShortAtmosAlbedo
                              ! 大気アルベド.
                              ! Albedo of air.
  logical , save     :: FlagSimpleMagFac
                              ! エアマス計算のためのフラグ
                              ! Flag for air-mass calculation
  logical , save     :: radiation_ca81_inited

  data radiation_ca81_inited /.false./

  integer , parameter:: nlogkmax = 30
!!$  integer , parameter:: nbmax    = 5
  real(DP), parameter:: kdfdlogk = 0.3d0
  real(DP), save     :: a_kdflogk(1:nlogkmax)
  real(DP), save     :: a_kdfh   (1:nlogkmax)
  real(DP), save     :: a_kdfk   (1:nlogkmax)
!!$  real(DP), save     :: aa_kdfhi  (1:nlogkmax,1:nbmax)

  ! Solar constant (radiative energy integrated over spectrum, W m-2) used for 
  ! Cho and Arking (1981).
  ! The solar flux data used by Cho and Arking (1981) is that presented by 
  ! Labs and Neckel (1968). The total solar radiative energy integrated over spectrum 
  ! presented by Labs and Neckel (1968) is 1.958 cal cm-2 min-1.
  ! This value is 1361 W m-2.
  ! Labs and Neckel (1981) paper can be viewed in 
  ! http://articles.adsabs.harvard.edu//full/1968ZA.....69....1L/0000001.000.html
  !
  real(DP), parameter :: solarconstCA81       = 1361.0d0

  ! See caption of TABLE 3 of Cho and Arking (1981).
  !
  real(DP), parameter :: totalfluxoverH2Oband = 552.17d0

  ! MEMO: Total wavenumber range is 2600 - 12040 cm-1
  !
  ! Unit of k in log(k) in this table is g-1 cm2.
  data a_kdflogk &
    & / -5.0d0,  -4.7d0,  -4.4d0,  -4.1d0,  -3.8d0,  &
    &   -3.5d0,  -3.2d0,  -2.9d0,  -2.6d0,  -2.3d0,  &
    &   -2.0d0,  -1.7d0,  -1.4d0,  -1.1d0,  -0.8d0,  &
    &   -0.5d0,  -0.2d0,   0.1d0,   0.4d0,   0.7d0,  &
    &    1.0d0,   1.3d0,   1.6d0,   1.9d0,   2.2d0,  &
    &    2.5d0,   2.8d0,   3.1d0,   3.4d0,   3.7d0   /

  ! Unit of h in this table is mW cm-2. 
  data a_kdfh &
    & / 25.032d0,  5.153d0,  6.198d0,  6.979d0,  7.023d0,  &
    &    6.762d0,  7.721d0,  8.118d0,  9.338d0, 10.287d0,  &
    &    9.902d0,  9.864d0,  9.798d0,  9.654d0,  9.020d0,  &
    &    8.061d0,  7.022d0,  6.259d0,  5.271d0,  4.439d0,  &
    &    3.507d0,  2.671d0,  2.000d0,  1.401d0,  0.957d0,  &
    &    0.638d0,  0.394d0,  0.268d0,  0.157d0,  0.176d0   /

!!$  ! Unit of h in this table is mW cm-2. 
!!$  data aa_kdfhi &
!!$    & /  8.605d0,  2.053d0,  2.380d0,  2.515d0,  1.985d0,  &
!!$    &    1.618d0,  1.615d0,  1.756d0,  2.471d0,  3.071d0,  &
!!$    &    3.346d0,  3.557d0,  3.278d0,  2.886d0,  2.397d0,  &
!!$    &    1.723d0,  1.223d0,  0.841d0,  0.558d0,  0.385d0,  &
!!$    &    0.237d0,  0.137d0,  0.078d0,  0.044d0,  0.022d0,  &
!!$    &    0.002d0,  0.0d0  ,  0.0d0  ,  0.0d0  ,  0.0d0  ,  & !  0.94 micron
!!$    !---
!!$    &    4.021d0,  0.432d0,  0.596d0,  0.598d0,  0.769d0,  &
!!$    &    0.830d0,  1.155d0,  1.832d0,  2.442d0,  2.686d0,  &
!!$    &    2.524d0,  2.501d0,  2.478d0,  2.478d0,  2.264d0,  &
!!$    &    1.856d0,  1.325d0,  0.910d0,  0.635d0,  0.434d0,  &
!!$    &    0.270d0,  0.166d0,  0.110d0,  0.081d0,  0.053d0,  &
!!$    &    0.024d0,  0.004d0,  0.0d0  ,  0,0d0  ,  0.0d0  ,  & !  1.14 micron
!!$    !---
!!$    &    5.246d0,  1.295d0,  1.440d0,  1.498d0,  1.533d0,  &
!!$    &    1.562d0,  1.573d0,  1.425d0,  1.478d0,  1.623d0,  &
!!$    &    1.366d0,  1.447d0,  1.838d0,  2.097d0,  2.340d0,  &
!!$    &    2.468d0,  2.428d0,  2.446d0,  2.112d0,  1.828d0,  &
!!$    &    1.420d0,  1.010d0,  0.701d0,  0.455d0,  0.297d0,  &
!!$    &    0.193d0,  0.112d0,  0.080d0,  0.034d0,  0.018d0,  & !  1.38 micron
!!$    !---
!!$    &    5.692d0,  0.764d0.  0.962d0,  1.168d0,  1.337d0,  &
!!$    &    1.493d0,  1.9350d,  1.722d0,  1.500d0,  1.523d0,  &
!!$    &    1.406d0,  1.197d0,  1.051d0,  1.051d0,  0.916d0,  &
!!$    &    1.003d0,  1.197d0,  1.288d0,  1.249d0,  1.054d0,  &
!!$    &    0.802d0,  0.613d0,  0.441d0,  0.285d0,  0.202d0,  &
!!$    &    0.133d0,  0.091d0,  0.064d0,  0.040d0,  0.024d0,  & !  1.87 micron
!!$    !---
!!$    &    0.081d0,  0.077d0,  0.172d0,  0.187d0,  0.236d0,  &
!!$    &    0.273d0,  0.318d0,  0.480d0,  0.589d0,  0.774d0,  &
!!$    &    0.835d0,  0.845d0,  0.933d0,  0.979d0,  0.998d0,  &
!!$    &    0.945d0,  0.823d0,  0.762d0,  0.710d0,  0.737d0,  &
!!$    &    0.777d0,  0.744d0,  0.670d0,  0.536d0,  0.383d0,  &
!!$    &    0.286d0,  0.187d0,  0.124d0,  0.084d0,  0.134     / !  2.7 micron



  character(*), parameter:: module_name = 'radiation_CA81'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: radiation_CA81.f90,v 1.1 2010-01-11 01:28:10 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  subroutine RadiationCA81Flux( &
    & xy_SurfAlbedo, xy_InAngle, xyz_QVap, xyr_Press, xyz_Press, & ! (in )
    & xyr_RadSFlux                                               & ! (out)
    & )


    ! USE statements
    !

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

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav, & ! $ g $ [m s-2].
                                 ! 
                                 ! Gravitational acceleration
                         PI      ! $ \pi $ .
                                 ! Circular constant

    real(DP), intent(in ):: xy_SurfAlbedo   (0:imax-1, 1:jmax)
    real(DP), intent(in ):: xy_InAngle      (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ):: xyz_QVap        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyr_Press       (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ):: xyz_Press       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyr_RadSFlux    (0:imax-1, 1:jmax, 0:kmax)


    !
    ! Work variables
    !
    real(DP):: RefPress

    real(DP):: xyz_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyr_H2OColDen   (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_H2OAbsAmt   (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyra_TransH2O   (0:imax-1, 1:jmax, 0:kmax, 1:nlogkmax)

    real(DP):: xy_MagFac       (0:imax-1, 1:jmax)
    real(DP):: xy_cosSZA       (0:imax-1, 1:jmax)

!!$    logical :: flag_dry_atmosphere

    integer :: i
    integer :: j
    integer :: k
    integer :: m


    if ( .not. radiation_ca81_inited ) then
      call RadiationCA81Init
    end if


!!$    ! Check for dry atmosphere
!!$    !
!!$    if ( all( xyz_QVap <= 0.0d0 ) ) then
!!$      flag_dry_atmosphere = .true.
!!$      write( 6, * ) 'Dry atmosphere'
!!$    else
!!$      flag_dry_atmosphere = .false.
!!$    end if


!!$    ! 大気アルベドの考慮
!!$    ! Taking atmospheric albedo into consideration
!!$    !
!!$    xy_IncomRadSFlux = xy_IncomRadSFlux * ( 1.0d0 - ShortAtmosAlbedo )

!!$    if ( flag_dry_atmosphere ) then
!!$      do k = 0, kmax
!!$        xyr_RadSFlux(:,:,k) = - xy_IncomRadSFlux(:,:) + ...
!!$      end do
!!$      return
!!$    end if


    RefPress = 300.0d2

    do k = 1, kmax
      xyz_H2ODelAbsAmt(:,:,k) = &
        &   ( xyz_Press(:,:,k) / RefPress )**H2OScaleIndex &
        & * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    xyr_H2OColDen(:,:,:) = 0.0d0
    do k = kmax-1, 0, -1
      xyr_H2OColDen(:,:,k) = xyr_H2OColDen(:,:,k+1) + xyz_H2ODelAbsAmt(:,:,k+1)
    end do


    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_InAngle(i,j) > 0.0d0 ) then
          xy_MagFac(i,j) = xy_InAngle(i,j)
          xy_cosSZA(i,j) = 1.0d0 / xy_InAngle(i,j)
        else
          xy_MagFac(i,j) = 0.0d0
          xy_cosSZA(i,j) = 0.0d0
        end if
      end do
    end do


    ! Calculation of flux
    !
    xyr_RadSFlux(:,:,:) = 0.0d0


    ! Downward flux
    !
    do k = 0, kmax
      xyr_H2OAbsAmt(:,:,k) = &
        &   xyr_H2OColDen(:,:,k) * xy_MagFac(:,:)
    end do
    do m = 1, nlogkmax
      xyra_TransH2O(:,:,:,m) = exp( - a_kdfk(m) * xyr_H2OAbsAmt )
    end do
    do m = 1, nlogkmax
      do k = 0, kmax
        xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k)            &
          & - xyra_TransH2O(:,:,k,m) * a_kdfh(m) * kdfdlogk  &
          &   * xy_cosSZA(:,:)
      end do
    end do


    ! Upward flux
    !
    do k = 0, kmax
      xyr_H2OAbsAmt(:,:,k) = &
        &   xyr_H2OColDen(:,:,0) * xy_MagFac(:,:) &
        & + ( xyr_H2OColDen(:,:,0) - xyr_H2OColDen(:,:,k) ) * DiffFactor
    end do
    do m = 1, nlogkmax
      xyra_TransH2O(:,:,:,m) = exp( - a_kdfk(m) * xyr_H2OAbsAmt )
    end do
    do m = 1, nlogkmax
      do k = 0, kmax
        xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k)            &
          & + xyra_TransH2O(:,:,k,m) * a_kdfh(m) * kdfdlogk  &
          &   * xy_cosSZA(:,:)                               &
          &   * xy_SurfAlbedo(:,:)
      end do
    end do


!!$    write( 6, * ) '***************************'
!!$    write( 6, * ) '***************************'
!!$    write( 6, * ) '***************************'
!!$    write( 6, * ) 'Short wave radiation out of H2O band is not added in radiation_CA81 module.'
!!$    write( 6, * ) '***************************'
!!$    write( 6, * ) '***************************'
!!$    write( 6, * ) '***************************'
!!$    ! Add flux over wavenumber range except for H2O band, that is treated by 
!!$    ! Cho and Arking (1981) scheme. 
!!$    !
!!$    do k = 0, kmax
!!$      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) &
!!$        & - ( solarconstCA81 - totalfluxoverH2Oband ) * xy_cosSZA(:,:)
!!$    end do

    ! 大気アルベドの考慮
    ! Taking atmospheric albedo into consideration
    !
    xyr_RadSFlux = xyr_RadSFlux * ( 1.0d0 - ShortAtmosAlbedo )


  end subroutine RadiationCA81Flux

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

  subroutine RadiationCA81NumKDFBin( nbin )

    integer, intent(out) :: nbin

    if ( .not. radiation_ca81_inited ) then
      call RadiationCA81Init
    end if

    nbin = nlogkmax

  end subroutine RadiationCA81NumKDFBin

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

  subroutine RadiationCA81H2ODelOptDep(        &
    & xyz_QVap, xyr_Press, xyz_Press, ikdfbin, & ! (in )
    & xyz_H2ODelOptDep, KDFWeight              & ! (out)
    & )


    ! USE statements
    !

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

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav, & ! $ g $ [m s-2].
                                 ! 
                                 ! Gravitational acceleration
                         PI      ! $ \pi $ .
                                 ! Circular constant

    real(DP), intent(in ):: xyz_QVap        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyr_Press       (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ):: xyz_Press       (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ):: ikdfbin
    real(DP), intent(out):: xyz_H2ODelOptDep(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: KDFWeight


    !
    ! Work variables
    !
    real(DP):: RefPress

    real(DP):: xyz_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax)

!!$    logical :: flag_dry_atmosphere

    integer :: k
    integer :: m


    if ( .not. radiation_ca81_inited ) then
      call RadiationCA81Init
    end if


!!$    ! Check for dry atmosphere
!!$    !
!!$    if ( all( xyz_QVap <= 0.0d0 ) ) then
!!$      flag_dry_atmosphere = .true.
!!$      write( 6, * ) 'Dry atmosphere'
!!$    else
!!$      flag_dry_atmosphere = .false.
!!$    end if


    RefPress = 300.0d2

    do k = 1, kmax
      xyz_H2ODelAbsAmt(:,:,k) = &
        &   ( xyz_Press(:,:,k) / RefPress )**H2OScaleIndex &
        & * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    m = ikdfbin
    xyz_H2ODelOptDep = xyz_H2ODelAbsAmt * a_kdfk( m )

    KDFWeight = a_kdfh(m) * kdfdlogk


  end subroutine RadiationCA81H2ODelOptDep

  !**************************************************************************

  subroutine RadiationCA81Init


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

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

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

!!$    ! ヒストリデータ出力
!!$    ! 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 /radiation_CA81_nml/ &
      & DiffFactor, ShortAtmosAlbedo
          !
          ! デフォルト値については初期化手続 "radiation_LH74#RadiationLH74Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_LH74#RadiationLH74Init" for the default values.
          !

    ! デフォルト値の設定
    ! Default values settings
    !
    DiffFactor       = 1.66d0

    ShortAtmosAlbedo = 0.2d0

    ! 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 = radiation_CA81_nml, & ! (out)
        & iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    H2OScaleIndex = 0.8d0


    ! Unit is changed of k from g-1 cm2 to kg-1 m2.
    !
!!$    a_kdflogk = a_kdflogk + log( 1.0d3 * 1.0d-4 )
    a_kdflogk = a_kdflogk + log10( 1.0d3 * 1.0d-4 )

    ! Unit is changed from mW cm-2 to W m-2.
    !
    a_kdfh    = a_kdfh    * 1.0d-3 * 1.0d4

    ! Calculation of k of table of k-distribution fnction
    !
!!$    a_kdfk    = exp( a_kdflogk )
    a_kdfk    = 10.0d0**a_kdflogk

!!$    ! Unit is changed from mW cm-2 to W m-2.
!!$    !
!!$    aa_kdfhi  = aa_kdfhi  * 1.0d0-3 * 1.0d4

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


    radiation_ca81_inited = .true.
  end subroutine RadiationCA81Init


end module radiation_CA81
