!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! physics_radiation_long.f90 
!
! History
!   2005/09/21 Yamada Yukiko     create
!

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) 温度 (整数)
       & xy_SurfTemp          , & ! (in) 地表面温度
       & xyr_TauQvap           , & ! (in) 光学的厚さ：水
       & xyr_TauDryAir           ) ! (in) 光学的厚さ：空気

    !==== Dependency
    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 dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Output
    !
    real(DBKIND), intent(out) :: &
       & xyr_RadLFlux(im,jm,km+1)          , & ! (out) 長波フラックス
       & xyro_DelRadLFlux(im,jm,km+1,0:1)      ! (out) 長波地表温度変化

    !==== Input
    !
    real(DBKIND), intent(in) :: &
       & xyz_Temp(im,jm,km)           , & ! (in) 温度 (整数)
       & xy_SurfTemp(im,jm)           , & ! (in) 地表面温度
       & xyr_TauQvap(im,jm,km+1)      , & ! (in) 光学的厚さ：水
       & xyr_TauDryAir(im,jm,km+1)        ! (in) 光学的厚さ：空気

    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "physics_radiation_long"

    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)    :: i, j, k

    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

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

    continue


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

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

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
















