| Class | rad_rte_two_stream_app |
| In: |
radiation/rad_rte_two_stream_app.f90
|
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.
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.
| !$ ! 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#radiation_DennouAGCM_nml
| 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 |
subroutine OLD_RadRTETwoStreamAppHomogAtm( 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 OLD_RadRTETwoStreamAppHomogAtm
| 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 )
| ||
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax) : | 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) |
Alias for RadRTETwoStreamAppWrapper
| 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 )
| ||
| xyr_OptDep( 0:imax-1, 1:jmax, 0:kmax ) : | 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 ) |
Alias for RadRTETwoStreamAppCore
| 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 |
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.
! メッセージ出力
! 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_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
! 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 : | |||
| 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 )
| ||
| xyr_OptDep( 0:imax-1, 1:jmax, 0:kmax ) : | 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 ) |
subroutine RadRTETwoStreamAppCore( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadUwFlux, xyr_RadDwFlux, 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_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 ) :: 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 )
!
! 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
! 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_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_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_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_RadDwFlux(:,j,k) = xyr_RadDwFlux(:,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
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_RadUwFlux(i,j,k) = 0.0d0
xyr_RadDwFlux(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 RadRTETwoStreamAppCore
| 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 )
| ||
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax) : | 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) |
subroutine RadRTETwoStreamAppWrapper( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadUwFlux, xyr_RadDwFlux )
! 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_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadDwFlux(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_RadUwFlux, &
!$OMP xyr_RadDwFlux)
!$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( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadUwFlux, xyr_RadDwFlux, js, je )
end do
!$OMP END DO
!$OMP END PARALLEL
deallocate( a_js )
deallocate( a_je )
end subroutine RadRTETwoStreamAppWrapper
| Constant : | |||
| module_name = ‘rad_rte_two_stream_app‘ : | character(*), parameter
|
| 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 ) |
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) |
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
| Constant : | |||
| version = ’$Name: dcpam5-20120229 $’ // ’$Id: rad_rte_two_stream_app.f90,v 1.3 2012-02-01 05:19:53 yot Exp $’ : | character(*), parameter
|