!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_radiation_long.f90 
!
! == TODO
!
!  BandNumber も入力できるようにしなくては!
! ==History
!   2005/09/21 Yamada Yukiko     create
!   2007-05-12 M. Ishiwatari interface 変更.  xyr_Temp を追加


module physics_radiation_long_mod

  implicit none

  private
  public :: physics_radiation_long

contains

  subroutine physics_radiation_long( &
    & xyr_RadLFlux          , & ! (out) 長波フラックス
    & xyro_DelRadLFlux      , & ! (out) 長波地表温度変化
    & xyz_Temp              , & ! (in) 温度 (整数)
    & xyr_Temp              , & ! (in) 温度 (半整数)
    & xy_SurfTemp          , & ! (in) 地表面温度
    & xyr_TauQvap           , & ! (in) 光学的厚さ：水
    & xyr_TauDryAir           ) ! (in) 光学的厚さ：空気

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

    implicit none
    real(DBKIND), intent(out) :: xyr_RadLFlux(im,jm,km+1) ! 長波フラックス
    real(DBKIND), intent(out) :: xyro_DelRadLFlux(im,jm,km+1,0:1)  
                                             ! 長波地表温度変化
    real(DBKIND), intent(in) :: xyz_Temp(im,jm,km) ! 温度 (整数)
    real(DBKIND), intent(in) :: xyr_Temp(im,jm,km+1) ! 温度 (半整数)
    real(DBKIND), intent(in) :: xy_SurfTemp(im,jm) ! 地表面温度
    real(DBKIND), intent(in) :: xyr_TauQvap(im,jm,km+1) !  光学的厚さ：水
    real(DBKIND), intent(in) :: xyr_TauDryAir(im,jm,km+1) ! 光学的厚さ：空気

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_radiation_long"
    integer(INTKIND)    :: i, j, k
            ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    real(DBKIND) :: &
         & xyr_Trans(im,jm,km+1)     , & ! 透過係数
         & xyr_Trans1(im,jm,km+1)    , & !  1/2 レベルからの透過係数
         & xyr_Trans2(im,jm,km+1)    , & ! 3/2 レベルからの透過係数
         & xyz_PiB(im,jm,km  )       , &  ! πＢ＝σＴ**4
         & xy_SurfPiB(im,jm)              ! 地表のπＢ

    integer(INTKIND)    :: kk , bn
    real(DBKIND)        :: BandWeightSum

!  (2007-5-20 石渡) ちょっとテストのため バンド数を 1 にしてみる.

    integer(INTKIND), parameter :: BandNumber = 4    ! 長波バンド数
!    integer(INTKIND), parameter :: BandNumber = 1    ! 長波バンド数
    real(DBKIND) :: AbsorpCoeffQvap(BandNumber) ! 水の吸収係数
    real(DBKIND) :: AbsorpCoeffDryAir(BandNumber) ! 空気の吸収係数
    real(DBKIND) :: BandWeight(BandNumber) ! バンドウェイト
    real(DBKIND) :: PathLengthFact     ! 光路長のファクタ

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

!  (2007-5-20 石渡) agcm5.3 のデフォルト放射(4 色バンドモデル)
    data AbsorpCoeffQvap / 8.0  , 1.0  , 0.1  , 0.0  /
    data AbsorpCoeffDryAir  / 0.0  , 0.0  , 0.0  , 5.E-5 /
    data BandWeight / 0.2  , 0.1  , 0.1  , 0.6 /
    PathLengthFact = 1.5 

    namelist /physics_radiation_long_nml/ &
      & AbsorpCoeffQvap, &
      & AbsorpCoeffDryAir, &
      & BandWeight, &
      & PathLengthFact

    continue

    !----------------------------------------------------------------
    !   開始処理
    !----------------------------------------------------------------
    call BeginSub(subname)

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


    !----------------------------------------------------------------
    !   放射計算
    !----------------------------------------------------------------

    ! ---- 1. バンドウェイトの設定 ----

    BandWeightSum = 0.0d0

    do bn = 1, BandNumber
       BandWeightSum =  BandWeightSum +  BandWeight(bn)
    end do

    do bn = 1, BandNumber
       BandWeight(bn) = BandWeight(bn) / BandWeightSum
    end do

    ! ---- 2. πＢの計算 ----
    xyz_PiB(:,:,:) = StB * ( xyz_Temp(:,:,:)**4 )
    xy_SurfPiB(:,:)  = StB * ( xy_SurfTemp(:,:)**4 )


    ! --------------------------------------------------
    do k = 1, km+1

       ! ---- 3. 透過関数計算 ----

       xyr_Trans = 0.0d0
       
       do bn = 1, BandNumber
          do kk = 1, km+1
             xyr_Trans(:,:,kk) = &
                  & xyr_Trans(:,:,kk) &
                  & + BandWeight(bn) &
                  & * EXP( - PathLengthFact &
                  &        * ( AbsorpCoeffQvap(bn) &
                  &            * ABS( xyr_TauQvap(:,:,kk) &
                  &                 - xyr_TauQvap(:,:,k)  ) &
                  &          + AbsorpCoeffDryAir(bn) &
                  &            * ABS( xyr_TauDryAir(:,:,kk) &
                  &                 - xyr_TauDryAir(:,:,k)  ) ) )
          end do
       end do
       

       ! ---- 4. 放射フラックス計算 ----
       xyr_RadLFlux(:,:,k) = xy_SurfPiB(:,:) * xyr_Trans(:,:,1) 
 
       do kk = 1, km
          xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) &
               & - xyz_PiB(:,:,kk) &
               &   * ( xyr_Trans(:,:,kk) - xyr_Trans(:,:,kk+1) ) 
       end do

    ! ---- 5.  補正項計算用透過関数  ----
       xyr_Trans1(:,:,k) = xyr_Trans(:,:,1)
       xyr_Trans2(:,:,k) = xyr_Trans(:,:,2)

    end do
    ! --------------------------------------------------

    ! ---- 6. 地表面補正用 ----
    do k = 1, km+1
       xyro_DelRadLFlux(:,:,k,0) = 4.0d0 * xy_SurfPiB(:,:) &
            & / xy_SurfTemp(:,:) * xyr_Trans1(:,:,k)
       
       xyro_DelRadLFlux(:,:,k,1) = 4.0d0 * xyz_PiB(:,:,1) &
            & / xyz_Temp(:,:,1) * ( xyr_Trans2(:,:,k) - xyr_Trans1(:,:,k) )
    end do


    !----------------------------------------------------------------
    !   終了処理
    !----------------------------------------------------------------
    call EndSub(subname)

  end subroutine physics_radiation_long


end module physics_radiation_long_mod