!= 
!
!= Solve radiative transfer equation in two stream approximation
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: rad_rte_two_stream_app.f90,v 1.7 2015/03/11 04:48:47 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module rad_rte_two_stream_app
  !
  != 
  !
  != Solve radiative transfer equation in two stream approximation
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! Solve radiative transfer equation in two stream approximation. 
  ! Analytic solution is used to calculate radiative flux in a homogeneous atmosphere 
  ! in which the single scattering albedo and the asymmetry factor are constant. 
  ! Radiative transfer equation is solved numerically with the method by Toon et al. 
  ! (1989) to calculate radiative flux in an inhomogeneous atmosphere. 
  !
  !
  !== References
  !
  !  Toon, O. B., C. P. McKay, and A. P. Ackerman, 
  !    Rapid calculation of radiative heating rates and photodissociation rates 
  !    in inhomogeneous multiple scattering atmospheres, 
  !    J. Geophys. Res., 94, 16287-16301, 1989.
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
!!$  ! NAMELIST#radiation_DennouAGCM_nml
  !

  !
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

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

  ! 物理・数学定数設定
  ! Physical and mathematical constants settings
  !
  use constants0, only: &
    & PI
                            ! $ \pi $ .
                            ! 円周率.  Circular constant

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

  implicit none

  private


  ! Local variables
  !
  real(DP), parameter :: DelTauThreshold = 1.0e-10_DP

!!$  integer, save :: NGaussQuad
  integer, parameter :: NGaussQuad = 8
  real(DP), save     :: a_GQP(1:NGaussQuad)
  real(DP), save     :: a_GQW(1:NGaussQuad)

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

  integer, parameter :: IDScatApproxEddington = 11
  integer, parameter :: IDScatApproxHemiMean  = 12

  public :: RadRTETwoStreamAppHomogAtm
!!$  public :: RadRTETwoStreamApp
  public :: RadRTETwoStreamAppSW
  public :: RadRTETwoStreamAppLW
  public :: RadRTETwoStreamAppInit


  ! INTERFACE 文 ; INTERFACE statements
  !
!!$  interface RadRTETwoStreamApp
!!$    module procedure RadRTETwoStreamAppWrapper, RadRTETwoStreamAppCore
!!$  end interface
!!$  interface RadRTETwoStreamApp
!!$    module procedure RadRTETwoStreamAppWrapper
!!$  end interface



  character(*), parameter:: module_name = 'rad_rte_two_stream_app'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_rte_two_stream_app.f90,v 1.7 2015/03/11 04:48:47 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

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

contains

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

  subroutine RadRTETwoStreamAppHomogAtm(                               &
    & xy_SurfAlbedo, SolarFluxAtTOA, xy_InAngle, SSA, AF, xyr_OptDep,  & ! (in )
    & xyr_RadSUwFlux, xyr_RadSDwFlux,                                  & ! (out)
    & FlagSemiInfAtm, FlagSL09                                         & ! (in ) optional
    & )

    ! Calculate radiative flux in a homogeneous scattering and absorbing atmosphere. 
    ! Analytical solution is used for calculation of radiative flux. 
    ! Radiative flux in a semi-infinite atmosphere is calculated if FlagSemiInfAtm 
    ! is .true.. If FlagSemiInfAtm is not given or is .false., radiative flux in a finite
    ! atmosphere (bounded by the surface) is calculated.
    !
    ! If FlagSL09 is .true., short wave radiative flux is calculated with the method by 
    ! Schneider and Liu (2009). 
    !
    ! See Meador and Weaver (19??), Toon et al. (1989), Liou (200?), and so on for 
    ! details of radiative transfer equation in this system.


    real(DP), intent(in ) :: xy_SurfAlbedo(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: SolarFluxAtTOA
    real(DP), intent(in ) :: xy_InAngle   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: SSA
    real(DP), intent(in ) :: AF
    real(DP), intent(in ) :: xyr_OptDep    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadSUwFlux(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadSDwFlux(0:imax-1, 1:jmax, 0:kmax)
    logical , intent(in ), optional :: FlagSemiInfAtm
    logical , intent(in ), optional :: FlagSL09

    !
    ! cosz     : cosine of solar zenith angle
    ! cosz2    : cosz squared
    !
    real(DP) :: xy_cosSZA     ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInv  ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )

    real(DP) :: SSAAdj
    real(DP) :: AFAdj
    real(DP) :: xyr_OptDepAdj(0:imax-1, 1:jmax, 0:kmax )

    real(DP) :: Lambda
    real(DP) :: LSigma
    real(DP) :: Gam1
    real(DP) :: Gam2
    real(DP) :: xy_Gam3   (0:imax-1, 1:jmax)
    real(DP) :: xy_Gam4   (0:imax-1, 1:jmax)
    real(DP) :: xyr_Trans (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_TMPVal(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_CUp   (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_CDo   (0:imax-1, 1:jmax, 0:kmax)

    real(DP) :: xy_k1           (0:imax-1, 1:jmax)
    real(DP) :: xy_k2           (0:imax-1, 1:jmax)
    real(DP) :: xyr_ExpLamOptDep(0:imax-1, 1:jmax, 0:kmax)

    real(DP) :: xy_DWRadSFluxAtTOA(0:imax-1, 1:jmax)

    logical  :: FlagSemiInfAtmLV
    logical  :: FlagSL09LV

    integer  :: i
    integer  :: j
    integer  :: k


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

    ! Check flags
    !
    FlagSL09LV = .false.
    if ( present( FlagSL09 ) ) then
      if ( FlagSL09 ) then
        FlagSL09LV = .true.
      end if
    end if
    FlagSemiInfAtmLV = .false.
    if ( present( FlagSemiInfAtm ) ) then
      if ( FlagSemiInfAtm ) then
        FlagSemiInfAtmLV = .true.
      end if
    end if
    if ( FlagSL09LV .and. ( .not. FlagSemiInfAtmLV ) ) then
      call MessageNotify( 'E', module_name, 'FlagSemiInfAtm has to be .true. when FlagSL09 is .true.' )
    end if


    if ( FlagSL09LV ) then
      SSAAdj        = SSA
      AFAdj         = AF
      xyr_OptDepAdj = xyr_OptDep
    else
      ! Delta-function adjustment
      !
      SSAAdj        =   ( 1.0_DP - AF**2 ) * SSA  / ( 1.0_DP - SSA * AF**2 )
      AFAdj         = AF / ( 1.0_DP + AF )
      xyr_OptDepAdj = ( 1.0_DP - SSA * AF**2 ) * xyr_OptDep
    end if


    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_InAngle(i,j) > 0.0_DP ) then
          xy_cosSZA     (i,j) = 1.0_DP / xy_InAngle(i,j)
          xy_cosSZAInv  (i,j) = xy_InAngle(i,j)
          xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
        else
          xy_cosSZA     (i,j) = 0.0_DP
          xy_cosSZAInv  (i,j) = 0.0_DP
          xy_cosSZAInvsq(i,j) = 0.0_DP
        end if
      end do
    end do


    if ( FlagSL09LV ) then
      ! Coefficients for Hemispheric mean approximation
      !
      Gam1    = 2.0_DP - SSAAdj * ( 1.0_DP + AFAdj )
      Gam2    = SSAAdj * ( 1.0_DP - AFAdj )
      xy_Gam3 = 1.0d100
      xy_Gam4 = 1.0d100
    else
      ! Coefficients for Eddington approximation
      !
      Gam1    =  ( 7.0_DP - SSAAdj * ( 4.0_DP + 3.0_DP * AFAdj ) ) / 4.0_DP
      Gam2    = -( 1.0_DP - SSAAdj * ( 4.0_DP - 3.0_DP * AFAdj ) ) / 4.0_DP
      xy_Gam3 =  ( 2.0_DP - 3.0_DP * AFAdj * xy_cosSZA )        / 4.0_DP
      xy_Gam4 = 1.0_DP - xy_Gam3
    end if


    Lambda = sqrt( Gam1**2 - Gam2**2 )
    LSigma = Gam2 / ( Gam1 + Lambda )

    do k = 0, kmax
      xyr_Trans(:,:,k) = exp( - xyr_OptDepAdj(:,:,k) * xy_cosSZAInv )
    end do


    if ( FlagSL09LV ) then
      xyr_CUp = 0.0_DP
      xyr_CDo = 0.0_DP
      xy_DWRadSFluxAtTOA = SolarFluxAtTOA * xy_CosSZA
    else
      do k = 0, kmax
        xyr_TMPVal(:,:,k) =                                       &
          &   SSAAdj * SolarFluxAtTOA * xyr_Trans(:,:,k)          &
          & / ( Lambda**2 - xy_cosSZAInvsq )
        xyr_CUp(:,:,k) = xyr_TMPVal(:,:,k)                        &
          &   * ( ( Gam1 - xy_cosSZAInv ) * xy_Gam3 + Gam2 * xy_Gam4 )
        xyr_CDo(:,:,k) = xyr_TMPVal(:,:,k)                        &
          &   * ( ( Gam1 + xy_cosSZAInv ) * xy_Gam4 + Gam2 * xy_Gam3 )
        xy_DWRadSFluxAtTOA = 0.0_DP
      end do
    end if


    ! A variable used in the following calculation
    !
    xyr_ExpLamOptDep = exp( Lambda * xyr_OptDepAdj )

    if ( FlagSemiInfAtmLV ) then
      xy_k1 = 0.0_DP
      xy_k2 = xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax)
    else
      xy_k2 = &
        & (                                                 &
        &     ( xy_SurfAlbedo * LSigma - 1.0_DP )           &
        &       * ( xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax) ) * xyr_ExpLamOptDep(:,:,0)  &
        &   + LSigma                                                                    &
        &       * ( - xyr_CUp(:,:,0)                                                    &
        &           + xy_SurfAlbedo * (   xyr_CDo(:,:,0)                                &
        &                               + SolarFluxAtTOA * xy_CosSZA * xyr_Trans(:,:,0) ) &
        &         )                                                                     &
        & ) &
        & / &
        & ( &
        &     ( xy_SurfAlbedo * LSigma - 1.0_DP ) * xyr_ExpLamOptDep(:,:,0) &
        &   - ( xy_SurfAlbedo - LSigma ) * LSigma / xyr_ExpLamOptDep(:,:,0) &
        & )

      xy_k1 = ( xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax) - xy_k2 ) / LSigma
    end if


    ! Calculate radiative flux
    !
