subroutine RadEarthSWV21Flux( xy_SurfAlbedo, xyz_DelAtmMass, xyz_DelH2OVapMass, xyz_DelH2OLiqMass, xyz_DelH2OSolMass, xyz_DelO3Mass, xyz_Press, xyz_Temp, xyz_CloudCover, xyr_RadSUwFlux, xyr_RadSDwFlux )
! USE statements
!
! 短波入射 (太陽入射)
! Short wave (insolation) incoming
!
use rad_short_income, only : RadShortIncome
!
! Solve radiative transfer equation in two stream approximation
!
use rad_rte_two_stream_app, only: RadRTETwoStreamApp
! Chou and Lee (1996) による短波放射モデル
! Short wave radiation model described by Chou and Lee (1996)
!
use rad_CL1996, only : RadCL1996NumBands , RadCL1996UVVISParams , RadCL1996IRH2ONumKDFBin, RadCL1996IRH2OKDFParams, RadCL1996ScaleH2OVapMass
! Chou et al (1998) による短波放射用雲モデル
! Cloud model for short wave radiation model described by Chou et al (1998)
!
use rad_C1998, only : RadC1998CalcCloudOptProp
!!$ ! 雲の設定
!!$ ! Set cloud
!!$ !
!!$ use set_cloud, only : SetCloudSmearCloudOptDep
real(DP), intent(in ) :: xy_SurfAlbedo (0:imax-1, 1:jmax )
real(DP), intent(in ) :: xyz_DelAtmMass (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_DelH2OVapMass(0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_DelH2OLiqMass(0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_DelH2OSolMass(0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_DelO3Mass (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_CloudCover (0:imax-1, 1:jmax, 1: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)
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_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_DelH2OVapMassScaled( 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 )
real(DP) :: xyr_RadUwFlux ( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: xyr_RadDwFlux ( 0:imax-1, 1:jmax, 0:kmax )
integer :: nkdf
integer :: ikdfbin
real(DP) :: KDFAbsCoef
real(DP) :: KDFWeight
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
! 初期化確認
! Initialization check
!
if ( .not. rad_Earth_SW_V2_1_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! 短波入射の計算
! Calculate short wave (insolation) incoming radiation
!
call RadShortIncome( xy_InAngle = xy_InAngle, DistFromStarScld = DistFromStarScld, DiurnalMeanFactor = DiurnalMeanFactor )
call RadCL1996NumBands( nbands1, nbands2 )
! Initialization
!
xyr_RadSUwFlux = 0.0_DP
xyr_RadSDwFlux = 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 RadC1998CalcCloudOptProp( 'Liquid', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudWatAF )
xyz_DelCloudWatOptDep = xyz_CloudExtCoef * xyz_DelH2OLiqMass
! This is commented out, because xyz_DelCloudWatOptDep is an averaged
! optical depth of a grid box (yot, 2011/11/20).
!!$ call SetCloudSmearCloudOptDep( &
!!$ & xyz_CloudCover, & ! (in )
!!$ & xyz_DelCloudWatOptDep & ! (inout)
!!$ & )
xyz_CloudWatSSA = 1.0_DP - xyz_CloudCoAlb
! Ice cloud optical properties
!
xyz_CloudREff = CloudIceREff
call RadC1998CalcCloudOptProp( 'Ice', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudIceAF )
xyz_DelCloudIceOptDep = xyz_CloudExtCoef * xyz_DelH2OSolMass
! This is commented out, because xyz_DelCloudWatOptDep is an averaged
! optical depth of a grid box (yot, 2011/11/20).
!!$ call SetCloudSmearCloudOptDep( &
!!$ & xyz_CloudCover, & ! (in )
!!$ & xyz_DelCloudIceOptDep & ! (inout)
!!$ & )
xyz_CloudIceSSA = 1.0_DP - xyz_CloudCoAlb
! UV and Visible optical properties and solar flux
!
call RadCL1996UVVISParams( 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 RadRTETwoStreamApp( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_TotOptDep, xyr_RadUwFlux, xyr_RadDwFlux )
xyr_RadSUwFlux = xyr_RadSUwFlux + xyr_RadUwFlux
xyr_RadSDwFlux = xyr_RadSDwFlux + xyr_RadDwFlux
end do
! * 1000 to 14286 cm-1 (0.70-10 micron):
! * absorption by H2O,
! * scattering by cloud droplets.
!
call RadCL1996ScaleH2OVapMass( xyz_Temp, xyz_DelH2OVapMass, xyz_Press, xyz_DelH2OVapMassScaled )
call RadCL1996IRH2ONumKDFBin( nkdf )
SolarFluxTOA = SolarConst / DistFromStarScld**2 * DiurnalMeanFactor
do l = nbands1+1, nbands1+nbands2
! Water cloud optical properties
!
xyz_CloudREff = CloudWatREff
call RadC1998CalcCloudOptProp( 'Liquid', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudWatAF )
xyz_DelCloudWatOptDep = xyz_CloudExtCoef * xyz_DelH2OLiqMass
! This is commented out, because xyz_DelCloudWatOptDep is an averaged
! optical depth of a grid box (yot, 2011/11/20).
!!$ call SetCloudSmearCloudOptDep( &
!!$ & xyz_CloudCover, & ! (in )
!!$ & xyz_DelCloudWatOptDep & ! (inout)
!!$ & )
xyz_CloudWatSSA = 1.0_DP - xyz_CloudCoAlb
! Ice cloud optical properties
!
xyz_CloudREff = CloudIceREff
call RadC1998CalcCloudOptProp( 'Ice', l, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudCoAlb, xyz_CloudIceAF )
xyz_DelCloudIceOptDep = xyz_CloudExtCoef * xyz_DelH2OSolMass
! This is commented out, because xyz_DelCloudWatOptDep is an averaged
! optical depth of a grid box (yot, 2011/11/20).
!!$ call SetCloudSmearCloudOptDep( &
!!$ & xyz_CloudCover, & ! (in )
!!$ & xyz_DelCloudIceOptDep & ! (inout)
!!$ & )
xyz_CloudIceSSA = 1.0_DP - xyz_CloudCoAlb
do ikdfbin = 1, nkdf
call RadCL1996IRH2OKDFParams( l, ikdfbin, KDFAbsCoef, KDFWeight )
xyz_H2ODelOptDep = KDFAbsCoef * xyz_DelH2OVapMassScaled
! Total optical parameter
!
xyz_DelTotOptDep = xyz_DelCloudWatOptDep + xyz_DelCloudIceOptDep + 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 RadRTETwoStreamApp( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_TotOptDep, xyr_RadUwFlux, xyr_RadDwFlux )
xyr_RadUwFlux = xyr_RadUwFlux * KDFWeight
xyr_RadDwFlux = xyr_RadDwFlux * KDFWeight
xyr_RadSUwFlux = xyr_RadSUwFlux + xyr_RadUwFlux
xyr_RadSDwFlux = xyr_RadSDwFlux + xyr_RadDwFlux
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 RadEarthSWV21Flux