subroutine RadTwoStreamAppHomogAtm( xy_SurfAlbedo, SolarFluxAtTOA, xy_InAngle, SSA, AF, xyr_OptDep, xyr_RadSFlux, 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.


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

    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_RadSFlux(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


    ! 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


    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
    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
          end if
        end do
      end do
    end do


  end subroutine RadTwoStreamAppHomogAtm

  subroutine RadTwoStreamAppCore( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadFlux, js, je )

    ! USE statements
    !

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

!!$  use pf_module , only : pfint_gq_array


    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 ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
    integer , intent(in ) :: IDScheme
    real(DP), intent(in ) :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ) :: xyr_OptDep    ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadFlux   ( 0:imax-1, 1:jmax, 0:kmax )

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


!!$    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 )

    !
    ! ssa2    : Delta-Function Adjusted Single-Scattering Albedo
    ! af2     : Delta-Function Adjusted Asymmetry Factor
    ! opdep2  : Delta-Function Adjusted Optical Depth
    !
    real(DP) :: xyz_SSA2  ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyz_AF2   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP) :: xyr_OpDep2( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: xyr_Trans2( 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 )

    !
    ! upfl      : Upward Flux
    ! dofl      : Downward Flux
    !
    real(DP) :: xyr_UwFlux( 0:imax-1, 1:jmax, 0:kmax )
    real(DP) :: xyr_DwFlux( 0:imax-1, 1:jmax, 0: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_DTau( 0:imax-1, 1:jmax, 1:kmax )

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


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

    !
    ! 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) :: 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) :: mu1
!!$  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) :: inteup( 0:imax-1, 1:jmax, 0:kmax )
!!$    real(DP) :: intedo( 0:imax-1, 1:jmax, 0:kmax )
!!$    real(DP) :: tmpg, tmph, tmpj, tmpk, alp1, alp2, sig1, sig2
!!$    real(DP) :: tmpcos
!!$!      integer(i4b), parameter :: quadn = 3
!!$    integer(i4b), parameter :: quadn = 5
!!$!      integer(i4b), parameter :: quadn = 8
!!$    real(DP) :: qucos ( quadn ), qudcos( quadn )


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


!!$      call gauleg( 0.0d0, 1.0d0, quadn, qucos, qudcos )

    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( IDScheme .eq. IDSchemeShortWave ) 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


!!$      gth(:,:,:) = 1.0d100

