subroutine RadiationDcpamESWV21Flux( xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_Temp, xyz_QH2OVap, xyz_QH2OLiq, xyz_QH2OSol, xyz_QO3, xyr_RadSFlux )
! USE statements
!
!
! Physical constants settings
!
use constants, only: Grav, PI ! $ \pi $ .
! Circular constant
use radiation_two_stream_app, only: RadiationTwoStreamApp
! 時刻管理
! Time control
!
use timeset, only: TimeN, EndTime, TimesetClockStart, TimesetClockStop
! Chou and Lee (1996) による短波放射モデル
! Short wave radiation model described by Chou and Lee (1996)
!
use radiation_CL1996, only : RadiationCL1996NumBands , RadiationCL1996UVVISParams , RadiationCL1996IRH2ONumKDFBin, RadiationCL1996IRH2OKDFParams, RadiationCL1996H2ODelAbsAmt
! 短波入射 (太陽入射)
! Short wave (insolation) incoming
!
use radiation_short_income, only : ShortIncoming
! Chou et al (1998) による短波放射用雲モデル
! Cloud model for short wave radiation model described by Chou et al (1998)
!
use radiation_C1998, only : RadiationC1998CloudParams
real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
real(DP), intent(in ) :: xyz_Press ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyr_Press ( 0:imax-1, 1:jmax, 0:kmax )
real(DP), intent(in ) :: xyz_Temp ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyz_QH2OVap ( 0:imax-1, 1:jmax, 1:kmax )
! $ q $ . 混合比. Mass mixing ratio of constituents (1)
real(DP), intent(in ) :: xyz_QH2OLiq ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyz_QH2OSol (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_QO3 ( 0:imax-1, 1:jmax, 1:kmax )
! O3 分布 (1)
! O3 distribution (1)
real(DP), intent(out) :: xyr_RadSFlux ( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: xy_InAngle (0:imax-1, 1:jmax)
! sec (入射角).
! sec (angle of incidence)
real(DP) :: DistFromStarScld
! Distance between the central star and the planet
real(DP) :: DiurnalMeanFactor
integer :: nbands1
integer :: nbands2
real(DP) :: UVVISFracSolarFlux
real(DP) :: UVVISO3AbsCoef
real(DP) :: UVVISRayScatCoef
real(DP) :: SolarFluxTOA
real(DP) :: xyz_DelAtmMass(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_DelO3Mass (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_SSA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_AF (0:imax-1, 1:jmax, 1:kmax)
real(DP), parameter :: RayScatSinScatAlb = 1.0d0 - 1.0d-10
real(DP), parameter :: RayScatAsymFact = 0.0d0
integer , parameter :: IDScheme = IDSchemeShortWave
real(DP) :: xyz_DelCloudWat ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_DelCloudIce ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_DelCloudWatOptDep( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_DelCloudIceOptDep( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_RayScatDelOptDep( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_O3AbsDelOptDep ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_DelTotOptDep( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyr_TotOptDep ( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: xyr_RadFlux ( 0:imax-1, 1:jmax, 0:kmax )
integer :: nkdf
integer :: ikdfbin
real(DP) :: KDFAbsCoef
real(DP) :: KDFWeight
real(DP) :: xyz_H2ODelAbsAmt( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_H2ODelOptDep( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_CloudREff (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudExtCoef(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudCoAlb (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudWatSSA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudWatAF (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudIceSSA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudIceAF (0:imax-1, 1:jmax, 1:kmax)
integer :: i
integer :: j
integer :: k
integer :: l
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! 初期化
! Initialization
!
if ( .not. radiation_dcpam_E_SW_V2_1_inited ) call RadiationDcpamESWV21Init
! 短波入射の計算
! Calculate short wave (insolation) incoming radiation
!
call ShortIncoming( xy_InAngle = xy_InAngle, DistFromStarScld = DistFromStarScld, DiurnalMeanFactor = DiurnalMeanFactor )
do k = 1, kmax
xyz_DelAtmMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k ) ) / Grav
end do
xyz_DelCloudWat = xyz_DelAtmMass * xyz_QH2OLiq
xyz_DelCloudIce = xyz_DelAtmMass * xyz_QH2OSol
xyz_DelO3Mass = xyz_DelAtmMass * xyz_QO3
call RadiationCL1996NumBands( nbands1, nbands2 )
xyr_RadSFlux = 0.0_DP
! * 14286 to 57143 cm-1 (0.175 to 0.70 micron):
! * Rayleigh scattering,
! * scattering by cloud droplets.
! * O3 absorption
!
do l = 1, nbands1
! Water cloud optical properties
!
xyz_CloudREff = CloudWatREff
call RadiationC1998CloudParams( 'Liquid', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudWatAF )
xyz_DelCloudWatOptDep = xyz_CloudExtCoef * xyz_DelCloudWat
xyz_CloudWatSSA = 1.0_DP - xyz_CloudCoAlb
! Ice cloud optical properties
!
xyz_CloudREff = CloudIceREff
call RadiationC1998CloudParams( 'Ice', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudIceAF )
xyz_DelCloudIceOptDep = xyz_CloudExtCoef * xyz_DelCloudIce
xyz_CloudIceSSA = 1.0_DP - xyz_CloudCoAlb
! UV and Visible optical properties and solar flux
!
call RadiationCL1996UVVISParams( l, UVVISFracSolarFlux, UVVISO3AbsCoef, UVVISRayScatCoef )
SolarFluxTOA = UVVISFracSolarFlux * SolarConst / DistFromStarScld**2 * DiurnalMeanFactor
! Rayleigh scattering
!
if ( FlagRayleighScattering ) then
xyz_RayScatDelOptDep = UVVISRayScatCoef * xyz_DelAtmMass
else
xyz_RayScatDelOptDep = 0.0d0
end if
! O3 absorption
!
xyz_O3AbsDelOptDep = UVVISO3AbsCoef * xyz_DelO3Mass
! Total optical parameter
!
xyz_DelTotOptDep = xyz_DelCloudWatOptDep + xyz_DelCloudIceOptDep + xyz_RayScatDelOptDep + xyz_O3AbsDelOptDep
!
xyr_TotOptDep(:,:,kmax) = 0.0d0
do k = kmax-1, 0, -1
xyr_TotOptDep(:,:,k) = xyr_TotOptDep(:,:,k+1) + xyz_DelTotOptDep(:,:,k+1)
end do
!
xyz_SSA = ( xyz_CloudWatSSA * xyz_DelCloudWatOptDep + xyz_CloudIceSSA * xyz_DelCloudIceOptDep + RayScatSinScatAlb * xyz_RayScatDelOptDep + 0.0d0 * xyz_O3AbsDelOptDep ) / ( xyz_DelTotOptDep + 1.0d-100 )
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_SSA(i,j,k) >= 1.0d0 ) then
xyz_SSA(i,j,k) = 1.0d0 - 1.0d-10
end if
end do
end do
end do
xyz_AF = ( xyz_CloudWatAF * xyz_CloudWatSSA * xyz_DelCloudWatOptDep + xyz_CloudIceAF * xyz_CloudIceSSA * xyz_DelCloudIceOptDep + RayScatAsymFact * RayScatSinScatAlb * xyz_RayScatDelOptDep + 0.0d0 * 0.0d0 * xyz_O3AbsDelOptDep ) / ( xyz_SSA * xyz_DelTotOptDep + 1.0d-100 )
call RadiationTwoStreamApp( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_TotOptDep, xyr_RadFlux )
xyr_RadSFlux = xyr_RadSFlux + xyr_RadFlux
end do
! * 1000 to 14286 cm-1 (0.70-10 micron):
! * absorption by H2O,
! * scattering by cloud droplets.
!
call RadiationCL1996H2ODelAbsAmt( xyz_Temp, xyz_QH2OVap, xyr_Press, xyz_Press, xyz_H2ODelAbsAmt )
call RadiationCL1996IRH2ONumKDFBin( nkdf )
SolarFluxTOA = SolarConst / DistFromStarScld**2 * DiurnalMeanFactor
do l = nbands1+1, nbands1+nbands2
! Water cloud optical properties
!
xyz_CloudREff = CloudWatREff
call RadiationC1998CloudParams( 'Liquid', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudWatAF )
xyz_DelCloudWatOptDep = xyz_CloudExtCoef * xyz_DelCloudWat
xyz_CloudWatSSA = 1.0_DP - xyz_CloudCoAlb
! Ice cloud optical properties
!
xyz_CloudREff = CloudIceREff
call RadiationC1998CloudParams( 'Ice', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudIceAF )
xyz_DelCloudIceOptDep = xyz_CloudExtCoef * xyz_DelCloudIce
xyz_CloudIceSSA = 1.0_DP - xyz_CloudCoAlb
do ikdfbin = 1, nkdf
call RadiationCL1996IRH2OKDFParams( l, ikdfbin, KDFAbsCoef, KDFWeight )
xyz_H2ODelOptDep = KDFAbsCoef * xyz_H2ODelAbsAmt
! Total optical parameter
!
xyz_DelTotOptDep = xyz_DelCloudWatOptDep + xyz_H2ODelOptDep
xyr_TotOptDep(:,:,kmax) = 0.0d0
do k = kmax-1, 0, -1
xyr_TotOptDep(:,:,k) = xyr_TotOptDep(:,:,k+1) + xyz_DelTotOptDep(:,:,k+1)
end do
xyz_SSA = ( xyz_CloudWatSSA * xyz_DelCloudWatOptDep + xyz_CloudIceSSA * xyz_DelCloudIceOptDep + 0.0d0 * xyz_H2ODelOptDep ) / ( xyz_DelTotOptDep + 1.0d-100 )
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_SSA(i,j,k) >= 1.0d0 ) then
xyz_SSA(i,j,k) = 1.0d0 - 1.0d-10
end if
end do
end do
end do
xyz_AF = ( xyz_CloudWatAF * xyz_CloudWatSSA * xyz_DelCloudWatOptDep + xyz_CloudIceAF * xyz_CloudIceSSA * xyz_DelCloudIceOptDep + 0.0d0 * 0.0d0 * xyz_H2ODelOptDep ) / ( xyz_ssa * xyz_DelTotOptDep + 1.0d-100 )
call RadiationTwoStreamApp( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_TotOptDep, xyr_RadFlux )
xyr_RadFlux = xyr_RadFlux * KDFWeight
xyr_RadSFlux = xyr_RadSFlux + xyr_RadFlux
end do
end do
!!$ i = 0
!!$ j = jmax / 2 + 1
!!$ do k = 1, kmax
!!$ write( 73, * ) xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), &
!!$ & xyz_Press(i,j,k)
!!$ end do
!!$ call flush( 73 )
!!$
!!$ i = 0
!!$ j = jmax / 2 + 1
!!$ do k = 1, kmax
!!$ write( 83, * ) &
!!$ & + ( xyr_RadSFlux(i,j,k-1) - xyr_RadSFlux(i,j,k) ) &
!!$ & / ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) &
!!$ & / 1004.6 * Grav, &
!!$ & xyz_Press(i,j,k)
!!$ end do
!!$ call flush( 83 )
!!$
!!$! write( 6, * ) '########## ', acos( 1.0d0 / xy_InAngle(i,j) ) * 180.0d0 / 3.141592d0
!!$
!!$
!!$ i = 0
!!$ j = jmax / 2 + 1
!!$ write( 93, * ) '# ', xy_SurfAlbedo(i,j)
!!$ write( 93, * ) '# ', 1.0_DP / xy_InAngle(i,j)
!!$ do k = 0, kmax
!!$ write( 93, * ) xyr_RadSFlux(i,j,k), xyr_Press(i,j,k)
!!$ end do
!!$ call flush( 93 )
!!$ stop
! 計算時間計測一時停止
! Pause measurement of computation time
!
call TimesetClockStop( module_name )
end subroutine RadiationDcpamESWV21Flux