subroutine RadMarsV1Flux( xy_SurfType, xy_SurfMajCompIce, xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_Temp, xyz_VirTemp, xy_SurfTemp, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )
! USE statements
!
! 物理・数学定数設定
! Physical and mathematical constants settings
!
use constants0, only: PI ! $ \pi $ .
! 円周率. Circular constant
use constants, only: GasRDry
! 雪と海氷の定数の設定
! Setting constants of snow and sea ice
!
use constants_snowseaice, only: CO2IceThreshold, CO2IceEmisS, CO2IceEmisN
! 座標データ設定
! Axes data settings
!
use axesset, only: y_Lat
! 時刻管理
! Time control
!
use timeset, only: TimeN, DelTime, TimesetClockStart, TimesetClockStop
! 短波入射 (太陽入射)
! Short wave (insolation) incoming
!
use rad_short_income, only : RadShortIncome
use rad_Mars_15m, only : rad15m_main
use set_Mars_dust, only : SetMarsDustSetDOD067
! 散乱を無視した放射伝達方程式
! Radiative transfer equation without considering scattering
!
use rad_rte_nonscat, only : RadRTENonScatWrapper
use rad_rte_two_stream_app, only : RadRTETwoStreamAppHomogAtm
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
integer , intent(in ) :: xy_SurfType (0:imax-1, 1:jmax)
real(DP), intent(in ) :: xy_SurfMajCompIce(0:imax-1, 1:jmax)
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_VirTemp (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax)
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), intent(out) :: xyr_RadLUwFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadLDwFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP), intent(out) :: xyra_DelRadLDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP), parameter :: DiffFact = 1.66_DP
real(DP) :: xyz_Rho (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: PlanetLonFromVE
real(DP) :: xyr_DOD067 (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_DOD (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyz_DelTrans(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyrr_Trans (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
real(DP) :: QeRat
real(DP) :: SSA
real(DP) :: AF
real(DP) :: xyr_RadLUwFluxComp (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_RadLDwFluxComp (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyra_DelRadLUwFluxComp(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP) :: xyra_DelRadLDwFluxComp(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP) :: MajCompIceThreshold
real(DP) :: MajCompIceEmis
real(DP) :: xy_SurfEmis(0:imax-1, 1:jmax)
real(DP) :: xy_InAngle (0:imax-1, 1:jmax)
real(DP) :: DistFromStarScld
real(DP) :: DiurnalMeanFactor
real(DP) :: SolarFluxTOA
real(DP) :: WNs
real(DP) :: WNe
integer, parameter :: NumGaussNode = 5
integer :: i
integer :: j
integer :: k
integer :: kk
!!$ real(DP) :: MaxError
! 初期化
! Initialization
!
if ( .not. rad_Mars_V1_inited ) call RadMarsV1Init
xyz_Rho = xyz_Press / ( GasRDry * xyz_VirTemp )
call RadShortIncome( xy_InAngle = xy_InAngle, DistFromStarScld = DistFromStarScld, DiurnalMeanFactor = DiurnalMeanFactor, PlanetLonFromVE = PlanetLonFromVE )
SolarFluxTOA = SolarConst / DistFromStarScld**2 * DiurnalMeanFactor
PlanetLonFromVE = PlanetLonFromVE * 180.0_DP / PI
! Dust optical depth
!
call SetMarsDustSetDOD067( PlanetLonFromVE, xyr_Press, xyz_Press, xyr_DOD067 )
! 短波放射
! Short wave radiation
!
!!$ QeRat = 0.9837_DP ! Ockert-Bell et al. (1997)
!!$ xyz_SSA = 0.86_DP
!!$ xyz_AF = 0.68_DP
!!$ QeRat = 1.0_DP ! Clancy and Lee (1991)
SSA = 0.92_DP
AF = 0.55_DP
call RadRTETwoStreamAppHomogAtm( xy_SurfAlbedo, SolarFluxTOA, xy_InAngle, SSA, AF, xyr_DOD067, xyr_RadSUwFlux, xyr_RadSDwFlux )
!!$ MaxError = 0.0_DP
!!$ do k = 0, kmax
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ MaxError = max( MaxError, &
!!$ & abs( OLD_xyr_RadSFlux(i,j,k) - xyr_RadSFlux(i,j,k) ) )
!!$ MaxError = max( MaxError, &
!!$ & abs( OLD_xyr_RadSFlux(i,j,k) - ( xyr_RadSUwFlux(i,j,k) - xyr_RadSDwFlux(i,j,k) ) ) )
!!$ end do
!!$ end do
!!$ end do
!!$ write( 6, * ) MaxError
!!$ write( 6, * ) MaxError, xyr_RadSUwFlux(0,jmax/2+1,kmax)
!!$ MaxError = 0.0_DP
!!$ do k = 0, kmax
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ MaxError = max( MaxError, &
!!$ & abs( OLD_xyr_RadLFlux(i,j,k) - xyr_RadLFlux(i,j,k) ) )
!!$ MaxError = max( MaxError, &
!!$ & abs( OLD_xyra_DelRadLFlux(i,j,k,0) - xyra_DelRadLFlux(i,j,k,0) ) )
!!$ MaxError = max( MaxError, &
!!$ & abs( OLD_xyra_DelRadLFlux(i,j,k,1) - xyra_DelRadLFlux(i,j,k,1) ) )
!!$ end do
!!$ end do
!!$ end do
!!$ write( 6, * ) MaxError
!!$ !
!!$ xyr_RadSFlux = OLD_xyr_RadSFlux
!!$ xyr_RadLFlux = OLD_xyr_RadLFlux
!!$ xyra_DelRadLFlux = OLD_xyra_DelRadLFlux
!!$ xyr_RadSUwFlux = OLD_xyr_RadSFlux
!!$ xyr_RadSDwFlux = 0.0_DP
! 長波放射
! Long wave radiation
!
xyr_RadLUwFlux = 0.0_DP
xyr_RadLDwFlux = 0.0_DP
xyra_DelRadLUwFlux = 0.0_DP
xyra_DelRadLDwFlux = 0.0_DP
! Surface emissivity
!
xy_SurfEmis = 1.0_DP
MajCompIceThreshold = CO2IceThreshold
do j = 1, jmax
if ( y_Lat(j) < 0.0_DP ) then
MajCompIceEmis = CO2IceEmisS
else
MajCompIceEmis = CO2IceEmisN
end if
do i = 0, imax-1
if ( xy_SurfType(i,j) > 0 ) then
if ( xy_SurfMajCompIce(i,j) > MajCompIceThreshold ) then
xy_SurfEmis(i,j) = MajCompIceEmis
else if ( xy_SurfMajCompIce(i,j) < 0.0_DP ) then
xy_SurfEmis(i,j) = xy_SurfEmis(i,j)
else
xy_SurfEmis(i,j) = ( MajCompIceEmis - xy_SurfEmis(i,j) ) / ( MajCompIceThreshold - 0.0_DP ) * ( xy_SurfMajCompIce(i,j) - 0.0_DP ) + xy_SurfEmis(i,j)
end if
end if
end do
end do
! Flux from 0 to 500 cm-1
!
WNs = 0.0d2
WNe = 500.0d2
QeRat = 0.17_DP ! Wavenumber averaged extinction coefficient
! ssa = 0.35d0
! af = 0.36d0
xyr_DOD = QeRat * xyr_DOD067
do k = 1, kmax
xyz_DelTrans(:,:,k) = exp( - DiffFact * ( xyr_DOD(:,:,k-1) - xyr_DOD(:,:,k) ) )
end do
!
do k = 0, kmax
do kk = k, k
xyrr_Trans(:,:,k,kk) = 1.0_DP
end do
do kk = k+1, kmax
xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk-1) * xyz_DelTrans(:,:,kk)
end do
end do
do k = 0, kmax
do kk = 0, k-1
xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,kk,k)
end do
end do
call RadRTENonScatWrapper( xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans, xyr_RadLUwFluxComp, xyr_RadLDwFluxComp, xyra_DelRadLUwFluxComp, xyra_DelRadLDwFluxComp, WNs, WNe, NumGaussNode )
xyr_RadLUwFlux = xyr_RadLUwFlux + xyr_RadLUwFluxComp
xyr_RadLDwFlux = xyr_RadLDwFlux + xyr_RadLDwFluxComp
xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadLUwFluxComp
xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadLDwFluxComp
! Flux from 500 to 850 cm-1
!
QeRat = 0.25_DP ! Wavenumber averaged extinction coefficient
SSA = 0.45_DP ! Wavenumber averaged single scattering albedo
call rad15m_main( TimeN, DelTime, xyz_Temp, xyr_Press, xyz_Press, xy_SurfTemp, xyz_Rho, xyr_DOD067, QeRat, SSA, xy_SurfEmis, xyr_RadLUwFluxComp, xyra_DelRadLUwFluxComp )
xyr_RadLDwFluxComp = 0.0_DP
xyra_DelRadLDwFluxComp = 0.0_DP
xyr_RadLUwFlux = xyr_RadLUwFlux + xyr_RadLUwFluxComp
xyr_RadLDwFlux = xyr_RadLDwFlux + xyr_RadLDwFluxComp
xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadLUwFluxComp
xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadLDwFluxComp
! Flux from 850 to 2000 cm-1
!
WNs = 850.0d2
WNe = 2000.0d2
QeRat = 0.41_DP ! Wavenumber averaged extinction coefficient
! ssa = 0.55d0
! af = 0.55d0
xyr_DOD = QeRat * xyr_DOD067
do k = 1, kmax
xyz_DelTrans(:,:,k) = exp( - DiffFact * ( xyr_DOD(:,:,k-1) - xyr_DOD(:,:,k) ) )
end do
!
do k = 0, kmax
do kk = k, k
xyrr_Trans(:,:,k,kk) = 1.0_DP
end do
do kk = k+1, kmax
xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk-1) * xyz_DelTrans(:,:,kk)
end do
end do
do k = 0, kmax
do kk = 0, k-1
xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,kk,k)
end do
end do
call RadRTENonScatWrapper( xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans, xyr_RadLUwFluxComp, xyr_RadLDwFluxComp, xyra_DelRadLUwFluxComp, xyra_DelRadLDwFluxComp, WNs, WNe, NumGaussNode )
xyr_RadLUwFlux = xyr_RadLUwFlux + xyr_RadLUwFluxComp
xyr_RadLDwFlux = xyr_RadLDwFlux + xyr_RadLDwFluxComp
xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadLUwFluxComp
xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadLDwFluxComp
! Output variables
!
call HistoryAutoPut( TimeN, 'DOD067', xyr_DOD067(:,:,0) )
end subroutine RadMarsV1Flux