!!$    else
!!$
!!$      stop 'sw != 1 is not supported.'
!!$
!!$         do ij = ijs, ije
!!$            cosz ( ij, 1 ) = 1.0d100
!!$            cosz2( ij, 1 ) = 1.0d100
!!$         end do
!!$
!!$         do ij = ijs, ije
!!$            gth( ij, 1, 1 ) = gt( ij, 1, 1 )
!!$         end do
!!$         do k = 2, km+1-1
!!$            do ij = ijs, ije
!!$               gth( ij, 1, k ) = ( gt( ij, 1, k-1 ) + gt( ij, 1, k ) ) * 0.5d0
!!$            end do
!!$         end do
!!$         do ij = ijs, ije
!!$            gth( ij, 1, km+1 ) = gts( ij, 1 )
!!$         end do
!!$!         call pfint_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$!              ijs, ije )
!!$         call pfint_gq_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$              ijs, ije )

    else
      stop 'Unexpected sw number in twostreamapp.f90'
    end if


    if( IDScheme .eq. IDSchemeShortWave ) then

      !
      ! Delta-Function Adjustment
      !
      do k = 1, kmax
        do j = js, je
          xyz_AF2 (:,j,k) = xyz_AF(:,j,k) / ( 1.0d0 + xyz_AF(:,j,k) )
          xyz_SSA2(:,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

!!$      opdep2(:,:,:) = ( 1.0d0 - xyz_ssa * xyz_af**2 ) * opdep(:,:,:)
      do k = 0, kmax
        do j = js, je
          xyr_OpDep2(:,j,k) = 0.0d0
        end do
      end do
      do k = kmax-1, 0, -1
        do j = js, je
          xyr_OpDep2(:,j,k) = xyr_OpDep2(:,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

      do k = 0, kmax
        do j = js, je
          xyr_Trans2(:,j,k) = exp( -xyr_OpDep2(:,j,k) * xy_cosSZAInv(:,j) )
        end do
      end do

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


!!$    else if( sw .eq. 2 ) then
!!$
!!$         !
!!$         ! In infrared, delta-function adjustment is not done. 
!!$         !
!!$         af2  = af
!!$         ssa2 = ssa
!!$
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               opdep2( ij, 1, k ) = opdep( ij, 1, k )
!!$            end do
!!$         end do
!!$
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               exp_opdep2( ij, 1, k ) = 1.0d100
!!$            end do
!!$         end do
!!$
!!$         !
!!$         ! Hemispheric mean approximation
!!$         !
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               gam1( ij, 1, k ) = 2.0d0 - ssa2 * ( 1.0d0 + af2 )
!!$               gam2( ij, 1, k ) = ssa2 * ( 1.0d0 - af2 )
!!$               gam3( ij, 1, k ) = 1.0d100
!!$               gam4( ij, 1, k ) = 1.0d100
!!$            end do
!!$         end do

    else
      stop 'Unexpected sw number in twostreamapp.f90'
    end if


    do k = 1, kmax
      do j = js, je
        xyz_DTau(:,j,k) = xyr_OpDep2(:,j,k-1) - xyr_OpDep2(:,j,k)
      end do
    end do
    !
    ! Avoiding singularity when dtau equal to zero 
    !
    do k = 1, kmax
      do j = js, je
        do i = 0, imax-1
          if( ( IDScheme .eq. IDSchemeLongWave ) .and. ( xyz_DTau(i,j,k) .lt. 1.0d-10 ) ) then
            xyz_DTau(i,j,k) = 0.0d0
          end if
        end do
      end do
    end do


    do k = 1, kmax
      do j = js, je
        xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) )
        xyz_LSigma(:,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_DTau(:,j,k) )
        xyaz_smalle(:,j,1,k) = xyz_LSigma(:,j,k) * xy_TMPVal(:,j) + 1.0_DP
        xyaz_smalle(:,j,2,k) = xyz_LSigma(:,j,k) * xy_TMPVal(:,j) - 1.0_DP
        xyaz_smalle(:,j,3,k) = xy_TMPVal(:,j) + xyz_LSigma(:,j,k)
        xyaz_smalle(:,j,4,k) = xy_TMPVal(:,j) - xyz_LSigma(:,j,k)
      end do
    end do


    if( IDScheme .eq. IDSchemeShortWave ) then
      do k = 1, kmax
        do j = js, je
          ! Lines below will be deleted in future (yot, 2010/09/14).
