!= Schneider and Liu (2009) の放射モデル
!
!= Radiation model by Schneider and Liu (2009)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: radiation_SL09.f90,v 1.3 2010-10-07 15:43:03 yot Exp $
! Tag Name::  $Name: dcpam5-20101015 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module radiation_SL09
  !
  != Schneider and Liu (2009) の放射モデル
  !
  != Radiation model by Schneider and Liu (2009)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! This is a radiation model described by Schneider and Liu (2009). 
  !
  !== References
  !
  !  Schneider, T. and J. Liu, 
  !    Formation of jets and equatorial superrotation on Jupiter, 
  !    J. Atmos. Sci., 69, 579, 2009.
  !
  !== 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_SL09_nml
  !

  ! USE statements
  !

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

  ! 格子点設定
  ! 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

  implicit none

  private


  ! Private variables
  !
  real(DP), save :: SWOptDepAtRefPress
  real(DP), save :: SWRefPress
  real(DP), save :: SWOrd
  real(DP), save :: LWOptDepAtRefPress
  real(DP), save :: LWRefPress
  real(DP), save :: LWOrd
  real(DP), save :: SolarConst


  ! 公開変数
  ! Public variables
  !
  logical, save, public:: radiation_SL09_inited = .false.
                              ! 初期設定フラグ.
                              ! Initialization flag

  public :: RadiationSL09Flux

  character(*), parameter:: module_name = 'radiation_SL09'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20101015 $' // &
    & '$Id: radiation_SL09.f90,v 1.3 2010-10-07 15:43:03 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

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

contains

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

  subroutine RadiationSL09Flux(                     &
    & xyr_Press, xyz_Temp,                          &
    & xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux  &
    & )


    ! USE statements
    !

    !
    ! Physical constants settings
    !
    use constants, only: &
      &                  PI       ! $ \pi $ .
                                  ! Circular constant

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only : y_Lat

    !
    ! Solve radiative transfer equation in two stream approximation
    !
    use radiation_two_stream_app, only: RadiationTwoStreamApp, IDSchemeShortWave, RadiationTwoStreamAppHomogAtm

    ! 放射関連ルーチン
    ! Routines for radiation calculation
    !
    use radiation_utils, only : RadiationRTEQNonScat

    real(DP), intent(in ) :: xyr_Press    ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(in ) :: xyz_Temp     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out) :: xyr_RadSFlux ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadLFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


    ! Work variables
    !
    real(DP) :: SolarFluxTOA
