Class rad_rte_two_stream_app
In: radiation/rad_rte_two_stream_app.f90

Solve radiative transfer equation in two stream approximation

Note that Japanese and English are described in parallel.

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

Methods

Included Modules

dc_types dc_message constants0 gridset gauss_quad

Public Instance methods

Subroutine :
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(in )
xyr_RadSUwFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
xyr_RadSDwFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
FlagSemiInfAtm :logical , intent(in ), optional
FlagSL09 :logical , intent(in ), optional

[Source]

  subroutine RadRTETwoStreamAppHomogAtm( xy_SurfAlbedo, SolarFluxAtTOA, xy_InAngle, SSA, AF, xyr_OptDep, xyr_RadSUwFlux, xyr_RadSDwFlux, FlagSemiInfAtm, FlagSL09 )

    ! 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 :

[Source]

  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, a_GQP, a_GQW )


    ! 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 :
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(in )
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) :real(DP), intent(out)

[Source]

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

    ! 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, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, xyr_PFInted = xyr_PFInted, xy_SurfPFInted = xy_SurfPFInted, xy_SurfDPFDTInted = xy_SurfDPFDTInted, xyra_DelRadUwFlux = xyra_DelRadUwFlux, xyra_DelRadDwFlux = xyra_DelRadDwFlux )


  end subroutine RadRTETwoStreamAppLW
Subroutine :
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) :real(DP), intent(in )
: sec (入射角). sec (angle of incidence)
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)

[Source]

  subroutine RadRTETwoStreamAppSW( xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, SolarFluxTOA, xy_InAngle, xyr_RadUwFlux, xyr_RadDwFlux )

    ! 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, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, SolarFluxTOA = SolarFluxTOA, xy_InAngle = xy_InAngle )


  end subroutine RadRTETwoStreamAppSW

Private Instance methods

DelTauThreshold
Constant :
DelTauThreshold = 1.0e-10_DP :real(DP), parameter
IDScatApproxEddington
Constant :
IDScatApproxEddington = 11 :integer, parameter
IDScatApproxHemiMean
Constant :
IDScatApproxHemiMean = 12 :integer, parameter
NGaussQuad
Constant :
NGaussQuad = 8 :integer, parameter
Subroutine :
IDScatApprox :integer , intent(in )
FlagTOAFlux :logical , intent(in )
FlagEmis :logical , intent(in )
FlagSrcFuncTech :logical , 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_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)
js :integer , intent(in )
je :integer , intent(in )
SolarFluxTOA :real(DP), intent(in ), optional
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(in ), optional
: sec (入射角). sec (angle of incidence)
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(in ), 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(out), optional

[Source]

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

    ! 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
        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
        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 :
IDScatApprox :integer , intent(in )
FlagTOAFlux :logical , intent(in )
FlagEmis :logical , intent(in )
FlagSrcFuncTech :logical , 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_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)
SolarFluxTOA :real(DP), intent(in ), optional
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(in ), optional
: sec (入射角). sec (angle of incidence)
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(in ), 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(out), optional

[Source]

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

    ! 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, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, js, je, SolarFluxTOA = SolarFluxTOA, xy_InAngle = xy_InAngle, xyr_PFInted = xyr_PFInted, xy_SurfPFInted = xy_SurfPFInted, xy_SurfDPFDTInted = xy_SurfDPFDTInted, xyra_DelRadUwFlux = xyra_DelRadUwFlux, xyra_DelRadDwFlux = xyra_DelRadDwFlux )

    end do


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


    deallocate( a_js )
    deallocate( a_je )



  end subroutine RadRTETwoStreamAppWrapper
a_GQP
Variable :
a_GQP(1:NGaussQuad) :real(DP), save
a_GQW
Variable :
a_GQW(1:NGaussQuad) :real(DP), save
module_name
Constant :
module_name = ‘rad_rte_two_stream_app :character(*), parameter
: モジュールの名称. Module name
rad_rte_two_stream_app_inited
Variable :
rad_rte_two_stream_app_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag
Subroutine :
mm :integer , intent(in )
jmx :integer , intent(in )
a( mm, jmx ) :real(DP), intent(in )
b( mm, jmx ) :real(DP), intent(in )
c( mm, jmx ) :real(DP), intent(in )
f( mm, jmx ) :real(DP), intent(inout)
ms :integer , intent(in )
me :integer , intent(in )

[Source]

  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 :
jmx :integer , intent(in )
a(jmx) :real(DP), intent(in )
b(jmx) :real(DP), intent(in )
c(jmx) :real(DP), intent(in )
f(jmx) :real(DP), intent(inout)

[Source]

  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
version
Constant :
version = ’$Name: $’ // ’$Id: rad_rte_two_stream_app.f90,v 1.7 2015/03/11 04:48:47 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version