!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! = Description
!         以下では多色に対応できてない.
!        水蒸気・乾燥空気それぞれ1 色まで.
!   正味フラックスだけではなく上向きフラックスと下向きフッラックスに
!   分けて出力したい場合もある.
!   agcm 5 の時は ENTRY で対処していた.
!
! == TODO
!
!  BandNumber も入力できるようにしなくては!
! ==History
!   2007-05-12 M. Ishiwatari agcm5 ソースからport

module physics_radiation_long_runaway_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)     ! 透過係数
    real(DBKIND) :: xyr_Trans1(im,jm,km+1)    ! 1/2 レベルからの透過係数
    real(DBKIND) :: xyr_Trans2(im,jm,km+1)    ! 3/2 レベルからの透過係数
    real(DBKIND) :: xyr_PiB(im,jm,km+1)      ! πＢ＝σＴ**4
    real(DBKIND) :: xy_SurfPiB(im,jm)         ! 地表のπＢ
    real(DBKIND) :: xyr_OptDepth(im,jm,km+1)  ! 光学的厚さ
    real(DBKIND) :: xyr_RadLUpFlux(im,jm,km+1)  ! 上向き長波フラックス
    real(DBKIND) :: xyr_RadLDownFlux(im,jm,km+1)  ! 下向き長波フラックス

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

    integer(INTKIND), parameter :: BandNumber = 4    ! 長波バンド数
    real(DBKIND) :: AbsorpCoeffQvap(BandNumber) ! 水の吸収係数
    real(DBKIND) :: AbsorpCoeffDryAir(BandNumber) ! 空気の吸収係数
    real(DBKIND) :: BandWeight(BandNumber) ! バンドウェイト
    real(DBKIND) :: PathLengthFact = 1.5    ! 光路長のファクタ
    real(DBKIND) :: CritDeltaOptDepth = 0.1 ! スキームを切替えるΔτ

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

    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 /

    namelist /physics_radiation_long_runaway_nml/ &
      & AbsorpCoeffQvap, &
      & AbsorpCoeffDryAir, &
      & BandWeight, &
      & PathLengthFact, &
      & CritDeltaOptDepth

    continue

    ! 開始処理
    call BeginSub(subname)

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

    ! バンドウェイトの設定
    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_PiB(:,:,:) = StB * ( xyr_Temp(:,:,:)**4 )
    xy_SurfPiB(:,:)  = StB * ( xy_SurfTemp(:,:)**4 )

    ! 光学的厚さ計算
    xyr_OptDepth = 0.0d0
    do bn = 1, BandNumber
      do k = 1, km+1
        xyr_OptDepth(:,:,k) = &
          & PathLengthFact &
          & *(   AbsorpCoeffQvap(bn) * xyr_TauQvap(:,:,k) &
          &    + AbsorpCoeffDryAir(bn) * xyr_TauDryAir(:,:,k) &
          &  )
      end do
    end do

    xyr_RadLUpFlux = 0.0d0
    xyr_RadLDownFlux = 0.0d0

    do k = 1, km+1

      ! 透過関数計算
      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

      ! 放射フラックス計算
      xyr_RadLUpFlux (:,:,k) = &
        &   xyr_PiB (:,:,k) &
        & + ( xy_SurfPiB(:,:) - xyr_PiB(:,:,1) ) &
        &   * xyr_Trans (:,:,1)
      xyr_RadLDownFlux (:,:,k) = &
        & - xyr_PiB (:,:,km+1) * xyr_Trans (:,:,km+1) &
        & + xyr_PiB (:,:,k)

      do kk = 2,k
        do j = 1,jm
          do i = 1,im
            if(ABS(xyr_OptDepth(i,j,kk) - xyr_OptDepth(i,j,kk-1) ) &
              & .GT. CritDeltaOptDepth ) then
              xyr_RadLUpFlux(i,j,k) = &
                &   xyr_RadLUpFlux(i,j,k) &
                & - (xyr_PiB(i,j,kk-1) - xyr_PiB(i,j,kk)) &
                &   /(xyr_OptDepth(i,j,kk-1) - xyr_OptDepth(i,j,kk)) &
                &   *(xyr_Trans(i,j,kk-1) - xyr_Trans(i,j,kk)) 
            else
              xyr_RadLUpFlux(i,j,k) = &
                &   xyr_RadLUpFlux(i,j,k) &
                & + (xyr_PiB(i,j,kk-1) - xyr_PiB(i,j,kk) ) &
                &   *(xyr_Trans(i,j,kk) + xyr_Trans (i,j,kk-1))/2.0d0
            end if
          end do
        end do
      end do

      do kk = k+1, km+1
        do j = 1,jm
          do i = 1,im
            if( ABS(xyr_OptDepth(i,j,kk) - xyr_OptDepth(i,j,kk-1))  &
              & .GT. CritDeltaOptDepth ) THEN
              xyr_RadLDownFlux(i,j,k) = &
                &   xyr_RadLDownFlux(i,j,k) &
                & - (xyr_PiB(i,j,kk-1) - xyr_PiB(i,j,kk)) &
                &   /(xyr_OptDepth(i,j,kk-1) - xyr_OptDepth(i,j,kk)) &
                &   *(xyr_Trans (i,j,kk-1) - xyr_Trans (i,j,kk)) 
            else
              xyr_RadLDownFlux(i,j,k) = &
                &   xyr_RadLDownFlux(i,j,k) &
                & - (xyr_PiB(i,j,kk-1) - xyr_PiB(i,j,kk)) &
                &   *(xyr_Trans(i,j,kk) + xyr_Trans(i,j,kk-1))/2.0d0
            end if
          end do
        end do
      end do

      xyr_RadLFlux(:,:,k) = xyr_RadLUpFlux(:,:,k) - xyr_RadLDownFlux(:,:,k) 

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

    end do

    ! 地表面補正用
    do k = 1, km+1
      do j = 1,jm
        do i = 1,im
          xyro_DelRadLFlux(i,j,k,0) = &
            & 4.0d0 * xy_SurfPiB(i,j)/xy_SurfTemp(i,j) &
            & *xyr_Trans1(i,j,k)
          xyro_DelRadLFlux(i,j,k,1) = &
            & -4.0d0 * xyr_PiB(i,j,1) / xyr_Temp(i,j,1) &
            &  *xyr_Trans1(i,j,k)
          if( ABS(xyr_OptDepth(i,j,2) - xyr_OptDepth(i,j,1) ) &
            & .GT. CritDeltaOptDepth ) then
            xyro_DelRadLFlux(i,j,k,1) = &
              &   xyro_DelRadLFlux(i,j,k,1) &
              & + xyr_PiB(i,j,1) / xyr_Temp(i,j,1) &
              &   /(xyr_OptDepth(i,j,1) - xyr_OptDepth(i,j,2)) &
              &   *( xyr_Trans1(i,j,k) - xyr_Trans2(i,j,k) ) 
          else 
            xyro_DelRadLFlux(i,j,k,1) = &
              &   xyro_DelRadLFlux(i,j,k,1) &
              & + xyr_PiB (i,j,1)/xyr_Temp(i,j,1) &
              &   *(xyr_Trans1(i,j,k)+xyr_Trans2(i,j,k))/2.0d0
          end if
        end do 
      end do 
    end do 

    !  終了処理
    call EndSub(subname)

  end subroutine physics_radiation_long

end module physics_radiation_long_runaway_mod
















