Class radiation_two_stream_app
In: radiation/radiation_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 gridset constants

Public Instance methods

IDSchemeLongWave
Constant :
IDSchemeLongWave = 2 :integer, parameter, public
IDSchemeShortWave
Constant :
IDSchemeShortWave = 1 :integer, parameter, public
RadiationTwoStreamApp( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadFlux )
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 )
SolarFluxTOA :real(DP), intent(in )
xy_SurfAlbedo( 0:imax-1, 1:jmax ) :real(DP), intent(in )
IDScheme :integer , intent(in )
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(in )
: sec (入射角). sec (angle of incidence)
xyr_OptDep( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(in )
xyr_RadFlux( 0:imax-1, 1:jmax, 0:kmax ) :real(DP), intent(out)

Alias for RadiationTwoStreamAppWrapper

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

Alias for RadiationTwoStreamAppCore

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

[Source]

  subroutine RadiationTwoStreamAppHomogAtm( 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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    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 RadiationTwoStreamAppHomogAtm

Private Instance methods

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

[Source]

  subroutine RadiationTwoStreamAppCore( 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

    !
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, PI, StB      ! $ \sigma_{SB} $ .
                                  ! ステファンボルツマン定数.
                                  ! Stefan-Boltzmann constant

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

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

[Source]

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

    ! USE statements
    !

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 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 RadiationTwoStreamAppCore( 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 RadiationTwoStreamAppWrapper
module_name
Constant :
module_name = ‘radiation_two_stream_app :character(*), parameter
: モジュールの名称. Module name
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: dcpam5-20101015 $’ // ’$Id: radiation_two_stream_app.f90,v 1.6 2010-10-07 15:38:11 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version