!!$    real(DP) :: QeRatio
!!$    real(DP) :: xyz_SSA      (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP) :: xyz_AF       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_SurfTemp  (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfAlbedo(0:imax-1, 1:jmax)
    real(DP) :: xy_InAngle   (0:imax-1, 1:jmax)
    real(DP) :: xy_CosSZA    (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfEmis  (0:imax-1, 1:jmax)
    real(DP) :: xyr_OptDep   (0:imax-1, 1:jmax, 0:kmax)

    real(DP) :: SSA
    real(DP) :: AF

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


    ! 初期化
    ! Initialization
    !
    if ( .not. radiation_SL09_inited ) call RadiationSL09Init


    ! Short wave radiation
    !
    xyr_OptDep = SWOptDepAtRefPress * ( xyr_Press / SWRefPress )**SWOrd

    SolarFluxTOA  = SolarConst / PI


    SSA = 0.8_DP
    AF  = 0.204_DP
    !   Af = 0 may be much better than 0.204 when Eddington approximation is used.
!!$    AF         = 0.0_DP


    do j = 1, jmax
      xy_CosSZA(:,j) = cos( y_Lat(j) )
    end do
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_CosSZA(i,j) > 0.0_DP ) then
          xy_InAngle(i,j) = 1.0_DP / xy_CosSZA(i,j)
        else
          xy_InAngle(i,j) = 0.0_DP
        end if
      end do
    end do

    !   Unused variable but this is required as an argument
    !
    xy_SurfAlbedo = 1.0d100

    call RadiationTwoStreamAppHomogAtm(                            &
      & xy_SurfAlbedo, SolarFluxTOA, xy_InAngle, SSA, AF, xyr_OptDep,  & ! (in )
      & xyr_RadSFlux,                                                    & ! (out)
      & FlagSemiInfAtm = .true., FlagSL09 = .true.                    & ! (in ) optional
      & )



    ! Long wave radiation
    !

    !   Although the surface temperature and surface emissivity are set, but are not used.
    !
    xy_SurfTemp = 1.0d100
    xy_SurfEmis = 1.0d100

    xyr_OptDep = LWOptDepAtRefPress * ( xyr_Press / LWRefPress )**LWOrd


    call RadiationRTEQNonScat(                    &
      & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyr_OptDep, & ! (in )
      & xyr_RadLFlux, xyra_DelRadLFlux,                 & ! (out)
      & xy_SurfUpRadLFluxBase = xyr_RadSFlux(:,:,0)     & ! (in ) optional
      & )


  end subroutine RadiationSL09Flux

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

  subroutine RadiationSL09SWFlux(                               &
    & SolarFluxTOA, xy_CosSZA, SSA, AF, xyr_OptDep, &
    & xyr_RadSFlux                                              &
    & )

    real(DP), intent(in ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_CosSZA(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: SSA
    real(DP), intent(in ) :: AF
    real(DP), intent(in ) :: xyr_OptDep   ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadSFlux ( 0:imax-1, 1:jmax, 0:kmax )

    ! Work variables
    !
    real(DP) :: BondAlbedo
    real(DP) :: Gamma
    integer  :: j, k


    BondAlbedo = &
      &   ( sqrt( 1.0_DP - SSA * AF ) - sqrt( 1.0_DP - SSA ) ) &
      & / ( sqrt( 1.0_DP - SSA * AF ) + sqrt( 1.0_DP - SSA ) )

    Gamma = 2.0_DP * sqrt( 1.0_DP - SSA ) * sqrt( 1.0_DP - SSA * AF )

    do k = 0, kmax
      do j = 1, jmax
        xyr_RadSFlux(:,j,k) =                   &
          & - SolarFluxTOA * xy_CosSZA(:,j) &
          & * ( 1.0_DP - BondAlbedo ) * exp( -Gamma * xyr_OptDep(:,j,k) )
      end do
    end do


  end subroutine RadiationSL09SWFlux

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

  subroutine RadiationSL09Init

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

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

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


    ! 宣言文 ; Declaration statements
    !

    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_SL09_nml/ &
      & SWOptDepAtRefPress, SWRefPress, SWOrd, &
      & LWOptDepAtRefPress, LWRefPress, LWOrd, &
      & SolarConst
          !
          ! デフォルト値については初期化手続 "radiation_SL09#RadiationSL09Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_SL09#RadiationSL09Init" for the default values.
          !


    ! デフォルト値の設定
    ! Default values settings
    !
    SWOptDepAtRefPress =  3.0_DP
    SWRefPress         =  3.0d5
    SWOrd              =  1.0_DP

    LWOptDepAtRefPress = 80.0_DP
    LWRefPress         =  3.0d5
    LWOrd              =  2.0_DP

    SolarConst         = 50.7_DP




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

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



    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'SWOptDepAtRefPress = %f', d = (/ SWOptDepAtRefPress /) )
    call MessageNotify( 'M', module_name, 'SWRefPress         = %f', d = (/ SWRefPress /) )
    call MessageNotify( 'M', module_name, 'SWOrd              = %f', d = (/ SWOrd /) )
    call MessageNotify( 'M', module_name, 'LWOptDepAtRefPress = %f', d = (/ LWOptDepAtRefPress /) )
    call MessageNotify( 'M', module_name, 'LWRefPress         = %f', d = (/ LWRefPress /) )
    call MessageNotify( 'M', module_name, 'LWOrd              = %f', d = (/ LWOrd /) )
    call MessageNotify( 'M', module_name, 'SolarConst         = %f', d = (/ SolarConst /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    radiation_SL09_inited = .true.

  end subroutine RadiationSL09Init

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

end module radiation_SL09
