!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2007. All rights reserved.
!---------------------------------------------------------------------

module physics_radiation_incoming_sr_mod
!= 物理過程 入射放射計算(固定日射分布)モジュール
!
!== 概要
!
!  * 式自体は, 日変化・季節変化有りの日射分布を計算するもの.
!    しかし, 現状では日付けを取得できないので
!    ソースコードにハードコードされた日付け・時刻の日射分布しか
!    計算できない.
!    同期回転惑星大気の計算はできるはず.
!
!== History
!   2007-05-02 Masaki Ishiwatari
!
!== TODO 
! * (2007-5-2 石渡)
!   * 日付け・時刻を問合せできるようにしなければならない.
!   * ソースドキュメントは英語で書くのだっけ?
!   * パラメータ値を可変に
!   * モジュール名長すぎる. phys_ にするという案が出ていたかな?
!   * 太陽定数は SolarConst とか
!   * x_Lon を追加した. interface 考え直すべき.
!     xy_InAngle (AGCM5 の SECZET) を出力する必要性は?
!   * 変数説明には次元まで書いた方が良いのでは?

  use type_mod,    only : REKIND, DBKIND, INTKIND, TOKEN, STRING

  implicit none

  private
  public :: physics_radiation_incoming

contains

  subroutine physics_radiation_incoming( &
    & xy_IncomRadSFlux        , & ! (out) 長波フラックス
    & xy_InAngle              , & ! (out) sec(入射角)
    & x_Lon                   , & !(in) 経度
    & y_Lat                    ) !(in) 緯度

    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use nmlfile_mod, only: nmlfile_init, nmlfile_open, nmlfile_close
    use constants_mod, only: PI  
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
    use dc_message, only: MessageNotify

    implicit none

    real(DBKIND), intent(out) :: xy_IncomRadSFlux(im,jm) ! 短波フラックス
    real(DBKIND), intent(out) :: xy_InAngle(im,jm) ! sec(入射角)
    real(DBKIND), intent(in) :: x_Lon(im) ! 経度
    real(DBKIND), intent(in) :: y_Lat(jm) ! 緯度

    character(STRING),  parameter:: subname = "physics_radiation_incoming_sr"
    integer(INTKIND)    :: i, j
 
    real(DBKIND) :: SolarCoeff = 1380.0d0 ! 太陽定数 [W/m^2]
    real(DBKIND) :: AtmosAlbedo = 0.0d0 ! 大気アルベド
    real(DBKIND) :: LatSubSolar = 180.0d0 ! 太陽直下点緯度 (単位:degree)
    real(DBKIND) :: LonSubSolar = 0.0d0 ! 太陽直下点経度 (単位:degree)
    real(DBKIND) :: SinDel ! 太陽入射角
    real(DBKIND) :: ANGHR ! 時角
    real(DBKIND) :: COSZET ! cos(入射角)

    logical                     :: nmlreadable
    integer(INTKIND)            :: nmlunit, nmlstat

    namelist /physics_radiation_incoming_sr_nml/ &
      & SolarCoeff, &
      & AtmosAlbedo, &
      & LonSubSolar, &
      & LatSubSolar

    continue

    call BeginSub(subname)
  
    ! read physics_radiation_incomfing_sr_nml
    call nmlfile_init
    call nmlfile_open(nmlunit, nmlreadable)
    if (nmlreadable) then
       read(nmlunit, nml=physics_radiation_incoming_sr_nml, iostat=nmlstat)
       call DbgMessage('Stat of NAMELIST physics_radiation_incoming_sr_nml Input is <%d>', &
            &           i=(/nmlstat/))
       write(0, nml=physics_radiation_incoming_sr_nml)
    else
       call DbgMessage('Not Read NAMELIST physics_radiation_incoming_sr_nml')
       call MessageNotify('W', subname, &
            & 'Can not Read NAMELIST physics_radiation_incoming_sr_nml. Force Use Default Value.')
    end if
    call nmlfile_close

    SinDel =  sin( LatSubSolar*PI/180.0d0  )

    do i = 1, im
      do j = 1, jm
        ANGHR  = (LonSubSolar - x_Lon(i))*PI/180.0d0
        COSZET = sin(y_Lat(j)*PI/180.0d0)*SinDel &
          & + cos(y_Lat(j)*PI/180.0d0)* sqrt(1.0d0 - SinDel**2) * cos(ANGHR)

        IF ( COSZET .GT. 0. ) THEN
          xy_IncomRadSFlux (i,j) = - SolarCoeff*(1.0d0 - AtmosAlbedo)*COSZET
          xy_InAngle (i,j) = 1.0d0 / COSZET
        ELSE
          xy_IncomRadSFlux  (i,j) = 0.0d0
          xy_InAngle (i,j) = 0.0d0
        ENDIF
      end do
    end do

    call EndSub(subname)

  end subroutine physics_radiation_incoming

end module physics_radiation_incoming_sr_mod