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

module physics_radiation_short_mod

  implicit none

  private
  public :: physics_radiation_short

contains

    subroutine physics_radiation_short( &
         & xyr_RadSFlux         , & ! (inout) 長波フラックス
         & xyr_TauQvap          , & ! (in) 光学的厚さ：水
         & xyr_TauDryAir        , & ! (in) 光学的厚さ：空気
         & xy_InAngle           , & ! (in) sec(入射角)
         & xy_SurfAlbedo          ) ! (in) 地表アルベド

    !==== Dependency
    use type_mod,    only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod, only: im, jm, km
    use constants_mod, only: PI  
    use dc_trace,    only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump

    implicit none

    !==== Output
    !
    real(DBKIND), intent(inout) :: &
         & xyr_RadSFlux(im,jm,km+1) ! (inout) 短波フラックス

    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & xyr_TauQvap(im,jm,km+1)   , & ! (in) 光学的厚さ：水
         & xyr_TauDryAir(im,jm,km+1) , & ! (in) 光学的厚さ：空気
         & xy_InAngle(im,jm)         , & ! (in) sec(入射角)
         & xy_SurfAlbedo(im,jm)          ! (in) 地表アルベド

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

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

    integer(INTKIND), parameter :: BandNumber = 1     ! 長波バンド数
    real(DBKIND) :: &
         & AbsorpCoeffQvap(BandNumber)    , &         ! 水の吸収係数
         & AbsorpCoeffDryAir(BandNumber)  , &         ! 空気の吸収係数
         & BandWeight(BandNumber)         , &         ! バンドウェイト
         & SECSCT                                     !  散乱のsecζ
    Data AbsorpCoeffQvap / 0.002  /
    Data AbsorpCoeffDryAir / 0.0  /
    Data BandWeight  / 1.0   /
    Data SECSCT / 1.66 /
    
    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


    ! ----------------------------------------------------------
    xyr_RadSFlux(:,:,1:km) = 0.0d0

    do bn = 1, BandNumber
       do k = 1, km+1    
          ! ----  2. 各レベルでの下向き透過 ----
          
          if (k .NE. km+1) then 
             
             xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) & 
                  & +  BandWeight(bn) * xyr_RadSFlux(:,:,km+1) &
                  &  * EXP( - xy_InAngle(:,:) &
                  &      * ( AbsorpCoeffQvap(bn) * xyr_TauQvap(:,:,k) &
                  &      + AbsorpCoeffDryAir(bn) * xyr_TauDryAir(:,:,k) ) &
                  &       )
          end if
          
          ! ----  3. 各レベルでの上向き透過 ----
          xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) & 
               & -  BandWeight(bn) * xyr_RadSFlux(:,:,km+1) &
               &  * EXP( - xy_InAngle(:,:) &
               &      * ( AbsorpCoeffQvap(bn) * xyr_TauQvap(:,:,1) &
               &      + AbsorpCoeffDryAir(bn) * xyr_TauDryAir(:,:,1) ) &
               &       ) &
               &  * xy_SurfAlbedo(:,:) &                        
               &  * EXP( - SECSCT &
               &      * ( AbsorpCoeffQvap(bn) * ( xyr_TauQvap(:,:,1) &
               &                                - xyr_TauQvap(:,:,k) ) &
               &      + AbsorpCoeffDryAir(bn) * ( xyr_TauDryAir(:,:,1) &
               &                               - xyr_TauDryAir(:,:,k) ) &
               &       ) ) 
       end do
    end do
    ! ----------------------------------------------------------

    ! ----  4. 吸収なしのとき ----

    if ( BandNumber .EQ. 0 ) then 

       do k = 1, km+1
          xyr_RadSFlux(:,:,k) = (1.0d0 - xy_SurfAlbedo(:,:) ) &
               &               * xyr_RadSFlux(:,:,km+1) 
       end do

    end if


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

  end subroutine physics_radiation_short


end module physics_radiation_short_mod
