!!$    do k = 0, kmax
!!$      xyr_RadSFlux(:,:,k) =                                                             &
!!$        &   ( 1.0_DP - LSigma )                                                         &
!!$        &     * ( xy_k1 * xyr_ExpLamOptDep(:,:,k) - xy_k2 / xyr_ExpLamOptDep(:,:,k) )   &
!!$        & + xyr_CUp(:,:,k) - xyr_CDo(:,:,k)
!!$    end do
    do k = 0, kmax
      xyr_RadSUwFlux(:,:,k) =                         &
        &            xy_k1 * xyr_ExpLamOptDep(:,:,k)  &
        & + LSigma * xy_k2 / xyr_ExpLamOptDep(:,:,k)  &
        & + xyr_CUp(:,:,k)
      xyr_RadSDwFlux(:,:,k) =                         &
        &   LSigma * xy_k1 * xyr_ExpLamOptDep(:,:,k)  &
        & +          xy_k2 / xyr_ExpLamOptDep(:,:,k)  &
        & + xyr_CDo(:,:,k)
    end do


    if ( .not. FlagSL09LV ) then
      !
      ! Add direct solar insolation
      !
!!$      do k = 0, kmax
!!$        xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) &
!!$          & - SolarFluxAtTOA * xyr_Trans(:,:,k) * xy_cosSZA
!!$      end do
      do k = 0, kmax
        xyr_RadSDwFlux(:,:,k) = xyr_RadSDwFlux(:,:,k) &
          & + SolarFluxAtTOA * xyr_Trans(:,:,k) * xy_cosSZA
      end do
    end if

    do k = 0, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if( xy_cosSZA(i,j) <= 0.0_DP ) then