!!$          xyz_cupb(:,j,k) =                                                       &
!!$            &   xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1)              &
!!$            &   * ( ( 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_cupt(:,j,k) =                                                       &
!!$            &   xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k  )              &
!!$            &   * ( ( 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_cdob(:,j,k) =                                                       &
!!$            &   xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1)              &
!!$            &   * ( ( 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_cdot(:,j,k) =                                                       &
!!$            &   xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k  )              &
!!$            &   * ( ( 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) )


          xy_TMPVal(:,j) = xyz_SSA2(:,j,k) * SolarFluxTOA * ( ( 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_cupb(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k-1)
          xyz_cupt(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k  )

          xy_TMPVal(:,j) = xyz_SSA2(:,j,k) * SolarFluxTOA * ( ( 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_cdob(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k-1)
          xyz_cdot(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k  )
        end do
      end do

!!$      else if( sw .eq. 2 ) then
!!$
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               !
!!$               ! Notice!
!!$               ! Avoiding singularity when dtau equal to zero. 
!!$               ! dtau occationally has much smaller value. 
!!$               ! When this occurs, b1 cannot be calculated correctly. 
!!$               !
!!$               if( dtau( ij, 1, k ) .ne. 0.0d0 ) then
!!$                  b0( ij, 1, k ) = pfinth( ij, 1, k )
!!$                  b1( ij, 1, k ) = ( pfinth( ij, 1, k+1 ) - pfinth( ij, 1, k ) ) / dtau( ij, 1, k )
!!$               else
!!$                  b0( ij, 1, k ) = 0.0d0
!!$                  b1( ij, 1, k ) = 0.0d0
!!$               end if
!!$
!!$            end do
!!$         end do
!!$
!!$         do k = 1, km
!!$            do ij = ijs, ije
!!$               mu1 = ( 1.0d0 - ssa2 ) / ( gam1( ij, 1, k ) - gam2( ij, 1, k ) )
!!$
!!$               cupb( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( dtau( ij, 1, k ) + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cupt( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( 0.0d0     + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cdob( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( dtau( ij, 1, k ) - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$               cdot( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$                    + b1( ij, 1, k ) &
!!$                    * ( 0.0d0     - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$            end do
!!$         end do
!!$
    else
      stop 'Unexpected sw number in twostreamapp.f90'
    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
    if( IDScheme .eq. IDSchemeShortWave ) then
      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_SurfAlbedo(i,j) * SolarFluxTOA * xyr_Trans2(i,j,0) * xy_cosSZA(i,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
      stop 'Unexpected sw number in twostreamapp.f90'
    end if


    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( IDScheme .eq. IDSchemeShortWave ) then

      k = 1
      do j = js, je
        do i = 0, imax-1
          xyr_UwFlux(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_DwFlux(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_UwFlux(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_DwFlux(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_Trans2(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_Trans2(i,j,0) * xy_cosSZA(i,j) )
!!$        stop


      !
      ! Addition of Direct Solar Insident
      !
      do k = 0, kmax
        do j = js, je
          xyr_DwFlux(:,j,k) = xyr_DwFlux(:,j,k) + SolarFluxTOA * xyr_Trans2(:,j,k) * xy_cosSZA(:,j)
        end do
      end do

!!$      else if( sw .eq. 2 ) then
!!$
!!$         ! Source function technique described by Toon et al. [1989] 
!!$         ! is used to calculated infrared heating rate. 
!!$         !
!!$         do k = 1, km+1
!!$            do ij = ijs, ije
!!$               upfl( ij, 1, k ) = 0.0d0
!!$               dofl( ij, 1, k ) = 0.0d0
!!$            end do
!!$         end do
!!$
!!$         do l = 1, quadn
!!$            tmpcos = qucos( l )
!!$
!!$            k = 1
!!$            do ij = ijs, ije
!!$               intedo( ij, 1, k ) = 0.0d0
!!$            end do
!!$            do k = 1, km
!!$               do ij = ijs, ije
!!$                  tmpj = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 2.0d0 )
!!$                  tmpk = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) &
!!$                       * ( 2.0d0 - lambda( ij, 1, k ) )
!!$                  sig1 = pix2 * ( b0( ij, 1, k ) &
!!$                       - b1( ij, 1, k ) &
!!$                       * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$                  sig2 = pix2 * b1( ij, 1, k )
!!$                  intedo( ij, 1, k+1 ) = intedo( ij, 1, k ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + tmpj / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$                       * ( 1.0d0 &
!!$                       - exp( -dtau( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$                       + tmpk / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$                       * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       -exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$                       + sig1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$                       + sig2 * ( tmpcos * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + dtau( ij, 1, k ) - tmpcos )
!!$               end do
!!$            end do
!!$
!!$            k = km+1
!!$            do ij = ijs, ije
!!$!               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 )
!!$            end do
!!$            do k = km, 1, -1
!!$               do ij = ijs, ije
!!$                  tmpg = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) &
!!$                       * ( 2.0d0 - lambda( ij, 1, k ) )
!!$                  tmph = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 2.0d0 )
!!$                  alp1 = pix2 * ( b0( ij, 1, k ) &
!!$                       + b1( ij, 1, k ) &
!!$                       * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$                  alp2 = pix2 * b1( ij, 1, k )
!!$                  inteup( ij, 1, k ) = inteup( ij, 1, k+1 ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       + tmpg / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$                       * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$                       - exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$                       + tmph / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$                       * ( 1.0d0 &
!!$                       -exp( -dtau( ij, 1, k ) &
!!$                       * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$                       + alp1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$                       + alp2 * ( tmpcos &
!!$                       - ( dtau( ij, 1, k ) + tmpcos ) &
!!$                       * exp( -dtau( ij, 1, k ) / tmpcos ) )
!!$               end do
!!$            end do
!!$
!!$            do k = 1, km+1
!!$               do ij = ijs, ije
!!$                  upfl( ij, 1, k ) = upfl( ij, 1, k ) &
!!$                       + inteup( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$                  dofl( ij, 1, k ) = dofl( ij, 1, k ) &
!!$                       + intedo( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$               end do
!!$            end do
!!$         end do
!!$
    else
      stop 'Unexpected sw number in twostreamapp.f90'
    end if


!!$      do ij = ijs, ije
!!$         gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$      end do
!!$      if( sw == 1 ) then
!!$        ! visible
!!$        do ij = ijs, ije
!!$          gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$        end do
!!$      else
!!$        ! ir
!!$        do ij = ijs, ije
!!$          gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$        end do
!!$      end if
!!$      do ij = ijs, ije
!!$        gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$      end do


!!$      do k = 1, km
!!$         do ij = ijs, ije
!!$            q( ij, 1, k ) &
!!$                 = ( ( upfl( ij, 1, k+1 ) - dofl( ij, 1, k+1 ) ) &
!!$                   - ( upfl( ij, 1, k   ) - dofl( ij, 1, k   ) ) ) &
!!$                   / ( gph ( ij, 1, k+1 ) - gph ( ij, 1, k   ) ) &
!!$                   * grav / cp
!!$         end do
!!$      end do


!!$      if( sw == 1 ) then
!!$        write( 6, * ) 'vis flux=', &
!!$!          -gdf((ijs+ije)/2,1), &
!!$!          -gdf((ijs+ije)/2,1) * ( 1.0d0 - albedo((ijs+ije)/2,1) ), &
!!$          -gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1), &
!!$!          upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1),
!!$          upfl((ijs+ije)/2,1,km+1), &
!!$          upfl((ijs+ije)/2,1,km+1) - gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1)
!!$      else
!!$        write( 6, * ) 'ir  flux=', &
!!$          -gdf((ijs+ije)/2,1), &
!!$          upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1), &
!!$          ' sig * Ts^4 = ', sboltz * gts((ijs+ije)/2,1)**4
!!$      end if

    do k = 0, kmax
      do j = js, je
        xyr_RadFlux(:,j,k) = xyr_UwFlux(:,j,k) - xyr_DwFlux(:,j,k)
      end do
    end do

    if( IDScheme .eq. IDSchemeShortWave ) then
      do k = 0, kmax
        do j = js, je
          do i = 0, imax-1
            if( xy_cosSZA(i,j) <= 0.0d0 ) then
              xyr_RadFlux(i,j,k) = 0.0d0
            end if
          end do
        end do
      end do
    end if


!!$      !
!!$      ! output variables
!!$      !
!!$      do ij = ijs, ije
!!$         goru( ij, 1 ) = upfl( ij, 1, 1 )
!!$         gord( ij, 1 ) = dofl( ij, 1, 1 )
!!$         gsru( ij, 1 ) = upfl( ij, 1, km+1 )
!!$         gsrd( ij, 1 ) = dofl( ij, 1, km+1 )
!!$      end do
!!$      if( sw .eq. 1 ) then
!!$         do ij = ijs, ije
!!$            if( cosz( ij, 1 ) .le. 0.0d0 ) then
!!$               goru( ij, 1 ) = 0.0d0
!!$               gord( ij, 1 ) = 0.0d0
!!$               gsru( ij, 1 ) = 0.0d0
!!$               gsrd( ij, 1 ) = 0.0d0
!!$            end if
!!$         end do
!!$      end if
!!$      do ij = ijs, ije
!!$         gor ( ij, 1 ) = goru( ij, 1 ) - gord( ij, 1 )
!!$         gsr ( ij, 1 ) = gsru( ij, 1 ) - gsrd( ij, 1 )
!!$      end do


  end subroutine RadTwoStreamAppCore

  subroutine RadTwoStreamAppWrapper( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadFlux )

    ! USE statements
    !

    ! OpenMP
    !
    !$ use omp_lib


    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 ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
    integer , intent(in ) :: IDScheme
    real(DP), intent(in ) :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP), intent(in ) :: xyr_OptDep    ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadFlux   ( 0:imax-1, 1:jmax, 0:kmax )


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

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

    integer :: n


    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        IDScheme, &
    !$OMP        xy_InAngle, &
    !$OMP        xyr_OptDep, &
    !$OMP        xyr_RadFlux)

    !$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 RadTwoStreamAppCore( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadFlux, js, je )

    end do


    !$OMP END DO
    !$OMP END PARALLEL


    deallocate( a_js )
    deallocate( a_je )



  end subroutine RadTwoStreamAppWrapper

  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