subroutine RadMarsV1Flux( xy_SurfType, xy_SurfMajCompIce, xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_Temp, xyr_Temp, xy_SurfTemp, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )
! USE statements
!
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 物理・数学定数設定
! Physical and mathematical constants settings
!
use constants0, only: PI ! $ \pi $ .
! 円周率. Circular constant
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, 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, SetMarsDustCalcDOD067
! 散乱を無視した放射伝達方程式
! Radiative transfer equation without considering scattering
!
use rad_rte_nonscat, only : RadRTENonScatWrapper
!
! Solve radiative transfer equation in two stream approximation
!
use rad_rte_two_stream_app, only: RadRTETwoStreamAppHomogAtm, RadRTETwoStreamAppLW
! プランク関数の計算
! Calculate Planck function
!
use planck_func, only : Integ_PF_GQ_Array3D, Integ_PF_GQ_Array2D, Integ_DPFDT_GQ_Array2D
! ヒストリデータ出力
! 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 ) :: xyr_Temp (0:imax-1, 1:jmax, 0: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) :: 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) :: xyz_SSA(0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_AF (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_LWSurfAlbedo(0:imax-1, 1:jmax)
real(DP) :: xyr_IntPF (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xy_SurfIntPF (0:imax-1, 1:jmax)
real(DP) :: xy_SurfIntDPFDT(0:imax-1, 1:jmax)
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) :: xyz_DustDensScledOptDep(0:imax-1, 1:jmax, 1:kmax)
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 ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
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 )
!!$ call SetMarsDustCalcDOD067( &
!!$ & xyr_Press, xyz_QDust, & ! (in)
!!$ & xyr_DOD067 & ! (out)
!!$ & )
! 短波放射
! 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
!
! Surface albedo for long wave is set.
!
xy_LWSurfAlbedo = 0.0_DP
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.35_DP
AF = 0.36_DP
xyz_SSA = SSA
xyz_AF = AF
xyr_DOD = QeRat * xyr_DOD067
!----------
! Modification of dust optical depth for use in non-scattering calculation
!!$ xyr_DOD = ( 1.0_DP - SSA ) * xyr_DOD
!
!!$ 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, & ! (in )
!!$ & xyr_RadLUwFluxComp, xyr_RadLDwFluxComp, & ! (out)
!!$ & xyra_DelRadLUwFluxComp, xyra_DelRadLDwFluxComp, & ! (out)
!!$ & WNs, WNe, NumGaussNode & ! (in ) optional
!!$ & )
!----------
!!$ i = 0
!!$ j = jmax/2+1
!!$ do k = 0, kmax
!!$ write( 70, * ) k, &
!!$ & xyr_RadLUwFluxComp(i,j,k), &
!!$ & xyr_RadLDwFluxComp(i,j,k), &
!!$ & xyra_DelRadLUwFluxComp(i,j,k,0), &
!!$ & xyra_DelRadLUwFluxComp(i,j,k,1), &
!!$ & xyra_DelRadLDwFluxComp(i,j,k,0), &
!!$ & xyra_DelRadLDwFluxComp(i,j,k,1)
!!$ end do
!!$ call flush( 70 )
! Integrate Planck function and temperature derivative of it
!
call Integ_PF_GQ_Array3D( WNs, WNe, NumGaussNode, 0, imax-1, 1, jmax, 0, kmax, xyr_Temp, xyr_IntPF )
call Integ_PF_GQ_Array2D( WNs, WNe, NumGaussNode, 0, imax-1, 1, jmax, xy_SurfTemp, xy_SurfIntPF )
! call Integ_DPFDT_GQ_Array2D( &
! & WNs, WNe, NumGaussNode, & ! (in )
! & 0, imax-1, 1, jmax, xyz_Temp(:,:,1), & ! (in )
! & xy_IntDPFDT1 & ! (out)
! & )
call Integ_DPFDT_GQ_Array2D( WNs, WNe, NumGaussNode, 0, imax-1, 1, jmax, xy_SurfTemp, xy_SurfIntDPFDT )
!
xyr_IntPF = PI * xyr_IntPF
xy_SurfIntPF = xy_SurfEmis * PI * xy_SurfIntPF
! xy_IntDPFDT1 = PI * xy_IntDPFDT1
xy_SurfIntDPFDT = xy_SurfEmis * PI * xy_SurfIntDPFDT
!
call RadRTETwoStreamAppLW( xyz_SSA, xyz_AF, xyr_DOD, xy_LWSurfAlbedo, xyr_IntPF, xy_SurfIntPF, xy_SurfIntDPFDT, xyr_RadLUwFluxComp, xyr_RadLDwFluxComp, xyra_DelRadLUwFluxComp, xyra_DelRadLDwFluxComp )
!!$ i = 0
!!$ j = jmax/2+1
!!$ do k = 0, kmax
!!$ write( 80, * ) k, &
!!$ & xyr_RadLUwFluxComp(i,j,k), &
!!$ & xyr_RadLDwFluxComp(i,j,k), &
!!$ & xyra_DelRadLUwFluxComp(i,j,k,0), &
!!$ & xyra_DelRadLUwFluxComp(i,j,k,1), &
!!$ & xyra_DelRadLDwFluxComp(i,j,k,0), &
!!$ & xyra_DelRadLDwFluxComp(i,j,k,1)
!!$ end do
!!$ call flush( 80 )
!!$ stop
!----------
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, 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.55_DP
AF = 0.55_DP
xyz_SSA = SSA
xyz_AF = AF
xyr_DOD = QeRat * xyr_DOD067
!----------
! Modification of dust optical depth for use in non-scattering calculation
!!$ xyr_DOD = ( 1.0_DP - SSA ) * xyr_DOD
!
!!$ 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, & ! (in )
!!$ & xyr_RadLUwFluxComp, xyr_RadLDwFluxComp, & ! (out)
!!$ & xyra_DelRadLUwFluxComp, xyra_DelRadLDwFluxComp, & ! (out)
!!$ & WNs, WNe, NumGaussNode & ! (in ) optional
!!$ & )
!----------
! Integrate Planck function and temperature derivative of it
!
call Integ_PF_GQ_Array3D( WNs, WNe, NumGaussNode, 0, imax-1, 1, jmax, 0, kmax, xyr_Temp, xyr_IntPF )
call Integ_PF_GQ_Array2D( WNs, WNe, NumGaussNode, 0, imax-1, 1, jmax, xy_SurfTemp, xy_SurfIntPF )
! call Integ_DPFDT_GQ_Array2D( &
! & WNs, WNe, NumGaussNode, & ! (in )
! & 0, imax-1, 1, jmax, xyz_Temp(:,:,1), & ! (in )
! & xy_IntDPFDT1 & ! (out)
! & )
call Integ_DPFDT_GQ_Array2D( WNs, WNe, NumGaussNode, 0, imax-1, 1, jmax, xy_SurfTemp, xy_SurfIntDPFDT )
!
xyr_IntPF = PI * xyr_IntPF
xy_SurfIntPF = xy_SurfEmis * PI * xy_SurfIntPF
! xy_IntDPFDT1 = PI * xy_IntDPFDT1
xy_SurfIntDPFDT = xy_SurfEmis * PI * xy_SurfIntDPFDT
!
call RadRTETwoStreamAppLW( xyz_SSA, xyz_AF, xyr_DOD, xy_LWSurfAlbedo, xyr_IntPF, xy_SurfIntPF, xy_SurfIntDPFDT, xyr_RadLUwFluxComp, xyr_RadLDwFluxComp, xyra_DelRadLUwFluxComp, xyra_DelRadLDwFluxComp )
!----------
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 )
do k = 1, kmax
xyz_DustDensScledOptDep(:,:,k) = ( xyr_DOD067(:,:,k-1) - xyr_DOD067(:,:,k) ) / ( xyr_Press (:,:,k-1) - xyr_Press (:,:,k) ) * Grav
end do
call HistoryAutoPut( TimeN, 'DustDensScledOptDep', xyz_DustDensScledOptDep )
end subroutine RadMarsV1Flux