!!$            xyr_RadSFlux(i,j,k) = 0.0_DP
            xyr_RadSUwFlux(i,j,k) = 0.0_DP
            xyr_RadSDwFlux(i,j,k) = 0.0_DP
          end if
        end do
      end do
    end do


  end subroutine RadRTETwoStreamAppHomogAtm

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

  subroutine RadRTETwoStreamAppSW(      &
    & xyz_SSA, xyz_AF,                  & ! (in)
    & xyr_OptDep,                       & ! (in)
    & xy_SurfAlbedo,                    & ! (in)
    & SolarFluxTOA, xy_InAngle,         & ! (in)
    & xyr_RadUwFlux, xyr_RadDwFlux      & ! (out)
    & )

    ! USE statements
    !

    real(DP), intent(in ) :: xyz_SSA       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyz_AF        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyr_OptDep   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
    real(DP), intent(in ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(out) :: xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax)

    ! Local variables
    !
    integer :: IDScatApprox
    logical :: FlagTOAFlux
    logical :: FlagEmis
    logical :: FlagSrcFuncTech

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


    IDScatApprox    = IDScatApproxEddington
    FlagTOAFlux     = .true.
    FlagEmis        = .false.
    FlagSrcFuncTech = .false.

    call RadRTETwoStreamAppWrapper(                     &
      & IDScatApprox, FlagTOAFlux, FlagEmis, FlagSrcFuncTech, & ! (in)
      & xyz_SSA, xyz_AF,                                & ! (in)
      & xyr_OptDep,                                     & ! (in)
      & xy_SurfAlbedo,                                  & ! (in)
      & xyr_RadUwFlux, xyr_RadDwFlux,                   & ! (out)
      & SolarFluxTOA = SolarFluxTOA, xy_InAngle = xy_InAngle        & ! (in) optional
      & )


  end subroutine RadRTETwoStreamAppSW

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

  subroutine RadRTETwoStreamAppLW(                     &
    & xyz_SSA, xyz_AF,                                 & ! (in)
    & xyr_OptDep,                                      & ! (in)
    & xy_SurfAlbedo,                                   & ! (in)
    & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
    & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
    & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
    & )

    ! USE statements
    !

    real(DP), intent(in ) :: xyz_SSA          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_AF           (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyr_OptDep       (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfAlbedo    (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xyr_PFInted      (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfPFInted   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_SurfDPFDTInted(0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyr_RadUwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadDwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP), intent(out) :: xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


    ! Local variables
    !
    integer :: IDScatApprox
    logical :: FlagTOAFlux
    logical :: FlagEmis
    logical :: FlagSrcFuncTech


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


    IDScatApprox    = IDScatApproxHemiMean
    FlagTOAFlux     = .false.
    FlagEmis        = .true.
    FlagSrcFuncTech = .true.

    call RadRTETwoStreamAppWrapper(                     &
      & IDScatApprox, FlagTOAFlux, FlagEmis, FlagSrcFuncTech,  & ! (in)
      & xyz_SSA, xyz_AF,                                & ! (in)
      & xyr_OptDep,                                     & ! (in)
      & xy_SurfAlbedo,                                  & ! (in)
      & xyr_RadUwFlux, xyr_RadDwFlux,                   & ! (out)
      & xyr_PFInted = xyr_PFInted, xy_SurfPFInted = xy_SurfPFInted, xy_SurfDPFDTInted = xy_SurfDPFDTInted,    & ! (in) optional
      & xyra_DelRadUwFlux = xyra_DelRadUwFlux, xyra_DelRadDwFlux = xyra_DelRadDwFlux     & ! (out) optional
      & )


  end subroutine RadRTETwoStreamAppLW

  !--------------------------------------------------------------------------------------
  ! NOTE:
  ! xyr_PFInted    = \pi B
  ! xy_SurfPFInted = \epsilon \pi B
  !--------------------------------------------------------------------------------------

  subroutine RadRTETwoStreamAppWrapper(                &
    & IDScatApprox, FlagTOAFlux, FlagEmis, FlagSrcFuncTech,  & ! (in)
    & xyz_SSA, xyz_AF,                                 & ! (in)
    & xyr_OptDep,                                      & ! (in)
    & xy_SurfAlbedo,                                   & ! (in)
    & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
    & SolarFluxTOA, xy_InAngle,                        & ! (in) optional
    & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in) optional
    & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out) optional
    & )

    ! USE statements
    !

    ! OpenMP
    !
    !$ use omp_lib


    integer , intent(in ) :: IDScatApprox
    logical , intent(in ) :: FlagTOAFlux
    logical , intent(in ) :: FlagEmis
    logical , intent(in ) :: FlagSrcFuncTech
    real(DP), intent(in ) :: xyz_SSA       ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyz_AF        ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(in ) :: xyr_OptDep   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
    real(DP), intent(out) :: xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax)

    real(DP), intent(in ), optional :: SolarFluxTOA
    real(DP), intent(in ), optional :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ), optional :: xyr_PFInted      (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ), optional :: xy_SurfPFInted   (0:imax-1, 1:jmax)
    real(DP), intent(in ), optional :: xy_SurfDPFDTInted(0:imax-1, 1:jmax)
    real(DP), intent(out), optional :: xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP), intent(out), optional :: xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)

    ! Local variables
    !
    integer :: js
    integer :: je

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

    integer :: n


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


    !
    ! Arguments are checked.
    !
    if ( FlagTOAFlux ) then
      if ( .not. present( SolarFluxTOA ) ) &
        & call MessageNotify( 'E', module_name, 'SolarFluxTOA has to be present.' )
      if ( .not. present( xy_InAngle ) ) &
        & call MessageNotify( 'E', module_name, 'xy_InAngle has to be present.' )
    else
      if ( present( SolarFluxTOA ) ) &
        & call MessageNotify( 'E', module_name, 'SolarFluxTOA need not be present.' )
      if ( present( xy_InAngle ) ) &
        & call MessageNotify( 'E', module_name, 'xy_InAngle need not be present.' )
    end if

    if ( FlagEmis ) then
      if ( .not. present( xyr_PFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xyr_PFInted has to be present.' )
      if ( .not. present( xy_SurfPFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfPFInted has to be present.' )
      if ( .not. present( xy_SurfDPFDTInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted has to be present.' )
      if ( .not. present( xyra_DelRadUwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux has to be present.' )
      if ( .not. present( xyra_DelRadDwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux has to be present.' )
    else
      if ( present( xyr_PFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xyr_PFInted need not be present.' )
      if ( present( xy_SurfPFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfPFInted need not be present.' )
      if ( present( xy_SurfDPFDTInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted need not be present.' )
      if ( present( xyra_DelRadUwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux need not be present.' )
      if ( present( xyra_DelRadDwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux need not be present.' )
    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(nthreads,a_js,a_je, &
!!$    !$OMP        xyz_SSA, xyz_AF, &
!!$    !$OMP        SolarFluxTOA, &
!!$    !$OMP        xy_SurfAlbedo, &
!!$    !$OMP        IDScatApprox, FlagTOAFlux, FlagEmis, FlagSrcFuncTech, &
!!$    !$OMP        xy_InAngle, &
!!$    !$OMP        xyr_OptDep, &
!!$    !$OMP        xyr_RadUwFlux, &
!!$    !$OMP        xyr_RadDwFlux, &
!!$    !$OMP        xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted, &
!!$    !$OMP        xyra_DelRadUwFlux, xyra_DelRadDwFlux)
!!$
!!$    !$OMP DO

    do n = 0, nthreads-1

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

      if ( js > je ) cycle

!!$      write( 6, * ) n, js, je

      call RadRTETwoStreamAppCore(              &
        & IDScatApprox, FlagTOAFlux, FlagEmis, FlagSrcFuncTech, & ! (in)
        & xyz_SSA, xyz_AF,                      & ! (in)
        & xyr_OptDep,                           & ! (in)
        & xy_SurfAlbedo,                        & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux,         & ! (out)
        & js, je,                               & ! (in)
        & SolarFluxTOA = SolarFluxTOA, xy_InAngle = xy_InAngle,      & ! (in) optional
        & xyr_PFInted = xyr_PFInted, xy_SurfPFInted = xy_SurfPFInted, xy_SurfDPFDTInted = xy_SurfDPFDTInted,    & ! (in) optional
        & xyra_DelRadUwFlux = xyra_DelRadUwFlux, xyra_DelRadDwFlux = xyra_DelRadDwFlux     & ! (out) optional
        & )

    end do


!!$    !$OMP END DO
!!$    !$OMP END PARALLEL


    deallocate( a_js )
    deallocate( a_je )



  end subroutine RadRTETwoStreamAppWrapper

  !--------------------------------------------------------------------------------------
  ! NOTE:
  ! xyr_PFInted    = \pi B
  ! xy_SurfPFInted = \epsilon \pi B
  !--------------------------------------------------------------------------------------

  subroutine RadRTETwoStreamAppCore(                         &
    & IDScatApprox, FlagTOAFlux, FlagEmis, FlagSrcFuncTech,  & ! (in)
    & xyz_SSA, xyz_AF,                                       & ! (in)
    & xyr_OptDep,                                            & ! (in)
    & xy_SurfAlbedo,                                         & ! (in)
    & xyr_RadUwFlux, xyr_RadDwFlux,                          & ! (out)
    & js, je,                                                & ! (in)
    & SolarFluxTOA, xy_InAngle,                              & ! (in) optional
    & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,        & ! (in) optional
    & xyra_DelRadUwFlux, xyra_DelRadDwFlux                   & ! (out) optional
    & )

    ! USE statements
    !

!!$  use pf_module , only : pfint_gq_array


    integer , intent(in ) :: IDScatApprox
    logical , intent(in ) :: FlagTOAFlux
    logical , intent(in ) :: FlagEmis
    logical , intent(in ) :: FlagSrcFuncTech
    real(DP), intent(in ) :: xyz_SSA      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_AF       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyr_OptDep   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfAlbedo(0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax)

    integer , intent(in ) :: js
    integer , intent(in ) :: je

    real(DP), intent(in ), optional :: SolarFluxTOA
    real(DP), intent(in ), optional :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ), optional :: xyr_PFInted      (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ), optional :: xy_SurfPFInted   (0:imax-1, 1:jmax)
    real(DP), intent(in ), optional :: xy_SurfDPFDTInted(0:imax-1, 1:jmax)
    real(DP), intent(out), optional :: xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP), intent(out), optional :: xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


!!$    real(DP), intent(in ) :: gt         ( 0:imax-1, 1:jmax, 1:kmax )
!!$    real(DP), intent(in ) :: gts        ( 0:imax-1, 1:jmax )
!!$    real(DP), intent(in ) :: gph        ( 0:imax-1, 1:jmax, 0:kmax )

!!$    real(DP), intent(in ) :: emis  ( 0:imax-1, 1:jmax )
!!$    real(DP), intent(in ) :: wn1, wn2
!!$    integer , intent(in ) :: divnum


!!$  real(DP)    , intent(out) :: &
!!$    gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$    gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )

    real(DP) :: xy_SolarFluxTOA(0:imax-1, 1:jmax)

    !
    ! SSAAdj    : Delta-Function Adjusted Single-Scattering Albedo
    ! AFAdj     : Delta-Function Adjusted Asymmetry Factor
    ! OpDepAdj  : Delta-Function Adjusted Optical Depth
    !
    real(DP) :: xyz_SSAAdj     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_AFAdj      ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyr_OpDepAdj   ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: xyr_TransDirAdj( 0:imax-1, 1:jmax, 0:kmax )

    !
    ! gam?    : Coefficients of Generalized Radiative Transfer Equation
    !
    real(DP) :: xyz_Gam1( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_Gam2( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_Gam3( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_Gam4( 0:imax-1, 1:jmax, 1:kmax )

    !
    ! cosz     : cosine of solar zenith angle
    ! cosz2    : cosz squared
    !
    real(DP) :: xy_cosSZA     ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInv  ( 0:imax-1, 1:jmax )
    real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )

    !
    ! temporary variables
    !
    real(DP) :: xyz_DelTau( 0:imax-1, 1:jmax, 1:kmax )

    real(DP) :: xyz_Lambda ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_LGamma ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyaz_smalle( 0:imax-1, 1:jmax, 4, 1:kmax )

    real(DP) :: xy_SurfSrc( 0:imax-1, 1:jmax )

    !
    ! CUpB      : upward C at bottom of layer
    ! CUpT      : upward C at top of layer
    ! CDoB      : downward C at bottom of layer
    ! CDoT      : downward C at top of layer
    !
    real(DP) :: xyz_CUpB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CUpT( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CDoB( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CDoT( 0:imax-1, 1:jmax, 1:kmax )
    !
    real(DP) :: xyz_CUpBDir( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CUpTDir( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CDoBDir( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CDoTDir( 0:imax-1, 1:jmax, 1:kmax )
    !
    real(DP) :: xyz_CUpBEmi( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CUpTEmi( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CDoBEmi( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_CDoTEmi( 0:imax-1, 1:jmax, 1:kmax )

    real(DP) :: aa_TridiagMtx1( 1:imax*jmax, 1:kmax*2 )
    real(DP) :: aa_TridiagMtx2( 1:imax*jmax, 1:kmax*2 )
    real(DP) :: aa_TridiagMtx3( 1:imax*jmax, 1:kmax*2 )
    real(DP) :: aa_Vec        ( 1:imax*jmax, 1:kmax*2 )

    real(DP) :: xy_TMPVal( 0:imax-1, 1:jmax )

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

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

!!$  real(DP) :: gth( im, jm, km+1 ), pfinth( im, jm, km+1 )
!!$  real(DP) :: b0( ijs:ije, 1, km ), b1( ijs:ije, 1, km )
!!$  real(DP) :: gemis


    real(DP) :: xyr_IUw( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: xyr_IDw( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: FactG
    real(DP) :: FactH
    real(DP) :: FactJ
    real(DP) :: FactK
    real(DP) :: Alp1
    real(DP) :: Alp2
    real(DP) :: Sig1
    real(DP) :: Sig2

    integer  :: i, j, k, l
    integer  :: ms, me


    ! Variables for debug
    !
!!$    real(DP) :: xyr_UwFluxDebug( 0:imax-1, 1:jmax, 0:kmax )
!!$    real(DP) :: xyr_DwFluxDebug( 0:imax-1, 1:jmax, 0:kmax )


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


    !
    ! Arguments are checked.
    !
    if ( FlagTOAFlux ) then
      if ( .not. present( SolarFluxTOA ) ) &
        & call MessageNotify( 'E', module_name, 'SolarFluxTOA has to be present.' )
      if ( .not. present( xy_InAngle ) ) &
        & call MessageNotify( 'E', module_name, 'xy_InAngle has to be present.' )
    else
      if ( present( SolarFluxTOA ) ) &
        & call MessageNotify( 'E', module_name, 'SolarFluxTOA need not be present.' )
      if ( present( xy_InAngle ) ) &
        & call MessageNotify( 'E', module_name, 'xy_InAngle need not be present.' )
    end if

    if ( FlagEmis ) then
      if ( .not. present( xyr_PFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xyr_PFInted has to be present.' )
      if ( .not. present( xy_SurfPFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfPFInted has to be present.' )
      if ( .not. present( xy_SurfDPFDTInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted has to be present.' )
      if ( .not. present( xyra_DelRadUwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux has to be present.' )
      if ( .not. present( xyra_DelRadDwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux has to be present.' )
    else
      if ( present( xyr_PFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xyr_PFInted need not be present.' )
      if ( present( xy_SurfPFInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfPFInted need not be present.' )
      if ( present( xy_SurfDPFDTInted ) ) &
        & call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted need not be present.' )
      if ( present( xyra_DelRadUwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux need not be present.' )
      if ( present( xyra_DelRadDwFlux ) ) &
        & call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux need not be present.' )
    end if


    do k = 1, kmax
      do j = js, je
        do i = 0, imax-1
          if ( xyz_SSA(i,j,k) >= 1.0d0 ) then
            call MessageNotify( 'E', module_name, 'Single Scattering Albedo = %f', &
              & d = (/ xyz_SSA(i,j,k) /) )
          end if
        end do
      end do
    end do


    if ( FlagTOAFlux ) then

      do j = js, je
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_cosSZA     (i,j) = 1.0d0 / xy_InAngle(i,j)
            xy_cosSZAInv  (i,j) = xy_InAngle(i,j)
            xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
          else
            xy_cosSZA     (i,j) = 0.0d0
            xy_cosSZAInv  (i,j) = 0.0d0
            xy_cosSZAInvsq(i,j) = 0.0d0
          end if
        end do
      end do

      do j = js, je
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_SolarFluxTOA(i,j) = SolarFluxTOA
          else
            xy_SolarFluxTOA(i,j) = 0.0_DP
          end if
        end do
      end do

    else

      do j = js, je
        xy_cosSZA     (:,j) = 1.0d100
        xy_cosSZAInv  (:,j) = 1.0d100
        xy_cosSZAInvsq(:,j) = 1.0d100
      end do

      do j = js, je
        do i = 0, imax-1
          xy_SolarFluxTOA(i,j) = 1.0d100
        end do
      end do

    end if


    !
    ! Delta-Function Adjustment
    !
    do k = 1, kmax
      do j = js, je
        xyz_AFAdj (:,j,k) = xyz_AF(:,j,k) / ( 1.0d0 + xyz_AF(:,j,k) )
        xyz_SSAAdj(:,j,k) =   ( 1.0d0 - xyz_AF(:,j,k)**2 ) * xyz_SSA(:,j,k) &
          &               / ( 1.0d0 - xyz_SSA(:,j,k) * xyz_AF(:,j,k)**2 )
      end do
    end do
    !
    do k = 0, kmax
      do j = js, je
        xyr_OpDepAdj(:,j,k) = 0.0d0
      end do
    end do
    do k = kmax-1, 0, -1
      do j = js, je
        xyr_OpDepAdj(:,j,k) =                                       &
          &   xyr_OpDepAdj(:,j,k+1)                                 &
          & + ( 1.0d0 - xyz_SSA(:,j,k+1) * xyz_AF(:,j,k+1)**2 )     &
          &   * ( xyr_OptDep(:,j,k) - xyr_OptDep(:,j,k+1) )
      end do
    end do


    if ( FlagTOAFlux ) then
      do k = 0, kmax
        do j = js, je
          xyr_TransDirAdj(:,j,k) = exp( -xyr_OpDepAdj(:,j,k) * xy_cosSZAInv(:,j) )
        end do
      end do
    else
      do k = 0, kmax
        do j = js, je
          xyr_TransDirAdj(:,j,k) = 1.0d100
        end do
      end do
    end if


    select case ( IDScatApprox )
    case ( IDScatApproxEddington )

      !
      ! Eddington approximation
      !
      do k = 1, kmax
        do j = js, je
          xyz_Gam1(:,j,k) =  ( 7.0d0 - xyz_SSAAdj(:,j,k) * ( 4.0d0 + 3.0d0 * xyz_AFAdj(:,j,k) ) ) / 4.0d0
          xyz_Gam2(:,j,k) = -( 1.0d0 - xyz_SSAAdj(:,j,k) * ( 4.0d0 - 3.0d0 * xyz_AFAdj(:,j,k) ) ) / 4.0d0
          xyz_Gam3(:,j,k) =  ( 2.0d0 - 3.0d0 * xyz_AFAdj(:,j,k) * xy_cosSZA(:,j) )              / 4.0d0
          xyz_Gam4(:,j,k) = 1.0d0 - xyz_Gam3(:,j,k)
        end do
      end do

      do k = 1, kmax
        do j = js, je
          do i = 0, imax-1
!!$            xyz_Mu1(i,j,k) = ( 1.0_DP - xyz_SSAAdj(i,j,k) ) / ( xyz_Gam1(i,j,k) - xyz_Gam2(i,j,k) )
            xyz_Mu1(i,j,k) = 0.5_DP
          end do
        end do
      end do

    case ( IDScatApproxHemiMean )

      !
      ! Treatment if delta-function adjustment is not performed.
      !
!!$      do k = 1, kmax
!!$        do j = js, je
!!$          xyz_AFAdj (:,j,k) = xyz_AF (:,j,k)
!!$          xyz_SSAAdj(:,j,k) = xyz_SSA(:,j,k)
!!$        end do
!!$      end do
!!$      do k = 0, kmax
!!$        do j = js, je
!!$          xyr_OpDepAdj(:,j,k) = xyr_OptDep(:,j,k)
!!$        end do
!!$      end do

      !
      ! Hemispheric mean approximation
      !
      do k = 1, kmax
        do j = js, je
          xyz_Gam1(:,j,k) = 2.0_DP - xyz_SSAAdj(:,j,k) * ( 1.0_DP + xyz_AFAdj(:,j,k) )
          xyz_Gam2(:,j,k) = xyz_SSAAdj(:,j,k) * ( 1.0_DP - xyz_AFAdj(:,j,k) )
!!$          xyz_Gam3(:,j,k) = 1.0d100
!!$          xyz_Gam4(:,j,k) = 1.0d100
          xyz_Gam3(:,j,k) =  ( 2.0d0 - 3.0d0 * xyz_AFAdj(:,j,k) * xy_cosSZA(:,j) )              / 4.0d0
          xyz_Gam4(:,j,k) = 1.0d0 - xyz_Gam3(:,j,k)
        end do
      end do

      do k = 1, kmax
        do j = js, je
          do i = 0, imax-1
!!$            xyz_Mu1(i,j,k) = ( 1.0_DP - xyz_SSAAdj(i,j,k) ) / ( xyz_Gam1(i,j,k) - xyz_Gam2(i,j,k) )
            xyz_Mu1(i,j,k) = 0.5_DP
          end do
        end do
      end do

    case default
      call MessageNotify( 'E', module_name, 'Unexpected IDScatApprox, %d', &
        & i = (/ IDScatApprox /) )
    end select


    do k = 1, kmax
      do j = js, je
        xyz_DelTau(:,j,k) = xyr_OpDepAdj(:,j,k-1) - xyr_OpDepAdj(:,j,k)
      end do
    end do


    if ( FlagEmis ) then
      !
      ! Avoiding singularity when dtau equal to zero 
      !
      do k = 1, kmax
        do j = js, je
          do i = 0, imax-1
            if( xyz_DelTau(i,j,k) < DelTauThreshold ) then
              xyz_DelTau(i,j,k) = 0.0d0
            end if
          end do
        end do
      end do
    end if


    do k = 1, kmax
      do j = js, je
        ! In very small parameter space of single scattering albedo and asymmetry 
        ! factor close to asymmetry factor of 1, xyz_Gam1**2 - xyz_Gam2**2 becomes 
        ! negative value. (yot, 2011/11/20)
!!$        xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) )
        xyz_Lambda(:,j,k) = sqrt( max( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k), 0.0_DP ) )
!!$        xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) + 1.0d-10 )

        xyz_LGamma(:,j,k) = xyz_Gam2(:,j,k) / ( xyz_Gam1(:,j,k) + xyz_Lambda(:,j,k) )
      end do
    end do

    do k = 1, kmax
      do j = js, je
        xy_TMPVal(:,j)       = exp( - xyz_Lambda(:,j,k) * xyz_DelTau(:,j,k) )
        xyaz_smalle(:,j,1,k) = xyz_LGamma(:,j,k) * xy_TMPVal(:,j) + 1.0_DP
        xyaz_smalle(:,j,2,k) = xyz_LGamma(:,j,k) * xy_TMPVal(:,j) - 1.0_DP
        xyaz_smalle(:,j,3,k) = xy_TMPVal(:,j) + xyz_LGamma(:,j,k)
        xyaz_smalle(:,j,4,k) = xy_TMPVal(:,j) - xyz_LGamma(:,j,k)
      end do
    end do


    if ( FlagTOAFlux ) then
      do k = 1, kmax
        do j = js, je
          xy_TMPVal(:,j) =                                                        &
            &   xyz_SSAAdj(:,j,k) * xy_SolarFluxTOA(:,j)                          &
            &   * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k)     &
            &       + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) )                         &
            &   / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
          xyz_CUpBDir(:,j,k) = xy_TMPVal(:,j) * xyr_TransDirAdj(:,j,k-1)
          xyz_CUpTDir(:,j,k) = xy_TMPVal(:,j) * xyr_TransDirAdj(:,j,k  )
          !
          xy_TMPVal(:,j) =                                                        &
            &   xyz_SSAAdj(:,j,k) * xy_SolarFluxTOA(:,j)                          &
            &   * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k)     &
            &       + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) )                         &
            &   / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
          xyz_CDoBDir(:,j,k) = xy_TMPVal(:,j) * xyr_TransDirAdj(:,j,k-1)
          xyz_CDoTDir(:,j,k) = xy_TMPVal(:,j) * xyr_TransDirAdj(:,j,k  )
        end do
      end do
    else
      do k = 1, kmax
        do j = js, je
          xyz_CUpBDir(:,j,k) = 0.0_DP
          xyz_CUpTDir(:,j,k) = 0.0_DP
          xyz_CDoBDir(:,j,k) = 0.0_DP
          xyz_CDoTDir(:,j,k) = 0.0_DP
        end do
      end do
    end if


    if ( FlagEmis ) then

      do k = 1, kmax
        do j = js, je
          do i = 0, imax-1
            !
            ! Notice!
            ! Avoiding singularity when dtau equal to zero. 
            ! dtau occationally has much smaller value. 
            ! When this occurs, b1 cannot be calculated correctly. 
            !
            if( xyz_DelTau(i,j,k) /= 0.0_DP ) then
              xyz_B0(i,j,k) = xyr_PFInted(i,j,k)
              xyz_B1(i,j,k) = ( xyr_PFInted(i,j,k-1) - xyr_PFInted(i,j,k) ) / xyz_DelTau(i,j,k)
            else
!!$              xyz_B0(i,j,k) = 0.0_DP ! replace with a line below 2015/02/22 based on discussion with Onishi-san
              xyz_B0(i,j,k) = xyr_PFInted(i,j,k)
              xyz_B1(i,j,k) = 0.0_DP
            end if
          end do
        end do
      end do

      do k = 1, kmax
        do j = js, je
          do i = 0, imax-1
!!$            xyz_CUpB(i,j,k) = 2.0_DP * PI * Mu1 &
            xyz_CUpBEmi(i,j,k) = 2.0_DP * xyz_Mu1(i,j,k) &
              & * (   xyz_B0(i,j,k) &
              &     + xyz_B1(i,j,k) &
              &       * ( xyz_DelTau(i,j,k) + 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
!!$            xyz_CUpT(i,j,k) = 2.0_DP * PI * Mu1 &
            xyz_CUpTEmi(i,j,k) = 2.0_DP * xyz_Mu1(i,j,k) &
              & * (   xyz_B0(i,j,k) &
              &     + xyz_B1(i,j,k) &
              &       * ( 0.0_DP            + 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
!!$            xyz_CDoB(i,j,k) = 2.0_DP * PI * Mu1 &
            xyz_CDoBEmi(i,j,k) = 2.0_DP * xyz_Mu1(i,j,k) &
              & * (   xyz_B0(i,j,k) &
              &     + xyz_B1(i,j,k) &
              &       * ( xyz_DelTau(i,j,k) - 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
!!$            xyz_CDoT(i,j,k) = 2.0_DP * PI * Mu1 &
            xyz_CDoTEmi(i,j,k) = 2.0_DP * xyz_Mu1(i,j,k) &
              & * (   xyz_B0(i,j,k) &
              &     + xyz_B1(i,j,k) &
              &       * ( 0.0_DP            - 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
          end do
        end do
      end do

    else

      do k = 1, kmax
        do j = js, je
          xyz_CUpBEmi(:,j,k) = 0.0_DP
          xyz_CUpTEmi(:,j,k) = 0.0_DP
          xyz_CDoBEmi(:,j,k) = 0.0_DP
          xyz_CDoTEmi(:,j,k) = 0.0_DP
        end do
      end do

    end if

    do k = 1, kmax
      do j = js, je
        xyz_CUpB(:,j,k) = xyz_CUpBDir(:,j,k) + xyz_CUpBEmi(:,j,k)
        xyz_CUpT(:,j,k) = xyz_CUpTDir(:,j,k) + xyz_CUpTEmi(:,j,k)
        xyz_CDoB(:,j,k) = xyz_CDoBDir(:,j,k) + xyz_CDoBEmi(:,j,k)
        xyz_CDoT(:,j,k) = xyz_CDoTDir(:,j,k) + xyz_CDoTEmi(:,j,k)
      end do
    end do


    if ( FlagEmis ) then
      do j = js, je
        do i = 0, imax-1
          xy_SurfSrc(:,j) = xy_SurfPFInted(:,j)
        end do
      end do
!!$      else if( sw .eq. 2 ) then
!!$         do ij = ijs, ije
!!$!            gemis = 1.0d0
!!$            gemis = emis( ij, 1 )
!!$            ea( (ij-ijs+1), l ) &
!!$                 = -CUpB( ij, 1, km ) + ( 1.0d0 - gemis ) * CDoB( ij, 1, km ) &
!!$                 + gemis * pi * pfinth( ij, 1, km+1 )
!!$         end do
    else
      do j = js, je
        do i = 0, imax-1
          xy_SurfSrc(:,j) = 0.0_DP
        end do
      end do
    end if

    if ( FlagTOAFlux ) then
      do j = js, je
        xy_SurfSrc(:,j) = xy_SurfSrc(:,j)                            &
          & + xy_SurfAlbedo(:,j) * xy_SolarFluxTOA(:,j) * xyr_TransDirAdj(:,j,0) * xy_cosSZA(:,j)
      end do
    end if



    k = 1
    l = 1
    do j = js, je
      do i = 0, imax-1
        aa_TridiagMtx1( i+imax*(j-1)+1, l ) =                                     &
          & 0.0_DP
        aa_TridiagMtx2( i+imax*(j-1)+1, l ) =                                     &
          &     xyaz_smalle(i,j,1,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,3,k)
        aa_TridiagMtx3( i+imax*(j-1)+1, l ) =                                     &
          & - ( xyaz_smalle(i,j,2,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,4,k) )
      end do
    end do
    do j = js, je
      do i = 0, imax-1
        aa_Vec( i+imax*(j-1)+1, l ) =                                              &
          & - xyz_CUpB(i,j,k)                                                      &
          & + xy_SurfAlbedo(i,j) * xyz_CDoB(i,j,k)                                 &
          & + xy_SurfSrc(i,j)
      end do
    end do


    do k = 1, kmax-1

      do j = js, je
        do i = 0, imax-1

          l = 2 * k     ! equation number
          !
          aa_TridiagMtx1( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,1,k  ) * xyaz_smalle(i,j,2,k+1)                       &
            & - xyaz_smalle(i,j,3,k  ) * xyaz_smalle(i,j,4,k+1)
          aa_TridiagMtx2( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,2,k  ) * xyaz_smalle(i,j,2,k+1)                       &
            & - xyaz_smalle(i,j,4,k  ) * xyaz_smalle(i,j,4,k+1)
          aa_TridiagMtx3( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,1,k+1) * xyaz_smalle(i,j,4,k+1)                       &
            & - xyaz_smalle(i,j,2,k+1) * xyaz_smalle(i,j,3,k+1)
          aa_Vec        ( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,2,k+1) * ( - xyz_CDoT(i,j,k  ) + xyz_CDoB(i,j,k+1) )  &
            & - xyaz_smalle(i,j,4,k+1) * ( - xyz_CUpT(i,j,k  ) + xyz_CUpB(i,j,k+1) )

          l = 2 * k + 1 ! equation number
          !
          aa_TridiagMtx1( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,2,k  ) * xyaz_smalle(i,j,3,k  )                       &
            & - xyaz_smalle(i,j,1,k  ) * xyaz_smalle(i,j,4,k  )
          aa_TridiagMtx2( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,1,k  ) * xyaz_smalle(i,j,1,k+1)                       &
            & - xyaz_smalle(i,j,3,k  ) * xyaz_smalle(i,j,3,k+1)
          aa_TridiagMtx3( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,3,k  ) * xyaz_smalle(i,j,4,k+1)                       &
            & - xyaz_smalle(i,j,1,k  ) * xyaz_smalle(i,j,2,k+1)
          aa_Vec        ( i+imax*(j-1)+1, l ) =                                       &
            &   xyaz_smalle(i,j,3,k  ) * ( - xyz_CDoT(i,j,k  ) + xyz_CDoB(i,j,k+1) )  &
            & - xyaz_smalle(i,j,1,k  ) * ( - xyz_CUpT(i,j,k  ) + xyz_CUpB(i,j,k+1) )
        end do
      end do
    end do


    k = kmax
    l = 2 * kmax
    do j = js, je
      do i = 0, imax-1
        aa_TridiagMtx1( i+imax*(j-1)+1, l ) =  xyaz_smalle(i,j,1,k)
        aa_TridiagMtx2( i+imax*(j-1)+1, l ) =  xyaz_smalle(i,j,2,k)
        aa_TridiagMtx3( i+imax*(j-1)+1, l ) = 0.0d0
        aa_Vec        ( i+imax*(j-1)+1, l ) = -xyz_CDoT(i,j,k) + 0.0d0
      end do
    end do

    ms = 0      + imax*(js-1)+1
    me = imax-1 + imax*(je-1)+1
    call tridiag( imax*jmax, 2*kmax, aa_TridiagMtx1, aa_TridiagMtx2, aa_TridiagMtx3, aa_Vec, ms, me )

    if ( .not. FlagSrcFuncTech ) then

      k = 1
      do j = js, je
        do i = 0, imax-1
          xyr_RadUwFlux(i,j,k-1) =                                    &
            &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
            & - aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,2,k) &
            & + xyz_CUpB(i,j,k)
          xyr_RadDwFlux(i,j,k-1) =                                    &
            &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
            & - aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,4,k) &
            & + xyz_CDoB(i,j,k)
        end do
      end do

      do k = 1, kmax
        do j = js, je
          do i = 0, imax-1
            xyr_RadUwFlux(i,j,k) =                                      &
              &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
              & + aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,4,k) &
              & + xyz_CUpT(i,j,k)
            xyr_RadDwFlux(i,j,k) =                                      &
              &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
              & + aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,2,k) &
              & + xyz_CDoT(i,j,k)
          end do
        end do
      end do





        ! Code for debug
        !
!!$        do k = 1, kmax
!!$          do j = js, je
!!$            do i = 0, imax-1
!!$              xyr_UwFluxDebug(i,j,k-1) =                                         &
!!$                &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
!!$                & - aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,2,k) &
!!$                & + xyz_CUpB(i,j,k)
!!$              xyr_DwFluxDebug(i,j,k-1) =                                         &
!!$                &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
!!$                & - aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,4,k) &
!!$                & + xyz_CDoB(i,j,k)
!!$            end do
!!$          end do
!!$        end do
!!$
!!$        k = kmax
!!$        do j = js, je
!!$          do i = 0, imax-1
!!$            xyr_UwFluxDebug(i,j,k) =                                       &
!!$              &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
!!$              & + aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,4,k) &
!!$              & + xyz_CUpT(i,j,k)
!!$            xyr_DwFluxDebug(i,j,k) =                                       &
!!$              &   aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
!!$              & + aa_Vec( i+imax*(j-1)+1, 2*k  ) * xyaz_smalle(i,j,2,k) &
!!$              & + xyz_CDoT(i,j,k)
!!$          end do
!!$        end do
!!$
!!$
!!$        i = imax/2
!!$        j = jmax/2
!!$        do k = kmax, 0, -1
!!$          write( 6, * ) k, xyr_UwFlux(i,j,k), xyr_UwFluxDebug(i,j,k), xyr_UwFlux(i,j,k) - xyr_UwFluxDebug(i,j,k)
!!$        end do
!!$        do k = kmax, 0, -1
!!$          write( 6, * ) k, xyr_DwFlux(i,j,k), xyr_DwFluxDebug(i,j,k), xyr_DwFlux(i,j,k) - xyr_DwFluxDebug(i,j,k)
!!$        end do
!!$        k = 0
!!$        write( 6, * ) k, xyr_UwFlux(i,j,k), xy_SurfAlbedo(i,j) * xyr_DwFlux(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_TransDirAdj(i,j,0) * xy_cosSZA(i,j), &
!!$          xyr_UwFlux(i,j,k) - ( xy_SurfAlbedo(i,j) * xyr_DwFlux(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_TransDirAdj(i,j,0) * xy_cosSZA(i,j) )
!!$        stop


    else

      ! Source function technique described by Toon et al. [1989] 
      ! is used to calculated infrared heating rate. 
      !
      do k = 0, kmax
        do j = js, je
          xyr_RadUwFlux(:,j,k) = 0.0_DP
          xyr_RadDwFlux(:,j,k) = 0.0_DP
        end do
      end do

      do k = 0, kmax
        do j = js, je
          xyra_DelRadUwFlux(:,j,k,0) = 0.0_DP
          xyra_DelRadUwFlux(:,j,k,1) = 0.0_DP
          xyra_DelRadDwFlux(:,j,k,0) = 0.0_DP
          xyra_DelRadDwFlux(:,j,k,1) = 0.0_DP
        end do
      end do

      do l = 1, NGaussQuad
        Mu = a_GQP( l )

        k = kmax
        do j = js, je
          xyr_IDw(:,j,k) = 0.0_DP
        end do
        do k = kmax, 1, -1
          do j = js, je
            do i = 0, imax-1
              FactJ = &
                & ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) + aa_Vec( i+imax*(j-1)+1, 2*k ) ) &
                &   * xyz_LGamma(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0_DP / xyz_Mu1(i,j,k) )
              FactK = &
                & ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) - aa_Vec( i+imax*(j-1)+1, 2*k ) ) &
                &   * ( 1.0_DP / xyz_Mu1(i,j,k) - xyz_Lambda(i,j,k) )
              Sig1  = &
                & 2.0_DP * (   xyz_B0(i,j,k) &
                &            - xyz_B1(i,j,k) * ( 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) - xyz_Mu1(i,j,k) ) )
              Sig2  = &
                & 2.0_DP * xyz_B1(i,j,k)
              xyr_IDw(i,j,k-1) = xyr_IDw(i,j,k) * exp( - xyz_DelTau(i,j,k) / Mu )           &
                & + FactJ / ( xyz_Lambda(i,j,k) * Mu + 1.0_DP )                             &
                &   * (   1.0_DP                                                            &
                &       - exp( - xyz_DelTau(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0d0 / Mu ) ) ) &
                & + FactK / ( xyz_Lambda(i,j,k) * Mu - 1.0_DP )                             &
                &   * (   exp( - xyz_DelTau(i,j,k) / Mu )                                   &
                &       - exp( - xyz_DelTau(i,j,k) * xyz_Lambda(i,j,k) ) )                  &
                & + Sig1 * ( 1.0_DP - exp( - xyz_DelTau(i,j,k) / Mu ) )                     &
                & + Sig2 * ( Mu * exp( - xyz_DelTau(i,j,k) / Mu ) + xyz_DelTau(i,j,k) - Mu )
            end do
          end do
        end do

        k = 0
        do j = js, je
!               gemis = 1.0d0
!!$          gemis = emis( ij, 1 )
!!$          inteup( ij, 1, k ) = ( 1.0d0 - gemis ) * intedo( ij, 1, km+1 ) &
!!$                    + gemis * pix2 * pfinth( ij, 1, km+1 )
          xyr_IUw(:,j,k) = xy_SurfAlbedo(:,j) * xyr_IDw(:,j,0) + 2.0_DP * xy_SurfPFInted(:,j)
        end do
        if ( FlagTOAFlux ) then
          ! Add of Direct Solar Insident
          do j = js, je
            xyr_IUw(:,j,k) = xyr_IUw(:,j,k) &
              & + xy_SurfAlbedo(:,j) &
              & * xy_SolarFluxTOA(:,j) * xyr_TransDirAdj(:,j,k) * xy_cosSZA(:,j) &
              & * 2.0_DP
          end do
        end if
        do k = 1, kmax
          do j = js, je
            do i = 0, imax-1
              FactG = &
                & ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) + aa_Vec( i+imax*(j-1)+1, 2*k ) ) &
                &   * ( 1.0_DP / xyz_Mu1(i,j,k) - xyz_Lambda(i,j,k) )
              FactH = &
                & ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) - aa_Vec( i+imax*(j-1)+1, 2*k ) ) &
                &   * xyz_LGamma(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0_DP / xyz_Mu1(i,j,k) )
              Alp1  = &
                & 2.0_DP * (   xyz_B0(i,j,k) &
                &            + xyz_B1(i,j,k) * ( 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) - xyz_Mu1(i,j,k) ) )
              Alp2  = &
                & 2.0_DP * xyz_B1(i,j,k)
              xyr_IUw(i,j,k) = xyr_IUw(i,j,k-1) * exp( - xyz_DelTau(i,j,k) / Mu )           &
                & + FactG / ( xyz_Lambda(i,j,k) * Mu - 1.0_DP )                             &
                &   * (   exp( - xyz_DelTau(i,j,k) / Mu )                                   &
                &       - exp( - xyz_DelTau(i,j,k) * xyz_Lambda(i,j,k) ) )                  &
                & + FactH / ( xyz_Lambda(i,j,k) * Mu + 1.0_DP )                             &
                &   * (   1.0_DP                                                            &
                &       - exp( - xyz_DelTau(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0d0 / Mu ) ) ) &
                & + Alp1 * ( 1.0_DP - exp( - xyz_DelTau(i,j,k) / Mu ) )                     &
                & + Alp2 * ( Mu - ( xyz_DelTau(i,j,k) + Mu ) * exp( - xyz_Deltau(i,j,k) / Mu ) )
            end do
          end do
        end do

        do k = 0, kmax
          do j = js, je
            xyr_RadUwFlux(:,j,k) = xyr_RadUwFlux(:,j,k) + Mu * xyr_IUw(:,j,k) * a_GQW( l )
            xyr_RadDwFlux(:,j,k) = xyr_RadDwFlux(:,j,k) + Mu * xyr_IDw(:,j,k) * a_GQW( l )
          end do
        end do

        do k = 0, kmax
          do j = js, je
            xyra_DelRadUwFlux(:,j,k,0) = xyra_DelRadUwFlux(:,j,k,0)             &
              & + Mu * 2.0_DP * xy_SurfDPFDTInted(:,j)                          &
              &   * exp( - ( xyr_OpDepAdj(:,j,0) - xyr_OpDepAdj(:,j,k) ) / Mu ) &
              &   * a_GQW( l )
          end do
        end do

      end do ! do l = 1, NGaussQuad

    end if


    if ( FlagTOAFlux ) then
      !
      ! Addition of Direct Solar Insident
      !
      do k = 0, kmax
        do j = js, je
          xyr_RadDwFlux(:,j,k) = xyr_RadDwFlux(:,j,k) &
            & + xy_SolarFluxTOA(:,j) * xyr_TransDirAdj(:,j,k) * xy_cosSZA(:,j)
        end do
      end do

    end if


  end subroutine RadRTETwoStreamAppCore

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

  !******************************************************************************
  !      subroutine tridiag
  !      tidiagonal solver
  !******************************************************************************
  !     a(j), b(j), and c(j) are, respectively, the subdiagonal, diagonal, 
  !     and superdiagonal entries in row j. 
  !     a(1) and c(jmx) need not be initialized. 
  !     The output is in f; a, b, and c are unchanged. 
  !******************************************************************************
  !     jmx    : dimensions of all the following arrays
  !     a      : sub (lower) diagonal
  !     b      : center diagonal
  !     c      : super (upper) diagonal
  !     f      : right hand side
  !******************************************************************************

  subroutine tridiag( mm, jmx, a, b, c, f, ms, me )

    integer , intent(in   ) :: mm, jmx
    real(DP), intent(in   ) :: a( mm, jmx ),b( mm, jmx )
    real(DP), intent(in   ) :: c( mm, jmx )
    real(DP), intent(inout) :: f( mm, jmx )
    integer , intent(in   ) :: ms, me


    ! Local variables
    !
    real(DP) :: q( mm, jmx ), p
    integer  :: j, m


    ! Forward elimination sweep
    !
    do m = ms, me
      q( m, 1 ) = - c( m, 1 ) / b( m, 1 )
      f( m, 1 ) =   f( m, 1 ) / b( m, 1 )
    end do

    do j = 2, jmx
      do m = ms, me
        p         = 1.0d0 / ( b( m, j ) + a( m, j ) * q( m, j-1 ) )
        q( m, j ) = - c( m, j ) * p
        f( m, j ) = ( f( m, j ) - a( m, j ) * f( m, j-1 ) ) * p
      end do
    end do

    ! Backward pass
    !
    do j = jmx - 1, 1, -1
      do m = ms, me
        f( m, j ) = f( m, j ) + q( m, j ) * f( m, j+1 )
      end do
    end do

  end subroutine tridiag

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

  subroutine tridiag1( jmx, a, b, c, f )

    integer , intent(in   ) :: jmx
    real(DP), intent(in   ) :: a(jmx),b(jmx)
    real(DP), intent(in   ) :: c(jmx)
    real(DP), intent(inout) :: f(jmx)


    ! Local variables
    !
    real(DP) :: q(jmx), p
    integer  :: j


    ! Forward elimination sweep
    !
    q( 1 ) = - c( 1 ) / b( 1 )
    f( 1 ) =   f( 1 ) / b( 1 )

    do j = 2, jmx
      p      = 1.0d0 / ( b( j ) + a( j ) * q( j-1 ) )
      q( j ) = - c( j ) * p
      f( j ) = ( f( j ) - a( j ) * f( j-1 ) ) * p
    end do

    ! Backward pass
    !
    do j = jmx - 1, 1, -1
      f( j ) = f( j ) + q( j ) * f( j+1 )
    end do

  end subroutine tridiag1

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

  subroutine RadRTETwoStreamAppInit

!!$    ! ファイル入出力補助
!!$    ! File I/O support
!!$    !
!!$    use dc_iounit, only: FileOpen
!!$
!!$    ! NAMELIST ファイル入力に関するユーティリティ
!!$    ! Utilities for NAMELIST file input
!!$    !
!!$    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ガウス重み, 分点の計算
    ! 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_two_stream_app_nml/ &
!!$      & NGaussQuad
          !
          ! デフォルト値については初期化手続 "rad_rte_two_stream_app#RadRTETwoStreamAppInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_rte_two_stream_app#RadRTETwoStreamAppInit" for the default values.
          !

    if ( rad_rte_two_stream_app_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
!!$    NGaussQuad = 8


    ! 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_two_stream_app_nml,    & ! (out)
!!$        & iostat = iostat_nml )                  ! (out)
!!$      close( unit_nml )
!!$
!!$      call NmlutilMsg( iostat_nml, module_name ) ! (in)
!!$    end if


!!$    allocate( a_GQP(1:NGaussQuad) )
!!$    allocate( a_GQW(1:NGaussQuad) )

    call GauLeg( &
      & 0.0_DP, 1.0_DP, NGaussQuad, & ! (in )
      & a_GQP, a_GQW                & ! (out)
      & )


    ! Initialization of modules used in this module
    !


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


    rad_rte_two_stream_app_inited = .true.


  end subroutine RadRTETwoStreamAppInit

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

!******************************************************************************
!!$
!!$    subroutine twostreamapp_vis( cossza, gph, q, gdf, &
!!$         dod067, qerat, ssa, af, sf, albedo, &
!!$         gor, goru, gord, gsr, gsru, gsrd, &
!!$         ijs, ije )
!!$
!!$      use matype
!!$      use maparam, only : im => imax, jm => jmax, km => kmax
!!$
!!$      real(dp)    , intent(in ) :: cossza( ijs:ije, 1 )
!!$      real(dp)    , intent(in ) :: gph   ( im, jm, km+1 )
!!$      real(dp)    , intent(out) :: q     ( im, jm, km   )
!!$      real(dp)    , intent(out) :: gdf   ( im, jm )
!!$      real(dp)    , intent(in ) :: dod067( im, jm, km+1 )
!!$      real(dp)    , intent(in ) :: qerat, ssa, af
!!$      real(dp)    , intent(in ) :: sf
!!$      real(dp)    , intent(in ) :: albedo( im, jm )
!!$
!!$      real(dp)    , intent(out) :: &
!!$           gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$           gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )
!!$
!!$      integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      real(dp)     :: gt  ( im, jm, km )
!!$      real(dp)     :: gts ( im, jm )
!!$      real(dp)     :: emis( im, jm )
!!$      real(dp)     :: wn1 = 1.0d100, wn2 = 1.0d100
!!$      integer(i4b) :: divnum = 3
!!$      integer(i4b) :: sw = 1
!!$
!!$      integer(i4b) :: ij, k
!!$
!!$
!!$      do k = 1, km
!!$         do ij = ijs, ije
!!$            gt  ( ij, 1, k ) = 1.0d100
!!$         end do
!!$      end do
!!$      do ij = ijs, ije
!!$         gts ( ij, 1 ) = 1.0d100
!!$         emis( ij, 1 ) = 1.0d100
!!$      end do
!!$
!!$      call twostreamapp( cossza, gt, gts, gph, q, gdf, &
!!$           dod067, qerat, ssa, af, sf, albedo, emis, wn1, wn2, divnum, sw, &
!!$           gor, goru, gord, gsr, gsru, gsrd, &
!!$           ijs, ije )
!!$
!!$
!!$    end subroutine twostreamapp_vis
!!$
!!$!******************************************************************************
!!$
!!$    subroutine twostreamapp_ir( gt, gts, gph, q, gdf, &
!!$         dod067, qerat, ssa, af, emis, wn1, wn2, divnum, &
!!$         gor, goru, gord, gsr, gsru, gsrd, &
!!$         ijs, ije )
!!$
!!$      use matype
!!$      use maparam, only : im => imax, jm => jmax, km => kmax
!!$
!!$      implicit none
!!$
!!$      real(dp)    , intent(in ) :: gt    ( im, jm, km   )
!!$      real(dp)    , intent(in ) :: gts   ( im, jm )
!!$      real(dp)    , intent(in ) :: gph   ( im, jm, km+1 )
!!$      real(dp)    , intent(out) :: q     ( im, jm,   km )
!!$      real(dp)    , intent(out) :: gdf   ( im, jm )
!!$      real(dp)    , intent(in ) :: dod067( im, jm, km+1 )
!!$      real(dp)    , intent(in ) :: qerat, ssa, af
!!$      real(dp)    , intent(in ) :: emis  ( im, jm )
!!$      real(dp)    , intent(in ) :: wn1, wn2
!!$      integer(i4b), intent(in ) :: divnum
!!$
!!$      real(dp)    , intent(out) :: &
!!$           gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$           gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )
!!$
!!$      integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      real(dp)     :: cossza( ijs:ije, 1 )
!!$      real(dp)     :: sf    = 1.0d100
!!$      real(dp)     :: albedo( im, jm )
!!$      integer(i4b) :: sw = 2
!!$
!!$      integer(i4b) :: ij
!!$
!!$
!!$      do ij = ijs, ije
!!$         albedo( ij, 1 ) = 1.0d0 - emis( ij, 1 )
!!$         cossza( ij, 1 ) = 1.0d100
!!$      end do
!!$
!!$      call twostreamapp( cossza, gt, gts, gph, q, gdf, &
!!$           dod067, qerat, ssa, af, sf, albedo, emis, wn1, wn2, divnum, sw, &
!!$           gor, goru, gord, gsr, gsru, gsrd, &
!!$           ijs, ije )
!!$
!!$
!!$    end subroutine twostreamapp_ir

end module rad_rte_two_stream_app