!= 謾ｾ蟆�ヵ繝ｩ繝�け繧ｹ (邁｡蜊倬聞豕｢繝舌Φ繝峨Δ繝�Ν)
!
!= Radiation flux (simple longwave band model)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: rad_simple_LW.f90,v 1.3 2013/05/25 06:47:33 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module rad_simple_LW
  !
  != 謾ｾ蟆�ヵ繝ｩ繝�け繧ｹ (邁｡蜊倬聞豕｢繝舌Φ繝峨Δ繝�Ν)
  !
  != Radiation flux (simple longwave band model)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 貂ｩ蠎ｦ, 豈疲ｹｿ, 豌怜悸縺九ｉ, 謾ｾ蟆�ヵ繝ｩ繝�け繧ｹ繧定ｨ育ｮ励☆繧区叛蟆�Δ繝�Ν縺ｧ縺�.
  !
  ! This is a radiation model that calculates radiation flux from
  ! temperature, specific humidity, and air pressure.
  !
  !== References
  !
  !
  !
  !== Procedures List
  !
  ! RadSimpleLWFlux     :: 謾ｾ蟆�ヵ繝ｩ繝�け繧ｹ縺ｮ險育ｮ�
!!$  ! RadSimpleLWFinalize :: 邨ゆｺ��逅� (繝｢繧ｸ繝･繝ｼ繝ｫ蜀�Κ縺ｮ螟画焚縺ｮ蜑ｲ繧贋ｻ倥￠隗｣髯､)
  ! ------------------- :: ------------
  ! RadSimpleLWFlux     :: Calculate radiation flux
!!$  ! RadSimpleLWFinalize :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
  ! NAMELIST#rad_Simple_LW_nml
  !

  ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
  !

  ! 譬ｼ蟄千せ險ｭ螳�
  ! Grid points settings
  !
  use gridset, only: imax, & ! 邨悟ｺｦ譬ｼ蟄千せ謨ｰ.
                             ! Number of grid points in longitude
    &                jmax, & ! 邱ｯ蠎ｦ譬ｼ蟄千せ謨ｰ.
                             ! Number of grid points in latitude
    &                kmax    ! 驩帷峩螻､謨ｰ.
                             ! Number of vertical level

  ! 遞ｮ蛻･蝙九ヱ繝ｩ繝｡繧ｿ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 蛟咲ｲｾ蠎ｦ螳滓焚蝙�. Double precision.
    &                 STRING, &  ! 譁�ｭ怜�.       Strings.
    &                 TOKEN      ! 繧ｭ繝ｼ繝ｯ繝ｼ繝�.   Keywords.

  ! NAMELIST 繝輔ぃ繧､繝ｫ蜈･蜉帙↓髢｢縺吶ｋ繝ｦ繝ｼ繝�ぅ繝ｪ繝�ぅ
  ! Utilities for NAMELIST file input
  !
  use namelist_util, only: MaxNmlArySize
                              ! NAMELIST 縺九ｉ隱ｭ縺ｿ霎ｼ繧驟榊�縺ｮ譛螟ｧ繧ｵ繧､繧ｺ.
                              ! Maximum size of arrays loaded from NAMELIST

  ! 繝｡繝�そ繝ｼ繧ｸ蜃ｺ蜉�
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 繝ｪ繧ｹ繧ｿ繝ｼ繝医ョ繝ｼ繧ｿ蜃ｺ蜉�
  ! Restart data output
  !
  use gtool_history, only: GT_HISTORY

  ! 螳｣險譁� ; Declaration statements
  !
  implicit none
  private

  ! 蜈ｬ髢区焔邯壹″
  ! Public procedure
  !
  public:: RadSimpleLWFlux
  public:: RadSimpleLWInit

  ! 蜈ｬ髢句､画焚
  ! Public variables
  !
  logical, save :: rad_simple_LW_inited = .false.
                              ! 蛻晄悄險ｭ螳壹ヵ繝ｩ繧ｰ.
                              ! Initialization flag

  ! 髱槫�髢句､画焚
  ! Private variables
  !
  logical, save:: Old_Flux_saved = .false.
                              ! 荳蠎ｦ險育ｮ励＠縺溘ヵ繝ｩ繝�け繧ｹ繧剃ｿ晏ｭ倥＠縺溘％縺ｨ繧堤､ｺ縺吶ヵ繝ｩ繧ｰ.
                              ! Flag for saving of flux calculated once

  integer , save:: nbmax
                              ! 髟ｷ豕｢繝舌Φ繝画焚.
                              ! Number of long wave band
  real(DP), save:: a_WNBnds(0:MaxNmlArySize)
                              !
                              ! Wavenumber bounds for bands
  real(DP), save:: a_AbsCoefDryCom(1:MaxNmlArySize)
                              ! $ \bar{k}_R $ . 遨ｺ豌励�蜷ｸ蜿惹ｿよ焚
                              ! Absorption coefficient of dry component
  real(DP), save:: a_PressScaleIndDryCom(1:MaxNmlArySize)
                              !
                              ! Pressure scaling index for dry component
  real(DP), save:: a_RefPressDryCom(1:MaxNmlArySize)
                              !
                              ! Reference pressure for dry component.
  real(DP), save:: a_AbsCoefH2OVap(1:MaxNmlArySize)
                              ! $ k_R $ . 豌ｴ縺ｮ蜷ｸ蜿惹ｿよ焚
                              ! Absorption coefficient of water vapor
  real(DP), save:: a_PressScaleIndH2OVap(1:MaxNmlArySize)
                              !
                              ! Pressure scaling index for water vapor
  real(DP), save:: a_RefPressH2OVap(1:MaxNmlArySize)
                              !
                              ! Reference pressure for water vapor
  real(DP), save:: DiffFact
                              !
                              ! Diffusivity factor
  real(DP), save:: RadActDryComMMR
                              !
                              ! Mass mixing ratio of radiative active dry component

  integer , save:: NumGaussNode

  character(*), parameter:: module_name = 'rad_simple_LW'
                              ! 繝｢繧ｸ繝･繝ｼ繝ｫ縺ｮ蜷咲ｧｰ.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_simple_LW.f90,v 1.3 2013/05/25 06:47:33 yot Exp $'
                              ! 繝｢繧ｸ繝･繝ｼ繝ｫ縺ｮ繝舌�繧ｸ繝ｧ繝ｳ
                              ! Module version

