!= 散乱を無視した放射伝達方程式
!
!= Radiative transfer equation without considering scattering
!
! Authors::   Yoshiyuki O. TAKAHASHI
! Version::   $Id: rad_rte_nonscat.f90,v 1.6 2014/06/29 07:48:29 yot Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module rad_rte_nonscat
  !
  != 散乱を無視した放射伝達方程式
  !
  != Radiative transfer equation without considering scattering
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! 
  !
  !== Procedures List
  !
!!$  ! RadDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadFluxOutput     :: 放射フラックスの出力
!!$  ! ------------            :: ------------
!!$  ! RadDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadFluxOutput     :: Output radiation fluxes
  !
  !== NAMELIST
  !
  ! NAMELIST#rad_rte_nonscat_nml
  !

  ! モジュール引用 ; USE statements
  !

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

  ! 物理・数学定数設定
  ! Physical and mathematical constants settings
  !
  use constants0, only: &
    & PI, &
                            ! $ \pi $ .
                            ! 円周率.  Circular constant
    & StB
                            ! $ \sigma_{SB} $ .
                            ! ステファンボルツマン定数.
                            ! Stefan-Boltzmann constant

  ! 物理定数設定
  ! Physical constants settings
  !
  use constants, only: &
    & Grav, &
                            ! $ g $ [m s-2]. 
                            ! 重力加速度. 
                            ! Gravitational acceleration
    & CpDry
                            ! $ C_p $ [J kg-1 K-1]. 
                            ! 乾燥大気の定圧比熱. 
                            ! Specific heat of air at constant pressure


  ! 格子点設定
  ! 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

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: RadRTENonScat
!!$  public :: RadRTENonScatAnotherForm
!!$  public :: RadRTENonScatMonoWithDiffFact
  public :: RadRTENonScatMonoSemiAnal
  public :: RadRTENonScatWrapper
  public :: RadRTENonScatInit


  ! 公開変数
  ! Public variables
  !
  logical, save, public:: rad_rte_nonscat_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag


  ! 非公開変数
  ! Private variables
  !
  integer , save              :: NumGaussNodeZAInt
  real(DP), save              :: DiffFact
  real(DP), save, allocatable :: a_CosZA( : )
  real(DP), save, allocatable :: a_GaussWeight( : )


  character(*), parameter:: module_name = 'rad_rte_nonscat'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_rte_nonscat.f90,v 1.6 2014/06/29 07:48:29 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

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

  subroutine RadRTENonScat(                                    &
    & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT,  & ! (in)
    & xyrr_Trans,                                              & ! (in)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                          & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !

    ! OpenMP
    !
    !$ use omp_lib


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP), intent(in ) :: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP), intent(in ) :: xy_IntDPFDT1     (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
    real(DP), intent(in ) :: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature
    real(DP), intent(in ) :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    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)
                              ! 長波地表温度変化. 
                              ! Upward longwave flux derivative
    real(DP), intent(out) :: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Downward longwave flux derivative


    ! 作業変数
    ! Work variables
    !
    integer :: js
    integer :: je

    integer :: nthreads
    integer, allocatable :: a_js(:)
    integer, allocatable :: a_je(:)

    integer :: n


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_rte_nonscat_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    nthreads = 1
    !$ nthreads  = omp_get_max_threads()