contains

  !--------------------------------------------------------------------------------------

  subroutine RadSimpleLWFlux(                                         &
    & xy_SurfAlbedo, xy_SurfEmis,                                     & ! (in)
    & xyr_Press, xyz_Press, xyz_Temp, xy_SurfTemp,                    & ! (in)
    & xyz_DelAtmMass, xyz_QH2OVap,                                    & ! (in)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                                 & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                          & ! (out)
    & )
    !
    ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ縺ｮ險育ｮ�
    !
    ! Calculate long wave flux
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    ! 迚ｩ逅��謨ｰ蟄ｦ螳壽焚險ｭ螳�
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI,               &   ! $ \pi $ .
                              ! 蜀�捉邇�.  Circular constant
      & StB                   ! $ \sigma_{SB} $ .
                              ! 繧ｹ繝�ヵ繧｡繝ｳ繝懊Ν繝��繝ｳ螳壽焚.
                              ! Stefan-Boltzmann constant

    ! 繝励Λ繝ｳ繧ｯ髢｢謨ｰ縺ｮ險育ｮ�
    ! Calculate Planck function
    !
    use planck_func, only : Integ_PF_GQ_Array3D, Integ_PF_GQ_Array2D, Integ_DPFDT_GQ_Array2D

    ! 謨｣荵ｱ繧堤┌隕悶＠縺滓叛蟆�ｼ晞＃譁ｹ遞句ｼ�
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScatMonoSemiAnal

    ! 螳｣險譁� ; Declaration statements
    !
    real(DP), intent(in):: xy_SurfAlbedo     (0:imax-1, 1:jmax)
                              !
                              ! Surface albedo
    real(DP), intent(in):: xy_SurfEmis       (0:imax-1, 1:jmax)
                              ! 蝨ｰ陦ｨ髱｢蟆��邇�.
                              ! Surface emissivity
    real(DP), intent(in):: xyr_Press         (0:imax-1, 1:jmax, 0:kmax)
                              ! $ P $ .     蝨ｧ蜉�. Pressure
    real(DP), intent(in):: xyz_Press         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ P $ .     蝨ｧ蜉�. Pressure
    real(DP), intent(in):: xyz_Temp          (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xy_SurfTemp       (0:imax-1, 1:jmax)
                              ! 蝨ｰ陦ｨ髱｢貂ｩ蠎ｦ.
                              ! Surface temperature
    real(DP), intent(in):: xyz_DelAtmMass    (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Atmospheric mass of layers
    real(DP), intent(in):: xyz_QH2OVap       (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Specific humidity
    real(DP), intent(out):: xyr_RadLUwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Upward longwave flux
    real(DP), intent(out):: xyr_RadLDwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Downward longwave flux
    real(DP), intent(out):: xyra_DelRadLUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !
    real(DP), intent(out):: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !

    ! 菴懈･ｭ螟画焚
    ! Work variables
    !
    real(DP) :: xyr_Temp    (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_IntPF   (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xy_SurfIntPF(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT0(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT1(0:imax-1, 1:jmax)

    real(DP) :: xyz_DelOptDepDryCom(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelOptDepH2OVap(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyr_RadUwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Upward longwave flux
    real(DP) :: xyr_RadDwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Downward longwave flux
    real(DP) :: xyra_DelRadUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !
    real(DP) :: xyra_DelRadDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !

    real(DP) :: WNs
    real(DP) :: WNe

    integer:: k               ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 豕｢髟ｷ縺ｫ縺､縺�※蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in wavenumber bands

    ! 螳溯｡梧枚 ; Executable statement
    !

    k = 0
    xyr_Temp(:,:,k) =                                    &
      &      ( xyz_Temp (:,:,k+2) - xyz_Temp (:,:,k+1) ) &
      & / log( xyz_Press(:,:,k+2) / xyz_Press(:,:,k+1) ) &
      & * log( xyr_Press(:,:,k  ) / xyz_Press(:,:,k+1) ) &
      & + xyz_Temp(:,:,k+1)
    do k = 1, kmax-1
      xyr_Temp(:,:,k) =                                &
        &      ( xyz_Temp (:,:,k+1) - xyz_Temp (:,:,k) ) &
        & / log( xyz_Press(:,:,k+1) / xyz_Press(:,:,k) ) &
        & * log( xyr_Press(:,:,k  ) / xyz_Press(:,:,k) ) &
        & + xyz_Temp(:,:,k)
    end do
    k = kmax
    xyr_Temp(:,:,k) = xyz_Temp(:,:,k)

    !   Initialization
    !
    xyr_RadLUwFlux     = 0.0_DP
    xyr_RadLDwFlux     = 0.0_DP
    xyra_DelRadLUwFlux = 0.0_DP
    xyra_DelRadLDwFlux = 0.0_DP
    !
    LOOP_BAND_RTE : do n = 1, nbmax

      ! $ \pi B $, $ \pi DBDT $ 縺ｮ險育ｮ�
      ! Calculate $ \pi B $ and $ \pi DBDT $
      !
      if ( nbmax == 1 ) then

        xy_SurfIntPF = xy_SurfEmis * StB * ( xy_SurfTemp**4 )
        xyr_IntPF    =               StB * ( xyr_Temp**4 )
        xy_IntDPFDT0 = xy_SurfEmis * 4.0_DP * xy_SurfIntPF / xy_SurfTemp
!!$        xy_IntDPFDT1 =               4.0_DP * xyz_IntPF(:,:,1) / xyz_Temp(:,:,1)
        xy_IntDPFDT1 = 0.0_DP

      else

        WNs = a_WNBnds(n-1)
        WNe = a_WNBnds(n  )
        call Integ_PF_GQ_Array3D(        &
          & WNs, WNe, NumGaussNode,      &
          & 0, imax-1, 1, jmax, 0, kmax, &
          & xyr_Temp,                    &
          & xyr_IntPF                    &
          & )
        call Integ_PF_GQ_Array2D(   &
          & WNs, WNe, NumGaussNode, &
          & 0, imax-1, 1, jmax,     &
          & xy_SurfTemp,            &
          & xy_SurfIntPF            &
          & )
!!$        call Integ_DPFDT_GQ_Array2D(             &
!!$          & WNs, WNe, NumGaussNode,              & ! (in )
!!$          & 0, imax-1, 1, jmax, xyz_Temp(:,:,1), & ! (in )
!!$          & xy_IntDPFDT1                         & ! (out)
!!$          & )
        xy_IntDPFDT1 = 0.0_DP
        call Integ_DPFDT_GQ_Array2D(         &
          & WNs, WNe, NumGaussNode,          & ! (in )
          & 0, imax-1, 1, jmax, xy_SurfTemp, & ! (in )
          & xy_IntDPFDT0                     & ! (out)
          & )

        xy_SurfIntPF = xy_SurfEmis * PI * xy_SurfIntPF
        xyr_IntPF    =               PI * xyr_IntPF
        xy_IntDPFDT0 = xy_SurfEmis * PI * xy_IntDPFDT0
        xy_IntDPFDT1 =               PI * xy_IntDPFDT1

      end if

      ! 蜈牙ｭｦ逧�字縺輔�險育ｮ�
      ! Calculate optical depth
      !
      xyz_DelOptDepDryCom = a_AbsCoefDryCom(n)                             &
        & * ( xyz_Press / a_RefPressDryCom(n) )**a_PressScaleIndDryCom(n)  &
        & * xyz_DelAtmMass                                                 &
        & * RadActDryComMMR
      xyz_DelOptDepH2OVap = a_AbsCoefH2OVap(n)                             &
        & * ( xyz_Press / a_RefPressH2OVap(n) )**a_PressScaleIndH2OVap(n)  &
        & * xyz_DelAtmMass * xyz_QH2OVap

      call RadRTENonScatMonoSemiAnal(                              &
        & xy_SurfAlbedo,                                           & ! (in)
        & xyr_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0,     & ! (in)
        & ( xyz_DelOptDepDryCom + xyz_DelOptDepH2OVap ),           & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux,                            & ! (out)
        & xyra_DelRadUwFlux, xyra_DelRadDwFlux                     & ! (out)
        & )

      xyr_RadLUwFlux     = xyr_RadLUwFlux     + xyr_RadUwFlux
      xyr_RadLDwFlux     = xyr_RadLDwFlux     + xyr_RadDwFlux
      xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadUwFlux
      xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadDwFlux

    end do LOOP_BAND_RTE

  end subroutine RadSimpleLWFlux

  !--------------------------------------------------------------------------------------

  subroutine OLD_RadSimpleLWFlux(                                      &
    & xy_SurfEmis, xyz_Temp, xyz_Press, xy_SurfTemp,                   & ! (in)
    & xyz_DelAtmMass, xyz_QH2OVap,                                     & ! (in)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                                  & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                           & ! (out)
    & )
    !
    ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ縺ｮ險育ｮ�
    !
    ! Calculate long wave flux
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    ! 迚ｩ逅��謨ｰ蟄ｦ螳壽焚險ｭ螳�
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI,               &   ! $ \pi $ .
                              ! 蜀�捉邇�.  Circular constant
      & StB                   ! $ \sigma_{SB} $ .
                              ! 繧ｹ繝�ヵ繧｡繝ｳ繝懊Ν繝��繝ｳ螳壽焚.
                              ! Stefan-Boltzmann constant

    ! 繝励Λ繝ｳ繧ｯ髢｢謨ｰ縺ｮ險育ｮ�
    ! Calculate Planck function
    !
    use planck_func, only : Integ_PF_GQ_Array3D, Integ_PF_GQ_Array2D, Integ_DPFDT_GQ_Array2D

    ! 謨｣荵ｱ繧堤┌隕悶＠縺滓叛蟆�ｼ晞＃譁ｹ遞句ｼ�
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScat

    ! 螳｣險譁� ; Declaration statements
    !
    real(DP), intent(in):: xy_SurfEmis       (0:imax-1, 1:jmax)
                              ! 蝨ｰ陦ｨ髱｢蟆��邇�.
                              ! Surface emissivity
    real(DP), intent(in):: xyz_Temp          (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     貂ｩ蠎ｦ. Temperature
    real(DP), intent(in):: xyz_Press         (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     貂ｩ蠎ｦ. Temperature

    real(DP), intent(in):: xy_SurfTemp       (0:imax-1, 1:jmax)
                              ! 蝨ｰ陦ｨ髱｢貂ｩ蠎ｦ.
                              ! Surface temperature
    real(DP), intent(in):: xyz_DelAtmMass    (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Atmospheric mass of layers
    real(DP), intent(in):: xyz_QH2OVap       (0:imax-1, 1:jmax, 1:kmax)
                              !
                              ! Specific humidity
    real(DP), intent(out):: xyr_RadLUwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Upward longwave flux
    real(DP), intent(out):: xyr_RadLDwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Downward longwave flux
    real(DP), intent(out):: xyra_DelRadLUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !
    real(DP), intent(out):: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !

    ! 菴懈･ｭ螟画焚
    ! Work variables
    !
    real(DP) :: xyz_DelOptDepDryCom(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelOptDepH2OVap(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_TransEachLayer (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyrr_Trans         (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 騾城℃菫よ焚.
                              ! Transmission coefficient
    real(DP) :: xyz_IntPF   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_SurfIntPF(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT0(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT1(0:imax-1, 1:jmax)

    real(DP) :: xyr_RadUwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Upward longwave flux
    real(DP) :: xyr_RadDwFlux     (0:imax-1, 1:jmax, 0:kmax)
                              ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ.
                              ! Downward longwave flux
    real(DP) :: xyra_DelRadUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !
    real(DP) :: xyra_DelRadDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 髟ｷ豕｢蝨ｰ陦ｨ貂ｩ蠎ｦ螟牙喧.
                              !

    real(DP) :: WNs
    real(DP) :: WNe

    integer:: k, kk           ! 驩帷峩譁ｹ蜷代↓蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 豕｢髟ｷ縺ｫ縺､縺�※蝗槭ｋ DO 繝ｫ繝ｼ繝礼畑菴懈･ｭ螟画焚
                              ! Work variables for DO loop in wavenumber bands

    ! 螳溯｡梧枚 ; Executable statement
    !

    !   Initialization
    !
    xyr_RadLUwFlux     = 0.0_DP
    xyr_RadLDwFlux     = 0.0_DP
    xyra_DelRadLUwFlux = 0.0_DP
    xyra_DelRadLDwFlux = 0.0_DP
    !
    LOOP_BAND_RTE : do n = 1, nbmax

      ! $ \pi B $, $ \pi DBDT $ 縺ｮ險育ｮ�
      ! Calculate $ \pi B $ and $ \pi DBDT $
      !
      if ( nbmax == 1 ) then

        xy_SurfIntPF = xy_SurfEmis * StB * ( xy_SurfTemp**4 )
        xyz_IntPF    =               StB * ( xyz_Temp**4 )
        xy_IntDPFDT0 = xy_SurfEmis * 4.0_DP * xy_SurfIntPF / xy_SurfTemp
        xy_IntDPFDT1 =               4.0_DP * xyz_IntPF(:,:,1) / xyz_Temp(:,:,1)

      else

        WNs = a_WNBnds(n-1)
        WNe = a_WNBnds(n  )
        call Integ_PF_GQ_Array3D(        &
          & WNs, WNe, NumGaussNode,      &
          & 0, imax-1, 1, jmax, 1, kmax, &
          & xyz_Temp,                    &
          & xyz_IntPF                    &
          & )
        call Integ_PF_GQ_Array2D(   &
          & WNs, WNe, NumGaussNode, &
          & 0, imax-1, 1, jmax,     &
          & xy_SurfTemp,            &
          & xy_SurfIntPF            &
          & )
        call Integ_DPFDT_GQ_Array2D(             &
          & WNs, WNe, NumGaussNode,              & ! (in )
          & 0, imax-1, 1, jmax, xyz_Temp(:,:,1), & ! (in )
          & xy_IntDPFDT1                         & ! (out)
          & )
        call Integ_DPFDT_GQ_Array2D(         &
          & WNs, WNe, NumGaussNode,          & ! (in )
          & 0, imax-1, 1, jmax, xy_SurfTemp, & ! (in )
          & xy_IntDPFDT0                     & ! (out)
          & )

        xy_SurfIntPF = xy_SurfEmis * PI * xy_SurfIntPF
        xyz_IntPF    =               PI * xyz_IntPF
        xy_IntDPFDT0 = xy_SurfEmis * PI * xy_IntDPFDT0
        xy_IntDPFDT1 =               PI * xy_IntDPFDT1

      end if

      ! 蜈牙ｭｦ逧�字縺輔�險育ｮ�
      ! Calculate optical depth
      !
      xyz_DelOptDepDryCom = a_AbsCoefDryCom(n)                             &
        & * ( xyz_Press / a_RefPressDryCom(n) )**a_PressScaleIndDryCom(n)  &
        & * xyz_DelAtmMass
      xyz_DelOptDepH2OVap = a_AbsCoefH2OVap(n)                             &
        & * ( xyz_Press / a_RefPressH2OVap(n) )**a_PressScaleIndH2OVap(n)  &
        & * xyz_DelAtmMass * xyz_QH2OVap

      ! 騾城℃髢｢謨ｰ縺ｮ險育ｮ�
      ! Calculate transmission functions
      !
      xyz_TransEachLayer = exp( - DiffFact * ( xyz_DelOptDepDryCom + xyz_DelOptDepH2OVap ) )
      do k = 0, kmax
        do kk = k, k
          xyrr_Trans(:,:,k,kk) = 1.0d0
        end do
        do kk = k+1, kmax
          xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk-1) * xyz_TransEachLayer(:,:,kk)
        end do
        do kk = 0, k-1
          xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,kk,k)
        end do
      end do

      call RadRTENonScat(                                      &
        & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0, & ! (in)
        & xyrr_Trans,                                          & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux,                        & ! (out)
        & xyra_DelRadUwFlux, xyra_DelRadDwFlux                 & ! (out)
        & )

      xyr_RadLUwFlux     = xyr_RadLUwFlux     + xyr_RadUwFlux
      xyr_RadLDwFlux     = xyr_RadLDwFlux     + xyr_RadDwFlux
      xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadUwFlux
      xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadDwFlux

    end do LOOP_BAND_RTE

  end subroutine OLD_RadSimpleLWFlux

  !--------------------------------------------------------------------------------------

  subroutine RadSimpleLWInit
    !
    ! rad_simple_LW 繝｢繧ｸ繝･繝ｼ繝ｫ縺ｮ蛻晄悄蛹悶ｒ陦後＞縺ｾ縺�.
    ! NAMELIST#rad_simple_LW_nml 縺ｮ隱ｭ縺ｿ霎ｼ縺ｿ縺ｯ縺薙�謇狗ｶ壹″縺ｧ陦後ｏ繧後∪縺�.
    !
    ! "rad_simple_LW" module is initialized.
    ! "NAMELIST#rad_simple_LW_nml" is loaded in this procedure.
    !

    ! 繝｢繧ｸ繝･繝ｼ繝ｫ蠑慕畑 ; USE statements
    !

    ! 蜃ｺ蜉帙ヵ繧｡繧､繝ｫ縺ｮ蝓ｺ譛ｬ諠��ｱ
    ! Basic information for output files
    !
    use fileset, only: &
      & FileTitle, &
                              ! 蜃ｺ蜉帙ョ繝ｼ繧ｿ繝輔ぃ繧､繝ｫ縺ｮ陦ｨ鬘�.
                              ! Title of output data files
      & FileSource, &
                              ! 繝��繧ｿ繝輔ぃ繧､繝ｫ菴懈�縺ｮ謇区ｮｵ.
                              ! Source of data file
      & FileInstitution
                              ! 繝��繧ｿ繝輔ぃ繧､繝ｫ繧呈怙邨ら噪縺ｫ螟画峩縺励◆邨�ｹ�/蛟倶ｺｺ.
                              ! Institution or person that changes data files for the last time

    ! 迚ｩ逅��謨ｰ蟄ｦ螳壽焚險ｭ螳�
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI                    ! $ \pi $ .
                              ! 蜀�捉邇�.  Circular constant

    ! 蠎ｧ讓吶ョ繝ｼ繧ｿ險ｭ螳�
    ! Axes data settings
    !
    use axesset, only: &
      & x_Lon, &
                              ! $ \lambda $ [rad.] . 邨悟ｺｦ. Longitude
      & x_Lon_Weight, &
                              ! $ \Delta \lambda $ [rad.] .
                              ! 邨悟ｺｦ蠎ｧ讓咎㍾縺ｿ.
                              ! Weight of longitude
      & y_Lat, &
                              ! $ \varphi $ [rad.] . 邱ｯ蠎ｦ. Latitude
      & y_Lat_Weight, &
                              ! $ \Delta \varphi $ [rad.] .
                              ! 邱ｯ蠎ｦ蠎ｧ讓咎㍾縺ｿ.
                              ! Weight of latitude
      & z_Sigma, &
                              ! $ \sigma $ 繝ｬ繝吶Ν (謨ｴ謨ｰ).
                              ! Full $ \sigma $ level
      & r_Sigma, &
                              ! $ \sigma $ 繝ｬ繝吶Ν (蜊頑紛謨ｰ).
                              ! Half $ \sigma $ level
      & z_DelSigma
                              ! $ \Delta \sigma $ (謨ｴ謨ｰ).
                              ! $ \Delta \sigma $ (Full)

    ! 譎ょ綾邂｡逅�
    ! Time control
    !
    use timeset, only: &
      & RestartTime           ! 繝ｪ繧ｹ繧ｿ繝ｼ繝磯幕蟋区凾蛻ｻ.
                              ! Retart time of calculation

    ! NAMELIST 繝輔ぃ繧､繝ｫ蜈･蜉帙↓髢｢縺吶ｋ繝ｦ繝ｼ繝�ぅ繝ｪ繝�ぅ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! 證ｦ縺ｨ譌･譎ゅ�蜿悶ｊ謇ｱ縺�
    ! Calendar and Date handler
    !
    use dc_calendar, only: DCCalConvertByUnit

    ! 繝輔ぃ繧､繝ｫ蜈･蜃ｺ蜉幄｣懷勧
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 遞ｮ蛻･蝙九ヱ繝ｩ繝｡繧ｿ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 讓呎ｺ門�蜉帙�陬�ｽｮ逡ｪ蜿ｷ. Unit number of standard output

    ! 譁�ｭ怜�謫堺ｽ�
    ! Character handling
    !
    use dc_string, only: toChar

    ! 謨｣荵ｱ繧堤┌隕悶＠縺滓叛蟆�ｼ晞＃譁ｹ遞句ｼ�
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScatInit

    ! 螳｣險譁� ; Declaration statements
    !
    integer:: unit_nml        ! NAMELIST 繝輔ぃ繧､繝ｫ繧ｪ繝ｼ繝励Φ逕ｨ陬�ｽｮ逡ｪ蜿ｷ.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 隱ｭ縺ｿ霎ｼ縺ｿ譎ゅ� IOSTAT.
                              ! IOSTAT of NAMELIST read

    real(DP) :: WNBnds             (0:MaxNmlArySize)
    real(DP) :: AbsCoefDryCom      (1:MaxNmlArySize)
    real(DP) :: PressScaleIndDryCom(1:MaxNmlArySize)
    real(DP) :: RefPressDryCom     (1:MaxNmlArySize)
    real(DP) :: AbsCoefH2OVap      (1:MaxNmlArySize)
    real(DP) :: PressScaleIndH2OVap(1:MaxNmlArySize)
    real(DP) :: RefPressH2OVap     (1:MaxNmlArySize)

    integer  :: n

    ! NAMELIST 螟画焚鄒､
    ! NAMELIST group name
    !
    namelist /rad_simple_LW_nml/ &
      & nbmax, &
      & WNBnds,                                             &
      & AbsCoefDryCom, PressScaleIndDryCom, RefPressDryCom, &
      & AbsCoefH2OVap, PressScaleIndH2OVap, RefPressH2OVap, &
      & RadActDryComMMR,                                    &
      & DiffFact,                                           &
      & NumGaussNode
          !
          ! 繝�ヵ繧ｩ繝ｫ繝亥､縺ｫ縺､縺�※縺ｯ蛻晄悄蛹匁焔邯� "rad_DennouAGCM#RadInit"
          ! 縺ｮ繧ｽ繝ｼ繧ｹ繧ｳ繝ｼ繝峨ｒ蜿ら�縺ｮ縺薙→.
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_DennouAGCM#RadInit" for the default values.
          !

    ! 螳溯｡梧枚 ; Executable statement
    !

    if ( rad_simple_LW_inited ) return

    ! 繝�ヵ繧ｩ繝ｫ繝亥､縺ｮ險ｭ螳�
    ! Default values settings
    !

    ! 髟ｷ豕｢繝輔Λ繝�け繧ｹ逕ｨ諠��ｱ
    ! Information for long wave flux
    !

    nbmax                  = 1

    WNBnds                 = -999.9_DP
    AbsCoefDryCom          = -999.9_DP
    PressScaleIndDryCom    = -999.9_DP
    RefPressDryCom         = -999.9_DP
    AbsCoefH2OVap          = -999.9_DP
    PressScaleIndH2OVap    = -999.9_DP
    RefPressH2OVap         = -999.9_DP

    AbsCoefDryCom      (1:nbmax) = (/ 5.0d-5 /)
    PressScaleIndDryCom(1:nbmax) = (/ 0.0d0  /)
    RefPressDryCom     (1:nbmax) = (/ 1.0d5  /)

    AbsCoefH2OVap      (1:nbmax) = (/ 1.0d-2 /)
    PressScaleIndH2OVap(1:nbmax) = (/ 0.0d0  /)
    RefPressH2OVap     (1:nbmax) = (/ 1.0d5  /)

    RadActDryComMMR       = 1.0_DP

    DiffFact        = 1.66_DP

    NumGaussNode    = 5

    ! NAMELIST 縺ｮ隱ｭ縺ｿ霎ｼ縺ｿ
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml,                    & ! (in)
        & nml = rad_simple_LW_nml,       & ! (out)
        & iostat = iostat_nml )            ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if

    if ( nbmax > MaxNmlArySize ) then
      call MessageNotify( 'E', module_name, 'nbmax = %d > %d', i = (/ nbmax, MaxNmlArySize /) )
    end if

    a_WNBnds              = WNBnds * 100.0_DP    ! Convert from cm-1 to m-1
!!$    a_AbsCoefDryCom = 5.0d-5
    a_AbsCoefDryCom       = AbsCoefDryCom
    a_PressScaleIndDryCom = PressScaleIndDryCom
    a_RefPressDryCom      = RefPressDryCom

!!$    a_AbsCoefH2OVap = 1.0d-2
    a_AbsCoefH2OVap = AbsCoefH2OVap
    a_PressScaleIndH2OVap = PressScaleIndH2OVap
    a_RefPressH2OVap      = RefPressH2OVap

    ! Initialization of modules used in this module
    !

    ! 謨｣荵ｱ繧堤┌隕悶＠縺滓叛蟆�ｼ晞＃譁ｹ遞句ｼ�
    ! Radiative transfer equation without considering scattering
    !
    call RadRTENonScatInit

    ! 蜊ｰ蟄� ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )

!!$    call MessageNotify( 'M', module_name, 'DelTime:' )
!!$    call MessageNotify( 'M', module_name, '  DelTime  = %f [%c]', &
!!$      & d = (/ DelTimeValue /), c1 = trim( DelTimeUnit ) )

    call MessageNotify( 'M', module_name, 'nbmax         = %d', i = (/ nbmax /) )

!!$    call MessageNotify( 'M', module_name, 'WNBnds              = (/ %*r /)', &
!!$      & r = real( a_WNBnds(0:nbmax) ), n = (/ nbmax+1 /) )
!!$    call MessageNotify( 'M', module_name, 'AbsCoefDryCom       = (/ %*r /)', &
!!$      & r = real( a_AbsCoefDryCom(1:nbmax) ), n = (/ nbmax /) )
!!$    call MessageNotify( 'M', module_name, 'PressScaleIndDryCom = (/ %*r /)', &
!!$      & r = real( a_PressScaleIndDryCom(1:nbmax) ), n = (/ nbmax /) )
!!$    call MessageNotify( 'M', module_name, 'RefPressDryCom      = (/ %*r /)', &
!!$      & r = real( a_RefPressDryCom(1:nbmax) ), n = (/ nbmax /) )
!!$    call MessageNotify( 'M', module_name, 'AbsCoefH2OVap       = (/ %*r /)', &
!!$      & r = real( a_AbsCoefH2OVap(1:nbmax) ), n = (/ nbmax /) )
!!$    call MessageNotify( 'M', module_name, 'PressScaleIndH2OVap = (/ %*r /)', &
!!$      & r = real( a_PressScaleIndH2OVap(1:nbmax) ), n = (/ nbmax /) )
!!$    call MessageNotify( 'M', module_name, 'RefPressH2OVap      = (/ %*r /)', &
!!$      & r = real( a_RefPressH2OVap(1:nbmax) ), n = (/ nbmax /) )

    do n = 1, nbmax
      call MessageNotify( 'M', module_name, '  %d : %f %f %f %f %f %f %f %f',       &
        & i = (/ n /),                                                              &
        & d = (/ a_WNBnds(n-1)*1.0e-2, a_WNBnds(n)*1.0e-2,                          &
        &        a_AbsCoefDryCom(n), a_PressScaleIndDryCom(n), a_RefPressDryCom(n), &
        &        a_AbsCoefH2OVap(n), a_PressScaleIndH2OVap(n), a_RefPressH2OVap(n) /) )
    end do

    call MessageNotify( 'M', module_name, 'RadActDryComMMR = %f', d = (/ RadActDryComMMR /) )
    call MessageNotify( 'M', module_name, 'DiffFact        = %f', d = (/ DiffFact /) )
!
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    rad_simple_LW_inited = .true.

  end subroutine RadSimpleLWInit

  !--------------------------------------------------------------------------------------

end module rad_simple_LW