!!$    !$ write( 6, * ) "Number of processors : ", omp_get_num_procs()
!!$    !$ write( 6, * ) "Number of threads    : ", nthreads

    allocate( a_js(0:nthreads-1) )
    allocate( a_je(0:nthreads-1) )

    do n = 0, nthreads-1

      if ( n == 0 ) then
        a_js(n) = 1
      else
        a_js(n) = a_je(n-1) + 1
      end if

      a_je(n) = a_js(n  ) + jmax / nthreads - 1
      if ( n + 1 <= mod( jmax, nthreads ) ) then
        a_je(n) = a_je(n) + 1
      end if

    end do

    !$OMP PARALLEL DEFAULT(PRIVATE) &
    !$OMP SHARED( &
    !$OMP         nthreads, a_js, a_je, &
    !$OMP         xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT,  &
    !$OMP         xyrr_Trans, &
    !$OMP         xyr_RadLUwFlux, xyr_RadLDwFlux,                          &
    !$OMP         xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   &
    !$OMP       )

    !$OMP DO

    do n = 0, nthreads-1

      js = a_js(n)
      je = a_je(n)

      if ( js > je ) cycle

      call RadRTENonScatCore(                                &
        & js, je,                                                  & ! (in)
        & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT,  & ! (in)
        & xyrr_Trans,                                              & ! (in)
        & xyr_RadLUwFlux, xyr_RadLDwFlux,                          & ! (out)
        & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   & ! (out)
        & )

    end do


    !$OMP END DO
    !$OMP END PARALLEL


    deallocate( a_js )
    deallocate( a_je )


  end subroutine RadRTENonScat

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

  subroutine RadRTENonScatCore(                                &
    & js, je,                                                  & ! (in)
    & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT,  & ! (in)
    & xyrr_Trans,                                              & ! (in)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                          & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !


    ! 宣言文 ; Declaration statements
    !

    integer , intent(in ) :: js
    integer , intent(in ) :: je
    real(DP), intent(in ) :: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP), intent(in ) :: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP), intent(in ) :: xy_IntDPFDT1     (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
    real(DP), intent(in ) :: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature
    real(DP), intent(in ) :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    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)
                              ! 長波地表温度変化. 
                              ! Upward longwave flux derivative
    real(DP), intent(out) :: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Downward longwave flux derivative


    ! 作業変数
    ! Work variables
    !
    integer:: j
    integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_rte_nonscat_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if



    ! 放射フラックス計算
    ! Calculate radiation flux
    !

    !   Initialization
    !
    do k = 0, kmax
      do j = js, je
        xyr_RadLDwFlux(:,j,k) = 0.0_DP
        xyr_RadLUwFlux(:,j,k) = 0.0_DP
      end do
    end do
    !
    !   Downward flux
    !
    do k = kmax, 0, -1

      do kk = kmax, k+1, -1
        do j = js, je
          xyr_RadLDwFlux(:,j,k) = xyr_RadLDwFlux(:,j,k)          &
            & + xyz_IntPF(:,j,kk)                                &
            & * ( xyrr_Trans(:,j,k,kk-1) - xyrr_Trans(:,j,k,kk) )
        end do
      end do

    end do
    !
    !   Upward flux
    !
    !     Set upward flux
    !
    do k = 0, kmax
      do j = js, je

        xyr_RadLUwFlux(:,j,k) = xy_SurfIntPF(:,j) * xyrr_Trans(:,j,k,0)

        do kk = 1, k
          xyr_RadLUwFlux(:,j,k) = xyr_RadLUwFlux(:,j,k)          &
            & - xyz_IntPF(:,j,kk)                                &
            & * ( xyrr_Trans(:,j,k,kk-1) - xyrr_Trans(:,j,k,kk) )
        end do

      end do
    end do


    ! 放射フラックスの変化率の計算
    ! Calculate rate of change of radiative flux
    !
    do k = 0, kmax
      do j = js, je
        xyra_DelRadLUwFlux(:,j,k,0) =                         &
          & xy_SurfIntDPFDT(:,j) * xyrr_Trans(:,j,k,0)
      end do
    end do
    k = 0
    do j = js, je
      xyra_DelRadLUwFlux(:,j,k,1) = 0.0_DP
    end do
    do k = 1, kmax
      do j = js, je
        xyra_DelRadLUwFlux(:,j,k,1) =                         &
          & - xy_IntDPFDT1(:,j)                               &
          &   * ( xyrr_Trans(:,j,k,0) - xyrr_Trans(:,j,k,1) )
      end do
    end do
    do k = 0, kmax
      do j = js, je
        xyra_DelRadLDwFlux(:,j,k,0) = 0.0_DP
      end do
    end do
    k = 0
    do j = js, je
      xyra_DelRadLDwFlux(:,j,k,1) =                         &
        & + xy_IntDPFDT1(:,j)                               &
        &   * ( xyrr_Trans(:,j,k,0) - xyrr_Trans(:,j,k,1) )
    end do
    do k = 1, kmax
      do j = js, je
        xyra_DelRadLDwFlux(:,j,k,1) = 0.0_DP
      end do
    end do


  end subroutine RadRTENonScatCore

  !--------------------------------------------------------------------------------------
  ! This is a test version, and now (2014/06/29), this is not used because of 
  ! not good performance. So, this will be deleted in future, probably. 
  ! This routine calculates radiative flux by solving a radiative transfer 
  ! equation without including scattering. 
  ! The difference between this routine and RadRTENonScat is that this 
  ! routine solves a radiative transfer equation of a form integrated by parts.

  subroutine RadRTENonScatAnotherForm(                         &
    & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT,  & ! (in)
    & xyrr_Trans,                                              & ! (in)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                          & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP), intent(in ) :: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP), intent(in ) :: xy_IntDPFDT1     (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
    real(DP), intent(in ) :: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature
    real(DP), intent(in ) :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    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)
                              ! 長波地表温度変化. 
                              ! Upward longwave flux derivative
    real(DP), intent(out) :: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Downward longwave flux derivative


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyr_IntPF(0:imax-1, 1:jmax, 0:kmax)

    integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_rte_nonscat_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if



    ! 
    ! Calculate integrated Planck function at layer interface
    !
    xyr_IntPF(:,:,0   ) = xyz_IntPF(:,:,1   )
    do k = 1, kmax-1
      xyr_IntPF(:,:,k) = ( xyz_IntPF(:,:,k) + xyz_IntPF(:,:,k+1) ) / 2.0_DP
    end do
    xyr_IntPF(:,:,kmax) = xyz_IntPF(:,:,kmax)


    ! 放射フラックス計算
    ! Calculate radiation flux
    !

    !   Initialization
    !
    xyr_RadLDwFlux = 0.0_DP
    xyr_RadLUwFlux = 0.0_DP
    !
    !   Downward flux
    !
    do k = kmax, 0, -1

      xyr_RadLDwFlux(:,:,k) =   &
        &   xyr_IntPF(:,:,k   ) &
        & - xyr_IntPF(:,:,kmax) * xyrr_Trans(:,:,k,kmax)
      !
      do kk = kmax, k+1, -1
        xyr_RadLDwFlux(:,:,k) = xyr_RadLDwFlux(:,:,k)                    &
          & - ( xyrr_Trans(:,:,k,kk-1) + xyrr_Trans(:,:,k,kk) ) / 2.0_DP &
          &   * ( xyr_IntPF(:,:,kk-1) - xyr_IntPF(:,:,kk) )
      end do

    end do
    !
    !   Upward flux
    !
    !     Set upward flux
    !
    do k = 0, kmax

      xyr_RadLUwFlux(:,:,k) =                                         &
        &   ( xy_SurfIntPF - xyr_IntPF(:,:,0) ) * xyrr_Trans(:,:,k,0) &
        & + xyr_IntPF(:,:,k)
      !
      do kk = 1, k
        xyr_RadLUwFlux(:,:,k) = xyr_RadLUwFlux(:,:,k)                    &
          & + ( xyrr_Trans(:,:,k,kk-1) + xyrr_Trans(:,:,k,kk) ) / 2.0_DP &
          &   * ( xyr_IntPF(:,:,kk-1) - xyr_IntPF(:,:,kk) )
      end do

    end do


    ! 放射フラックスの変化率の計算
    ! Calculate rate of change of radiative flux
    !
    do k = 0, kmax
      xyra_DelRadLUwFlux(:,:,k,0) =                         &
        & xy_SurfIntDPFDT * xyrr_Trans(:,:,k,0)
    end do
    k = 0
    xyra_DelRadLUwFlux(:,:,k,1) = 0.0_DP
    k = 1
    xyra_DelRadLUwFlux(:,:,k,1) =                         &
      &   xy_IntDPFDT1                                    &
      &   * ( - xyrr_Trans(:,:,k,0) + 1.0_DP / 2.0_DP     &
      &       + ( xyrr_Trans(:,:,k,0) + xyrr_Trans(:,:,k,1) ) / 2.0_DP / 2.0_DP )
    do k = 2, kmax
      xyra_DelRadLUwFlux(:,:,k,1) =                         &
        &   xy_IntDPFDT1                                    &
        &   * ( - xyrr_Trans(:,:,k,0)                       &
        &       + ( xyrr_Trans(:,:,k,0) + xyrr_Trans(:,:,k,1) ) / 2.0_DP / 2.0_DP &
        &       + ( xyrr_Trans(:,:,k,1) + xyrr_Trans(:,:,k,2) ) / 2.0_DP / 2.0_DP )
    end do
    do k = 0, kmax
      xyra_DelRadLDwFlux(:,:,k,0) = 0.0_DP
    end do
    if ( kmax <= 1 ) then
      k = 0
      xyra_DelRadLDwFlux(:,:,k,1) =                         &
        & + xy_IntDPFDT1                                    &
        &   * (   1.0_DP - xyrr_Trans(:,:,k,k+1) )
      k = 1
      xyra_DelRadLDwFlux(:,:,k,1) = 0.0_DP
    else
      k = 0
      xyra_DelRadLDwFlux(:,:,k,1) =                         &
        & + xy_IntDPFDT1                                    &
        &   * (   1.0_DP                                      &
        &       - ( xyrr_Trans(:,:,k,k  ) + xyrr_Trans(:,:,k,k+1) ) / 2.0_DP / 2.0_DP &
        &       - ( xyrr_Trans(:,:,k,k+1) + xyrr_Trans(:,:,k,k+2) ) / 2.0_DP / 2.0_DP )
      k = 1
      xyra_DelRadLDwFlux(:,:,k,1) =                         &
        & + xy_IntDPFDT1                                    &
        &   * (   1.0_DP / 2.0_DP                           &
        &       - ( xyrr_Trans(:,:,k,k  ) + xyrr_Trans(:,:,k,k+1) ) / 2.0_DP / 2.0_DP )
    end if
    do k = 2, kmax
      xyra_DelRadLDwFlux(:,:,k,1) = 0.0_DP
    end do


  end subroutine RadRTENonScatAnotherForm

  !--------------------------------------------------------------------------------------
  ! This calculation method may not be good. Numerical integration over zenith 
  ! angle should be performed.
  subroutine RadRTENonScatMonoWithDiffFact(                    &
    & DiffFact,                                                & ! (in)
    & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT,  & ! (in)
    & xyz_OptDep,                                              & ! (in)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                          & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: DiffFact
                              ! Diffusivity factor
    real(DP), intent(in ) :: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP), intent(in ) :: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP), intent(in ) :: xy_IntDPFDT1     (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
    real(DP), intent(in ) :: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature
    real(DP), intent(in ) :: xyz_OptDep   (0:imax-1, 1:jmax, 1:kmax)
                              ! 光学的厚さ. 
                              ! Optical depth
    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)
                              ! 長波地表温度変化. 
                              ! Upward longwave flux derivative
    real(DP), intent(out) :: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Downward longwave flux derivative


    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyz_TransEachLayer(0:imax-1, 1:jmax, 1:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    real(DP) :: xyr_Trans0(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transmission coefficient from surface to layer interfaces
    real(DP) :: xyr_Trans1(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transmission coefficient from top of the lowest layer to layer interfaces

    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_rte_nonscat_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if



    ! 放射フラックス計算
    ! Calculate radiation flux
    !

    xyz_TransEachLayer = exp( - DiffFact * xyz_OptDep )
    !
    xyr_Trans0(:,:,0) = 1.0_DP
    do k = 1, kmax
      xyr_Trans0(:,:,k) = xyr_Trans0(:,:,k-1) * xyz_TransEachLayer(:,:,k)
    end do
    xyr_Trans1(:,:,0) = xyz_TransEachLayer(:,:,1)
    xyr_Trans1(:,:,1) = 1.0_DP
    do k = 2, kmax
      xyr_Trans1(:,:,k) = xyr_Trans0(:,:,k-1) * xyz_TransEachLayer(:,:,k)
    end do

    !   Initialization
    !
    xyr_RadLDwFlux = 0.0_DP
    xyr_RadLUwFlux = 0.0_DP
    !
    !   Downward flux
    !
    k = kmax
    xyr_RadLDwFlux(:,:,k) = 0.0_DP
    do k = kmax-1, 0, -1
      xyr_RadLDwFlux(:,:,k) =                                     &
        &   xyr_RadLDwFlux(:,:,k+1) * xyz_TransEachLayer(:,:,k+1) &
        & + xyz_IntPF(:,:,k+1) * ( 1.0_DP - xyz_TransEachLayer(:,:,k+1) )
    end do
    !
    !   Upward flux
    !
    k = 0
    xyr_RadLUwFlux(:,:,k) = xy_SurfIntPF
    do k = 1, kmax
      xyr_RadLUwFlux(:,:,k) =                                   &
        &   xyr_RadLUwFlux(:,:,k-1) * xyz_TransEachLayer(:,:,k) &
        & - xyz_IntPF(:,:,k) * ( xyz_TransEachLayer(:,:,k) - 1.0_DP )
    end do


    ! 放射フラックスの変化率の計算
    ! Calculate rate of change of radiative flux
    !
    do k = 0, kmax
      xyra_DelRadLUwFlux(:,:,k,0) =                         &
        & xy_SurfIntDPFDT * xyr_Trans0(:,:,k)
    end do
    k = 0
    xyra_DelRadLUwFlux(:,:,k,1) = 0.0_DP
    do k = 1, kmax
      xyra_DelRadLUwFlux(:,:,k,1) =                         &
        & - xy_IntDPFDT1 * ( xyr_Trans0(:,:,k) - xyr_Trans1(:,:,k) )
    end do
    do k = 0, kmax
      xyra_DelRadLDwFlux(:,:,k,0) = 0.0_DP
    end do
    k = 0
    xyra_DelRadLDwFlux(:,:,k,1) =                         &
      & + xy_IntDPFDT1 * ( xyr_Trans0(:,:,k) - xyr_Trans1(:,:,k) )
    do k = 1, kmax
      xyra_DelRadLDwFlux(:,:,k,1) = 0.0_DP
    end do


  end subroutine RadRTENonScatMonoWithDiffFact

  !--------------------------------------------------------------------------------------
  ! ###This memo should be checked.###
  ! It should be noted that IntPF is Planck function integrated over 
  ! wavenumber and azimuthal angle and divided by 2. 
!!$  ! So, pi is multiplied after integrated over wavenumber. 
  !
  subroutine RadRTENonScatMonoSemiAnal(                        &
    & xy_SurfAlbedo,                                           & ! (in)
    & xyr_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT,  & ! (in)
    & xyz_OptDep,                                              & ! (in)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                          & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                   & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xy_SurfAlbedo    (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xyr_IntPF        (0:imax-1, 1:jmax, 0:kmax)
                              ! Integrated Planck function
    real(DP), intent(in ) :: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP), intent(in ) :: xy_IntDPFDT1     (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
    real(DP), intent(in ) :: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature
    real(DP), intent(in ) :: xyz_OptDep   (0:imax-1, 1:jmax, 1:kmax)
                              ! 光学的厚さ. 
                              ! Optical depth
    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)
                              ! 長波地表温度変化. 
                              ! Upward longwave flux derivative
    real(DP), intent(out) :: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Downward longwave flux derivative


    ! 作業変数
    ! Work variables
    !
    real(DP) :: CosZA
    real(DP) :: GaussWeight

    real(DP) :: xyz_TransEachLayer0(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_TransEachLayer (0:imax-1, 1:jmax, 1:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    real(DP) :: xyr_Trans0(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transmission coefficient from surface to layer interfaces
    real(DP) :: xyr_Trans1(0:imax-1, 1:jmax, 0:kmax)
                              ! 
                              ! Transmission coefficient from top of the lowest layer to layer interfaces

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

    real(DP) :: xyr_UwFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_DwFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyra_DelUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP) :: xyra_DelDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)

    real(DP) :: IntFact

    integer:: i
    integer:: j
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: n


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_rte_nonscat_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if



    ! 放射フラックス計算
    ! Calculate radiation flux
    !

    xyz_TransEachLayer0 = exp( - xyz_OptDep )
    !  This is ad hoc treatment to avoid underflow.
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_TransEachLayer0(i,j,k) < 1.0e-50_DP ) then
            xyz_TransEachLayer0(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do

    do k = 1, kmax
      xyz_DPFDOptDep(:,:,k) =                      &
        &   ( xyr_IntPF(:,:,k-1) - xyr_IntPF(:,:,k) ) &
        & / max( xyz_OptDep(:,:,k), 1.0d-100 )
    end do

    if ( NumGaussNodeZAInt > 0 ) then
      ! for case with Gaussian quadrature
      IntFact = 2.0_DP
    else
      ! for two-stream approximation or case with diffusivity factor
      IntFact = 1.0_DP
    end if


    !   Initialization
    !
    xyr_RadLUwFlux     = 0.0_DP
    xyr_RadLDwFlux     = 0.0_DP
    xyra_DelRadLUwFlux = 0.0_DP
    xyra_DelRadLDwFlux = 0.0_DP


    !   Loop for Gaussian quadrature
    !
    loop_gq : do n = 1, max( NumGaussNodeZAInt, 1 )

      ! Preparetion
      !
      if ( NumGaussNodeZAInt > 0 ) then
        CosZA = a_CosZA(n)
      else
        CosZA = 1.0_DP / DiffFact
      end if
      !
      xyz_TransEachLayer = ( xyz_TransEachLayer0 )**(1.0d0/CosZA)
      !
      xyr_Trans0(:,:,0) = 1.0_DP
      do k = 1, kmax
        xyr_Trans0(:,:,k) = xyr_Trans0(:,:,k-1) * xyz_TransEachLayer(:,:,k)
      end do
      xyr_Trans1(:,:,0) = xyz_TransEachLayer(:,:,1)
      xyr_Trans1(:,:,1) = 1.0_DP
      do k = 2, kmax
        xyr_Trans1(:,:,k) = xyr_Trans0(:,:,k-1) * xyz_TransEachLayer(:,:,k)
      end do

      !
      !   Downward flux
      !
      k = kmax
      xyr_DwFlux(:,:,k) = 0.0_DP
      do k = kmax-1, 0, -1
        xyr_DwFlux(:,:,k) =                                          &
          &   ( xyr_DwFlux(:,:,k+1) - IntFact * xyr_IntPF(:,:,k+1) ) &
          & * xyz_TransEachLayer(:,:,k+1)                            &
          & + IntFact * xyr_IntPF(:,:,k)                             &
          & - IntFact * CosZA * xyz_DPFDOptDep(:,:,k+1)              &
          &   * ( 1.0_DP - xyz_TransEachLayer(:,:,k+1) )
      end do
      !
      !   Upward flux
      !
      k = 0
      xyr_UwFlux(:,:,k) =                       &
        &   IntFact * xy_SurfIntPF              &
        & + xy_SurfAlbedo * xyr_DwFlux(:,:,0)
      do k = 1, kmax
        xyr_UwFlux(:,:,k) =                                          &
          &   ( xyr_UwFlux(:,:,k-1) - IntFact * xyr_IntPF(:,:,k-1) ) &
          & * xyz_TransEachLayer(:,:,k)                              &
          & + IntFact * xyr_IntPF(:,:,k)                             &
          & + IntFact * CosZA * xyz_DPFDOptDep(:,:,k)                &
          &   * ( 1.0_DP - xyz_TransEachLayer(:,:,k) )
      end do


      ! 放射フラックスの変化率の計算
      ! Calculate rate of change of radiative flux
      !
      do k = 0, kmax
        xyra_DelUwFlux(:,:,k,0) = IntFact * xy_SurfIntDPFDT * xyr_Trans0(:,:,k)
      end do
      do k = 0, kmax
        xyra_DelDwFlux(:,:,k,0) = 0.0_DP
      end do
      xyra_DelUwFlux(:,:,:,1) = 0.0_DP
      xyra_DelDwFlux(:,:,:,1) = 0.0_DP


      ! Sum over zenith angle
      !
      if ( NumGaussNodeZAInt > 0 ) then
        GaussWeight = a_GaussWeight( n )

        xyr_RadLUwFlux = xyr_RadLUwFlux &
          & + xyr_UwFlux * CosZA * GaussWeight
        xyr_RadLDwFlux = xyr_RadLDwFlux &
          & + xyr_DwFlux * CosZA * GaussWeight

        xyra_DelRadLUwFlux = xyra_DelRadLUwFlux &
          & + xyra_DelUwFlux * CosZA * GaussWeight
        xyra_DelRadLDwFlux = xyra_DelRadLDwFlux &
          & + xyra_DelDwFlux * CosZA * GaussWeight
      else
        xyr_RadLUwFlux = xyr_UwFlux
        xyr_RadLDwFlux = xyr_DwFlux

        xyra_DelRadLUwFlux = xyra_DelUwFlux
        xyra_DelRadLDwFlux = xyra_DelDwFlux
      end if

    end do loop_gq


  end subroutine RadRTENonScatMonoSemiAnal

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

  subroutine RadRTENonScatWrapper(                    &
    & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans, & ! (in )
    & xyr_RadLUwFlux, xyr_RadLDwFlux,                 & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux,         & ! (out)
    & WNs, WNe, NumGaussNode                          & ! (in ) optional
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !


    ! プランク関数の計算
    ! Calculate Planck function
    !
    use planck_func, only :                            &
      & Integ_PF_GQ_Array3D, Integ_PF_GQ_Array2D, Integ_DPFDT_GQ_Array2D

    ! 宣言文 ; Declaration statements
    !

    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 ) :: xy_SurfEmis (0:imax-1, 1:jmax)
                              ! 惑星表面射出率. 
                              ! Surface emissivity
    real(DP), intent(in ) :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    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)
                              ! 長波地表温度変化. 
                              ! Upward longwave flux derivative
    real(DP), intent(out) :: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Downward longwave flux derivative
    real(DP), intent(in ), optional :: WNs
    real(DP), intent(in ), optional :: WNe
    integer , intent(in ), optional :: NumGaussNode


    ! 作業変数
    ! Work variables
    !
    real(DP):: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP):: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP):: xy_IntDPFDT1     (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
    real(DP):: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature


    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_rte_nonscat_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Check arguments
    !
    if ( present( WNs ) .or. present( WNe ) .or. present( NumGaussNode ) ) then
      if ( .not. ( present( WNs ) .and. present( WNe ) .and. present( NumGaussNode ) ) ) then
        call MessageNotify( 'E', module_name, &
          & 'All of WNs, WNe, and NumGaussNode have to be present.' )
      end if
    end if


    if ( present( WNs ) ) then
      ! Case for non-grey atmosphere
      !

      ! Integrate Planck function and temperature derivative of it
      !
      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_SurfIntDPFDT                  & ! (out)
        & )

      xyz_IntPF       =               PI * xyz_IntPF
      xy_SurfIntPF    = xy_SurfEmis * PI * xy_SurfIntPF
      xy_IntDPFDT1    =               PI * xy_IntDPFDT1
      xy_SurfIntDPFDT = xy_SurfEmis * PI * xy_SurfIntDPFDT

    else

      ! Case for grey atmosphere
      !
      xyz_IntPF       =                        StB * xyz_Temp**4
      xy_SurfIntPF    = xy_SurfEmis          * StB * xy_SurfTemp**4
      xy_IntDPFDT1    =               4.0_DP * StB * xyz_Temp(:,:,1)**3
      xy_SurfIntDPFDT = xy_SurfEmis * 4.0_DP * StB * xy_SurfTemp**3

    end if


    call RadRTENonScat(                                         &
      & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT, & ! (in)
      & xyrr_Trans,                                             & ! (in)
      & xyr_RadLUwFlux, xyr_RadLDwFlux,                         & ! (out)
      & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux                  & ! (out)
      & )


  end subroutine RadRTENonScatWrapper

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

  subroutine RadRTENonScatInit
    !
    ! rad_rte_nonscat モジュールの初期化を行います. 
    ! NAMELIST#rad_rte_nonscat_nml の読み込みはこの手続きで行われます. 
    !
    ! "rad_rte_nonscat" module is initialized. 
    ! "NAMELIST#rad_rte_nonscat_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: toChar

    ! ガウス重み, 分点の計算
    ! Calculate Gauss node and Gaussian weight
    !
    use gauss_quad, only : GauLeg

    ! 宣言文 ; Declaration statements
    !

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read


    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /rad_rte_nonscat_nml/ &
      & DiffFact,                  &
      & NumGaussNodeZAInt
          !
          ! デフォルト値については初期化手続 "rad_rte_nonscat#RadRTENonScatInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_rte_nonscat#RadRTENonScatInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( rad_rte_nonscat_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    NumGaussNodeZAInt = 3
    DiffFact          = 1.66_DP


    ! 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_rte_nonscat_nml,      & ! (out)
        & iostat = iostat_nml )             ! (out)
      close( unit_nml )

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


    if ( NumGaussNodeZAInt > 0 ) then
      allocate( a_CosZA      ( NumGaussNodeZAInt ) )
      allocate( a_GaussWeight( NumGaussNodeZAInt ) )
      call GauLeg( 0.0_DP, 1.0_DP, NumGaussNodeZAInt, a_CosZA, a_GaussWeight )
    end if


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'NumGaussNodeZAInt = %d', i = (/ NumGaussNodeZAInt /) )
    call MessageNotify( 'M', module_name, 'DiffFact          = %f', d = (/ DiffFact /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    rad_rte_nonscat_inited = .true.

  end subroutine RadRTENonScatInit

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

end module rad_rte_nonscat