| Class | radiation_RL78 |
| In: |
radiation/radiation_RL78.f90
|
Note that Japanese and English are described in parallel.
長波放射モデル.
This is a model of long wave radiation.
Roewe, D., and K.-N. Liou, Influence of cirrus clouds on the infrared cooling rate in the troposphere and lower stratosphere, J. Appl. Met., 17, 92-106, 1978.
| !$ ! 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) |
| Subroutine : | |
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_QVap(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_QO3(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_Press(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) |
| xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(out) |
subroutine RadiationRL78Flux( xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, xy_SurfTemp, xyr_RadLFlux, xyra_DelRadLFlux )
! USE statements
!
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
!
! Physical constants settings
!
use constants, only: Grav, PI ! $ \pi $ .
! Circular constant
! 時刻管理
! Time control
!
use timeset, only: TimeN ! ステップ $ t $ の時刻.
! Time of step $ t $.
! 日付および時刻の取り扱い
! Date and time handler
!
use dc_date, only: operator(-), operator(>=), operator(+), operator(==), toChar
use read_time_series, only: SetValuesFromTimeSeriesWrapper
use set_cloud, only : SetCloudLW
real(DP), intent(in ):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_QVap (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_QO3 (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_Press (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xy_SurfTemp (0:imax-1, 1:jmax)
real(DP), intent(out):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out):: xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
!
! Work variables
!
real(DP):: RefPress
real(DP):: xyz_QCO2 (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_H2ODelAbsAmt (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_H2ODelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_O3DelScaledAbsAmt (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyrr_TransCO2 (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
real(DP):: xyz_CloudDelOptDep (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_TransCloudOneLayer(0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyrr_TransCloud (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
real(DP):: xyrr_Trans (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
real(DP):: xyz_PlankFunc (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xy_SurfPlankFunc (0:imax-1, 1:jmax)
real(DP):: xy_DPFDT0 (0:imax-1, 1:jmax)
real(DP):: xy_DPFDT1 (0:imax-1, 1:jmax)
logical :: flag_dry_atmosphere
integer :: k, kk, k1, k2
integer :: n
!!$ integer :: i, j
if ( sw_fs ) then
call RadiationRoeweLiouInit
end if
!!$ xyz_QCO2(:,:,:) = 382.0d-6 * CO2MolWeight / MeanMolWeight
xyz_QCO2(:,:,:) = VMRCO2 * CO2MolWeight / MeanMolWeight
!!$ xyz_QO3 (:,:,:) = 0.0d0
!!$ xyz_QVap(:,:,:) = 0.0d0
!!$ xyz_QCO2(:,:,:) = 300.0d-6 * CO2MolWeight / MeanMolWeight
!!$ xyz_QCO2(:,:,:) = 0.0d0
!!$ call SetRefAtmPro( &
!!$ & 5, .true., xyz_Press, &
!!$ & xyz_QO3 &
!!$ & )
!!$ if ( .true. ) then
if ( ( TimeN - PrevTimeSave >= IntTimeSave ) .or. ( .not. FlagTransSaved ) ) then
!!$ write( 6, * ) 'CalcTrans'
if ( .not. FlagTransSaved ) then
PrevTimeSave = TimeN
else
PrevTimeSave = PrevTimeSave + IntTimeSave
end if
FlagTransSaved = .true.
! Check for dry atmosphere
!
if ( all( xyz_QVap <= 0.0d0 ) ) then
flag_dry_atmosphere = .true.
!!$ write( 6, * ) 'Dry atmosphere'
else
flag_dry_atmosphere = .false.
end if
!!$ write( 6, * ) 'The validity of scaling absorber amount has to be checked.'
do k = 1, kmax
xyz_H2ODelAbsAmt(:,:,k) = xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
end do
RefPress = 1.0d5
do k = 1, kmax
xyz_H2ODelScaledAbsAmt(:,:,k) = ( xyz_Press(:,:,k) / RefPress )**H2OScaleIndex * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
!!$ xyz_H2ODelScaledAbsAmt(:,:,k) = &
!!$ & ( xyz_Press(:,:,k) / RefPress ) &
!!$ & * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
!!$ xyz_H2ODelScaledAbsAmt(:,:,k) = &
!!$ & ( xyz_Press(:,:,k) / RefPress ) &
!!$ & * xyz_H2ODelAbsAmt(:,:,k)
end do
do k = 1, kmax
xyz_CO2DelScaledAbsAmt(:,:,k) = ( xyz_Press(:,:,k) / RefPress )**CO2ScaleIndex * xyz_QCO2(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
!!$ xyz_CO2DelScaledAbsAmt(:,:,k) = &
!!$ & ( xyz_Press(:,:,k) / RefPress ) &
!!$ & * xyz_QCO2(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
end do
do k = 1, kmax
xyz_O3DelScaledAbsAmt(:,:,k) = ( xyz_Press(:,:,k) / RefPress )**O3ScaleIndex * xyz_QO3(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
!!$ xyz_O3DelScaledAbsAmt(:,:,k) = &
!!$ & xyz_QO3(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
end do
!-----
! Transmission in CO2 15 micron band
call RadiationRoeweLiouCalcCO2Trans( xyz_CO2DelScaledAbsAmt, xyrr_TransCO2 )
!-----
! Transmission by cloud
call SetCloudLW( xyz_CloudDelOptDep )
if ( all( xyz_CloudDelOptDep <= 0.0d0 ) ) then
xyrr_TransCloud(:,:,:,:) = 1.0d0
else
xyz_TransCloudOneLayer = exp( - xyz_CloudDelOptDep * DiffFactor )
do k1 = 0, kmax
k2 = k1
xyrr_TransCloud(:,:,k1,k2) = 1.0d0
do k2 = k1+1, kmax
xyrr_TransCloud(:,:,k1,k2) = xyrr_TransCloud(:,:,k1,k2-1) * xyz_TransCloudOneLayer(:,:,k2)
end do
end do
do k1 = 0, kmax
do k2 = 0, k1-1
xyrr_TransCloud(:,:,k1,k2) = xyrr_TransCloud(:,:,k2,k1)
end do
end do
end if
!-----
do n = 1, nbmax
if ( n <= 5 ) then
! H2O band
if ( flag_dry_atmosphere ) then
xyrr_Trans = 1.0d0
else
call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans )
end if
else if ( n <= 8 ) then
! H2O band & CO2 band
if ( flag_dry_atmosphere ) then
xyrr_Trans = 1.0d0
else
call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans )
end if
xyrr_Trans = xyrr_Trans * xyrr_TransCO2
else if ( n == 9 ) then
! H2O band
if ( flag_dry_atmosphere ) then
xyrr_Trans = 1.0d0
else
call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans )
end if
else if ( n == 10 ) then
! H2O band & continuum
if ( flag_dry_atmosphere ) then
xyrr_Trans = 1.0d0
else
call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans )
call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
end if
else if ( n == 11 ) then
! H2O continuum
if ( flag_dry_atmosphere ) then
xyrr_Trans = 1.0d0
else
xyrr_Trans = 1.0d0
call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
end if
else if ( n <= 21 ) then
! O3 band & H2O continuum
call RadiationRoeweLiouCalcBandTrans( n, xyz_O3DelScaledAbsAmt, xyrr_Trans )
if ( flag_dry_atmosphere ) then
xyrr_Trans = xyrr_Trans
else
call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
end if
else if ( n == 22 ) then
! H2O continuum
if ( flag_dry_atmosphere ) then
xyrr_Trans = 1.0d0
else
xyrr_Trans = 1.0d0
call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
end if
else if ( n <= 31 ) then
! H2O band
if ( flag_dry_atmosphere ) then
xyrr_Trans = 1.0d0
else
call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans )
end if
else if ( n == 32 ) then
xyrr_Trans = 1.0d0
else
write( 6, * ) 'Unexpected number of band, ', n, '.'
stop
end if
xyrr_Trans = xyrr_Trans * xyrr_TransCloud
xyrra_TransSave(:,:,:,:,n) = xyrr_Trans
end do
!!$ else
!!$
!!$ write( 6, * ) 'NO CalcTrans'
end if
xyr_RadLFlux (:,:,:) = 0.0d0
xyra_DelRadLFlux(:,:,:,:) = 0.0d0
do n = 1, nbmax
xyrr_Trans = xyrra_TransSave(:,:,:,:,n)
if ( all( xyrr_Trans(:,:,0,kmax) == 1.0d0 ) ) then
!!$ call CalcIntegratedPF2D( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, &
!!$ imax, jmax, xy_SurfTemp, xy_SurfPlankFunc )
call CalcIntegratedPFWithTable2D( n, imax, jmax, xy_SurfTemp, xy_SurfPlankFunc )
do k = 0, kmax
xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) + PI * xy_SurfPlankFunc(:,:) ! * xyrr_Trans(:,:,0,k)
end do
!!$ call Integ_DPFDT_GQ_Array2D( &
!!$ & a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, xy_SurfTemp, & ! (in )
!!$ & xy_DPFDT0 & ! (out)
!!$ & )
call CalcIntegratedPFWithTable2D( n, imax, jmax, xy_SurfTemp, xy_DPFDT0, .true. )
do k = 0, kmax
xyra_DelRadLFlux(:,:,k,0) = xyra_DelRadLFlux(:,:,k,0) + PI * xy_DPFDT0(:,:) ! * xyrr_Trans(:,:,0,k)
end do
else
!!$ call CalcIntegratedPF2D( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, &
!!$ imax, jmax, xy_SurfTemp, xy_SurfPlankFunc )
!!$ call CalcIntegratedPF ( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, &
!!$ imax, jmax, kmax, xyz_Temp, xyz_PlankFunc )
call CalcIntegratedPFWithTable2D( n, imax, jmax, xy_SurfTemp, xy_SurfPlankFunc )
call CalcIntegratedPFWithTable3D( n, imax, jmax, kmax, xyz_Temp, xyz_PlankFunc )
do k = 0, kmax
xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) + PI * xy_SurfPlankFunc(:,:) * xyrr_Trans(:,:,0,k)
do kk = 1, kmax
xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k) - PI * xyz_PlankFunc(:,:,kk) * ( xyrr_Trans(:,:,k,kk-1) - xyrr_Trans(:,:,k,kk) )
end do
end do
!!$ call Integ_DPFDT_GQ_Array2D( &
!!$ & a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, xy_SurfTemp, & ! (in )
!!$ & xy_DPFDT0 & ! (out)
!!$ & )
!!$ call Integ_DPFDT_GQ_Array2D( &
!!$ & a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, xyz_Temp(:,:,1), & ! (in )
!!$ & xy_DPFDT1 & ! (out)
!!$ & )
call CalcIntegratedPFWithTable2D( n, imax, jmax, xy_SurfTemp, xy_DPFDT0, .true. )
call CalcIntegratedPFWithTable2D( n, imax, jmax, xyz_Temp(:,:,1), xy_DPFDT1, .true. )
do k = 0, kmax
xyra_DelRadLFlux(:,:,k,0) = xyra_DelRadLFlux(:,:,k,0) + PI * xy_DPFDT0(:,:) * xyrr_Trans(:,:,0,k)
xyra_DelRadLFlux(:,:,k,1) = xyra_DelRadLFlux(:,:,k,1) - PI * xy_DPFDT1(:,:) * ( xyrr_Trans(:,:,k,0) - xyrr_Trans(:,:,k,1) )
end do
end if
end do
end subroutine RadiationRL78Flux
| Subroutine : | |
| wn1 : | real(DP), intent(in ) |
| wn2 : | real(DP), intent(in ) |
| num : | integer , intent(in ) |
| im : | integer , intent(in ) |
| jm : | integer , intent(in ) |
| km : | integer , intent(in ) |
| temp( im, jm, km ) : | real(DP), intent(in ) |
| pfinted( im, jm, km ) : | real(DP), intent(out) |
subroutine CalcIntegratedPF( wn1, wn2, num, im, jm, km, temp, pfinted )
real(DP), intent(in ) :: wn1,wn2
integer , intent(in ) :: num
integer , intent(in ) :: im, jm, km
real(DP), intent(in ) :: temp ( im, jm, km )
real(DP), intent(out) :: pfinted( im, jm, km )
!
! local variables
!
real(DP):: x( num ), w( num )
integer :: i
integer :: j
integer :: k
integer :: l
call gauleg( wn1, wn2, num, x, w )
pfinted(:,:,:) = 0.0d0
do l = 1, num
do k = 1, km
do j = 1, jm
do i = 1, im
pfinted(i,j,k) = pfinted(i,j,k) + pf( x( l ), temp(i,j,k) ) * w( l )
end do
end do
end do
end do
end subroutine CalcIntegratedPF
| Subroutine : | |
| wn1 : | real(DP), intent(in ) |
| wn2 : | real(DP), intent(in ) |
| num : | integer , intent(in ) |
| im : | integer , intent(in ) |
| jm : | integer , intent(in ) |
| temp( im, jm ) : | real(DP), intent(in ) |
| pfinted( im, jm ) : | real(DP), intent(out) |
subroutine CalcIntegratedPF2D( wn1, wn2, num, im, jm, temp, pfinted )
real(DP), intent(in ) :: wn1,wn2
integer , intent(in ) :: num
integer , intent(in ) :: im, jm
real(DP), intent(in ) :: temp ( im, jm )
real(DP), intent(out) :: pfinted( im, jm )
!
! local variables
!
real(dp) :: temp3d( im, jm, 1 ), pfinted3d( im, jm, 1 )
temp3d(:,:,1) = temp(:,:)
call CalcIntegratedPF( wn1, wn2, num, im, jm, 1, temp3d, pfinted3d )
pfinted(:,:) = pfinted3d(:,:,1)
end subroutine CalcIntegratedPF2D
| Subroutine : | |
| iband : | integer , intent(in ) |
| im : | integer , intent(in ) |
| jm : | integer , intent(in ) |
| xy_temp(im, jm) : | real(DP), intent(in ) |
| xy_IntegPF(im, jm) : | real(DP), intent(out) |
| flag_DPFDT : | logical , intent(in ), optional |
subroutine CalcIntegratedPFWithTable2D( iband, im, jm, xy_temp, xy_IntegPF, flag_DPFDT )
! USE statements
!
integer , intent(in ) :: iband
integer , intent(in ) :: im, jm
real(DP), intent(in ) :: xy_temp (im, jm)
real(DP), intent(out) :: xy_IntegPF(im, jm)
logical , intent(in ), optional :: flag_DPFDT
!
! local variables
!
real(DP) :: xyz_Temp (im, jm, 1)
real(DP) :: xyz_IntegPF(im, jm, 1)
xyz_Temp(:,:,1) = xy_Temp(:,:)
call CalcIntegratedPFWithTable3D( iband, im, jm, 1, xyz_temp, xyz_IntegPF, flag_DPFDT )
xy_IntegPF(:,:) = xyz_IntegPF(:,:,1)
end subroutine CalcIntegratedPFWithTable2D
| Subroutine : | |
| iband : | integer , intent(in ) |
| im : | integer , intent(in ) |
| jm : | integer , intent(in ) |
| km : | integer , intent(in ) |
| xyz_temp(im, jm, km) : | real(DP), intent(in ) |
| xyz_IntegPF(im, jm, km) : | real(DP), intent(out) |
| flag_DPFDT : | logical , intent(in ), optional |
subroutine CalcIntegratedPFWithTable3D( iband, im, jm, km, xyz_temp, xyz_IntegPF, flag_DPFDT )
! USE statements
!
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
integer , intent(in ) :: iband
integer , intent(in ) :: im, jm, km
real(DP), intent(in ) :: xyz_temp (im, jm, km)
real(DP), intent(out) :: xyz_IntegPF(im, jm, km)
logical , intent(in ), optional :: flag_DPFDT
!
! local variables
!
logical , save :: flag_table_inited
real(DP), save :: TableTempMin
real(DP), save :: TableTempMax
real(DP), save :: TableTempIncrement
!!$ integer , save :: ntmax
!!$ real(DP), save, allocatable :: a_TableTemp (:)
real(DP), save, allocatable :: aa_TableIPF (:,:)
real(DP), save, allocatable :: aa_TableIDPFDT(:,:)
real(DP) , allocatable :: GQP(:)
real(DP) , allocatable :: GQW(:)
logical :: local_flag_DPFDT
integer :: xyz_TempIndex(im, jm, km)
integer :: i
integer :: j
integer :: k
integer :: l
integer :: m
integer :: n
data flag_table_inited /.false./
if ( .not. flag_table_inited ) then
TableTempMin = 50.0d0
TableTempMax = 400.0d0
TableTempIncrement = 0.1d0
ntmax = ( TableTempMax - TableTempMin ) / TableTempIncrement + 1
allocate( a_TableTemp (1:ntmax) )
allocate( aa_TableIPF (1:ntmax, 1:nbmax) )
allocate( aa_TableIDPFDT(1:ntmax, 1:nbmax) )
do m = 1, ntmax
a_TableTemp(m) = TableTempMin + TableTempIncrement * ( m - 1 )
end do
aa_TableIPF (:,:) = 0.0d0
aa_TableIDPFDT(:,:) = 0.0d0
allocate( GQP(1:NGaussQuad) )
allocate( GQW(1:NGaussQuad) )
do n = 1, nbmax
call gauleg( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, GQP, GQW )
do m = 1, ntmax
do l = 1, NGaussQuad
aa_TableIPF (m,n) = aa_TableIPF (m,n) + pf ( GQP(l), a_TableTemp(m) ) * GQW(l)
aa_TableIDPFDT(m,n) = aa_TableIDPFDT(m,n) + DPFDT( GQP(l), a_TableTemp(m) ) * GQW(l)
end do
end do
end do
deallocate( GQP )
deallocate( GQW )
flag_table_inited = .true.
end if
do k = 1, km
do j = 1, jm
do i = 1, im
if ( ( xyz_Temp(i,j,k) < a_TableTemp(1) ) .or. ( xyz_Temp(i,j,k) > a_TableTemp(ntmax) ) ) then
call MessageNotify( 'E', module_name, 'Temperature is not appropriate, Temp(%d,%d,%d) = %f.', i = (/i, j, k/), d = (/xyz_Temp(i,j,k)/) )
end if
xyz_TempIndex(i,j,k) = int( ( xyz_Temp(i,j,k) - TableTempMin ) / TableTempIncrement ) + 2
if ( xyz_TempIndex(i,j,k) == 1 ) then
xyz_TempIndex(i,j,k) = 2
else if ( xyz_TempIndex(i,j,k) >= ntmax ) then
xyz_TempIndex(i,j,k) = ntmax - 1
end if
!!$ xyz_TempIndex(i,j,k) = ntmax-1
!!$ search_index: do m = 2, ntmax-1
!!$ if ( a_TableTemp(m) >= xyz_Temp(i,j,k) ) then
!!$ xyz_TempIndex(i,j,k) = m
!!$ exit search_index
!!$ end if
!!$ end do search_index
end do
end do
end do
local_flag_DPFDT = .false.
if ( present( flag_DPFDT ) ) then
if ( flag_DPFDT ) then
local_flag_DPFDT = .true.
end if
end if
if ( .not. local_flag_DPFDT ) then
do k = 1, km
do j = 1, jm
do i = 1, im
m = xyz_TempIndex(i,j,k)
!!$ xyz_IntegPF(i,j,k) = &
!!$ & ( aa_TableIPF( m, iband ) - aa_TableIPF( m-1, iband ) ) &
!!$ & / ( a_TableTemp( m ) - a_TableTemp( m-1 ) ) &
!!$ & * ( xyz_Temp(i,j,k) - a_TableTemp( m-1 ) ) &
!!$ & + aa_TableIPF( m-1, iband )
xyz_IntegPF(i,j,k) = aa_TableIPF(m-1,iband) * ( xyz_Temp (i,j,k) - a_TableTemp( m ) ) * ( xyz_Temp (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m-1 ) - a_TableTemp( m ) ) * ( a_TableTemp( m-1 ) - a_TableTemp( m+1 ) ) ) + aa_TableIPF(m ,iband) * ( xyz_Temp (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m ) - a_TableTemp( m+1 ) ) ) + aa_TableIPF(m+1,iband) * ( xyz_Temp (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp (i,j,k) - a_TableTemp( m ) ) / ( ( a_TableTemp( m+1 ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m+1 ) - a_TableTemp( m ) ) )
end do
end do
end do
else
do k = 1, km
do j = 1, jm
do i = 1, im
m = xyz_TempIndex(i,j,k)
!!$ xyz_IntegPF(i,j,k) = &
!!$ & ( aa_TableIDPFDT( m, iband ) - aa_TableIDPFDT( m-1, iband ) ) &
!!$ & / ( a_TableTemp ( m ) - a_TableTemp ( m-1 ) ) &
!!$ & * ( xyz_Temp(i,j,k) - a_TableTemp( m-1 ) ) &
!!$ & + aa_TableIDPFDT( m-1, iband )
xyz_IntegPF(i,j,k) = aa_TableIDPFDT(m-1,iband) * ( xyz_Temp (i,j,k) - a_TableTemp( m ) ) * ( xyz_Temp (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m-1 ) - a_TableTemp( m ) ) * ( a_TableTemp( m-1 ) - a_TableTemp( m+1 ) ) ) + aa_TableIDPFDT(m ,iband) * ( xyz_Temp (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m ) - a_TableTemp( m+1 ) ) ) + aa_TableIDPFDT(m+1,iband) * ( xyz_Temp (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp (i,j,k) - a_TableTemp( m ) ) / ( ( a_TableTemp( m+1 ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m+1 ) - a_TableTemp( m ) ) )
end do
end do
end do
end if
end subroutine CalcIntegratedPFWithTable3D
| Function : | |
| DPFDT : | real(DP) |
| WN : | real(DP), intent(in ) |
| Temp : | real(DP), intent(in ) |
real(DP) function DPFDT( WN, Temp )
! USE statements
!
real(DP), intent(in ) :: WN
real(DP), intent(in ) :: Temp
real(DP), parameter :: c = 2.99792458d8 , planc = 6.6260755d-34 , boltz = 1.380658d-23
real(DP) :: ExpTerm
real(DP) :: PF
ExpTerm = exp( planc * c * ( WN + 1.0d-10 ) / ( boltz * Temp ) )
PF = 2.0d0 * planc * c * c * WN * WN * WN / ( ExpTerm - 1.0d0 )
DPFDT = 1.0d0 / ( 2.0d0 * c * WN * WN * boltz ) * ( pf / Temp )**2 * ExpTerm
end function DPFDT
| Subroutine : | |
| WN : | real(DP), intent(in ) |
| xy_Temp(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| xy_DPFDT(0:imax-1, 1:jmax) : | real(DP), intent(out) |
subroutine DPFDT_Array2D( WN, xy_Temp, xy_DPFDT )
! USE statements
!
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
real(DP), intent(in ) :: WN
real(DP), intent(in ) :: xy_Temp (0:imax-1, 1:jmax)
real(DP), intent(out) :: xy_DPFDT(0:imax-1, 1:jmax)
real(DP), parameter :: c = 2.99792458d8 , planc = 6.6260755d-34 , boltz = 1.380658d-23
real(DP) :: xy_ExpTerm(0:imax-1, 1:jmax)
real(DP) :: xy_PF (0:imax-1, 1:jmax)
xy_ExpTerm(:,:) = exp( planc * c * ( WN + 1.0d-10 ) / ( boltz * xy_Temp ) )
xy_PF(:,:) = 2.0d0 * planc * c * c * WN * WN * WN / ( xy_ExpTerm(:,:) - 1.0d0 )
xy_DPFDT(:,:) = 1.0d0 / ( 2.0d0 * c * WN * WN * boltz ) * ( xy_pf(:,:) / xy_Temp(:,:) )**2 * xy_ExpTerm(:,:)
end subroutine DPFDT_ARRAY2D
| Variable : | |||
| IntTimeSave : | type(DC_DIFFTIME), save
|
| Subroutine : | |
| WN1 : | real(DP), intent(in ) |
| WN2 : | real(DP), intent(in ) |
| Num : | integer , intent(in ) |
| xy_Temp(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| xy_DPFDTInted(0:imax-1, 1:jmax) : | real(DP), intent(out) |
subroutine Integ_DPFDT_GQ_Array2D( WN1, WN2, Num, xy_Temp, xy_DPFDTInted )
! USE statements
!
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
real(DP), intent(in ) :: WN1
real(DP), intent(in ) :: WN2
integer , intent(in ) :: Num
real(DP), intent(in ) :: xy_Temp (0:imax-1, 1:jmax)
real(DP), intent(out) :: xy_DPFDTInted(0:imax-1, 1:jmax)
!
! local variables
!
real(DP):: GP( Num )
real(DP):: GW( Num )
real(DP):: xy_DPFDT(0:imax-1, 1:jmax)
integer :: l
call gauleg( WN1, WN2, Num, GP, GW )
xy_DPFDTInted(:,:) = 0.0d0
do l = 1, num
call DPFDT_Array2D( GP(l), xy_Temp, xy_DPFDT )
xy_DPFDTInted(:,:) = xy_DPFDTInted(:,:) + xy_DPFDT(:,:) * GW(l)
end do
end subroutine Integ_DPFDT_GQ_Array2D
| Variable : | |||
| PrevTimeSave : | type(DC_DIFFTIME), save
|
| Subroutine : | |
| iband : | integer , intent(in ) |
| xyz_Temp(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_QVap(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_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) : | real(DP), intent(inout) |
subroutine RadiationRLCalcH2OContTrans( iband, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
! USE statements
!
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
integer , intent(in ):: iband
real(DP), intent(in ):: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_QVap (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_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(inout):: xyrr_Trans (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!
! Work variables
!
real(DP):: H2OContConstant
real(DP):: a_H2OContParam(3)
real(DP):: RefTemp
real(DP):: gamma
real(DP):: xyz_PressH2O (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyz_DelOptDep (0:imax-1, 1:jmax, 1:kmax)
real(DP):: xyrr_H2OContOptDep(0:imax-1, 1:jmax, 0:kmax, 0:kmax)
integer :: k1
integer :: k2
integer :: k3
integer :: n
if ( ( iband < 9 ) .and. ( iband > 11 ) ) then
write( 6, * ) 'Unexpected number of band in RadiationRLCalcH2OContTrans.'
stop
end if
! These constants are obtained from Roberts et al. (1976) [Applied Optics]
a_H2OContParam(1) = 1.25d-22 ! mol-1 cm2 atm-1
a_H2OContParam(2) = 1.67d-19 ! mol-1 cm2 atm-1
a_H2OContParam(3) = 7.87d-3 ! (cm)
! unit conversion
a_H2OContParam(1) = a_H2OContParam(1) / ( H2OMolWeight / AvogNum ) * 1.0d-4 / 101325.0d0 ! (kg-1 m2 Pa-1)
a_H2OContParam(2) = a_H2OContParam(2) / ( H2OMolWeight / AvogNum ) * 1.0d-4 / 101325.0d0 ! (kg-1 m2 Pa-1)
a_H2OContParam(3) = a_H2OContParam(3) * 1.0d-2 ! (m)
n = iband
! This constant is deribed from Roberts et al. (1976) [Applied Optics]
!!$ gamma = 0.005d0
gamma = 0.0d0
xyz_PressH2O = xyz_Press * xyz_QVap * MeanMolWeight / H2OMolWeight
RefTemp = 1800.0d0
H2OContConstant = a_H2OContParam(1) + a_H2OContParam(2) * exp( - a_H2OContParam(3) * ( a_BandParam(1,n) + a_BandParam(2,n) ) * 0.5d0 )
xyz_DelOptDep = H2OContConstant * exp( RefTemp * ( 1.0d0 / xyz_Temp - 1.0d0 / 296.0d0 ) ) * ( xyz_PressH2O + gamma * ( xyz_Press - xyz_PressH2O ) ) * xyz_H2ODelAbsAmt
xyrr_H2OContOptDep(:,:,:,:) = 0.0d0
do k1 = 0, kmax
do k2 = k1, kmax
do k3 = k1+1, k2
xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) + xyz_DelOptDep(:,:,k3)
end do
xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) * DiffFactor
xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k1,k2) * exp( - xyrr_H2OContOptDep(:,:,k1,k2) )
end do
end do
do k1 = 0, kmax
do k2 = 0, k1-1
xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
end do
end do
do k1 = 0, kmax
k2 = k1
xyrr_Trans(:,:,k1,k2) = 1.0d0
end do
end subroutine RadiationRLCalcH2OContTrans
| Subroutine : | |
| iband : | integer , intent(in ) |
| xyz_DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) : | real(DP), intent(out) |
subroutine RadiationRoeweLiouCalcBandTrans( iband, xyz_DelScaledAbsAmt, xyrr_Trans )
! USE statements
!
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
integer , intent(in ):: iband
real(DP), intent(in ):: xyz_DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(out):: xyrr_Trans (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!
! Work variables
!
real(DP):: xyrr_AbsAmt (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
integer :: k1
integer :: k2
integer :: k3
integer :: n
n = iband
xyrr_AbsAmt(:,:,:,:) = 0.0d0
do k1 = 0, kmax
do k2 = k1, kmax
do k3 = k1+1, k2
xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) + xyz_DelScaledAbsAmt(:,:,k3)
end do
xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) * DiffFactor
xyrr_Trans(:,:,k1,k2) = - a_BandParam(3,n) * xyrr_AbsAmt(:,:,k1,k2) / sqrt( 1.0d0 + a_BandParam(3,n) / a_BandParam(4,n) * xyrr_AbsAmt(:,:,k1,k2) )
xyrr_Trans(:,:,k1,k2) = exp( xyrr_Trans(:,:,k1,k2) )
end do
end do
do k1 = 0, kmax
do k2 = 0, k1-1
xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
end do
end do
do k1 = 0, kmax
k2 = k1
xyrr_Trans(:,:,k1,k2) = 1.0d0
end do
end subroutine RadiationRoeweLiouCalcBandTrans
| Subroutine : | |
| xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) : | real(DP), intent(out) |
subroutine RadiationRoeweLiouCalcCO2Trans( xyz_CO2DelScaledAbsAmt, xyrr_Trans )
! USE statements
!
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
!!$ integer , intent(in ):: iband
real(DP), intent(in ):: xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(out):: xyrr_Trans (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!
! Work variables
!
real(DP):: xyrr_AbsAmt (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
integer :: k1
integer :: k2
integer :: k3
!!$ integer :: n
!!$ n = iband
!!$ n = 2
xyrr_AbsAmt(:,:,:,:) = 0.0d0
do k1 = 0, kmax
do k2 = k1, kmax
do k3 = k1+1, k2
xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) + xyz_CO2DelScaledAbsAmt(:,:,k3)
end do
xyrr_AbsAmt(:,:,k1,k2) = xyrr_AbsAmt(:,:,k1,k2) * DiffFactor
xyrr_Trans(:,:,k1,k2) = - a_CO2BandParam(3) * xyrr_AbsAmt(:,:,k1,k2) / sqrt( 1.0d0 + a_CO2BandParam(3) / a_CO2BandParam(4) * xyrr_AbsAmt(:,:,k1,k2) )
xyrr_Trans(:,:,k1,k2) = exp( xyrr_Trans(:,:,k1,k2) )
end do
end do
do k1 = 0, kmax
do k2 = 0, k1-1
xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
end do
end do
do k1 = 0, kmax
k2 = k1
xyrr_Trans(:,:,k1,k2) = 1.0d0
end do
end subroutine RadiationRoeweLiouCalcCO2Trans
| Subroutine : |
This procedure input/output NAMELIST#radiation_RL78_nml .
subroutine RadiationRoeweLiouInit
!
! Grid points settings
!
use gridset, only: imax, jmax, kmax !
! Number of vertical level
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! 日付および時刻の取り扱い
! Date and time handler
!
use dc_date, only: DCDiffTimeCreate
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid
real(DP) :: DelTimeCalcTransValue
character(STRING) :: DelTimeCalcTransUnit
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
logical :: FlagCheckLoopExit
real(DP) :: xy_TempTMP (0:imax-1, 1:jmax)
real(DP) :: xy_PF (0:imax-1, 1:jmax)
real(DP) :: xy_DPFDT (0:imax-1, 1:jmax)
real(DP) :: xy_PFTable (0:imax-1, 1:jmax)
real(DP) :: xy_DPFDTTable(0:imax-1, 1:jmax)
real(DP) :: ErrorPFInteg
real(DP), parameter:: ThresholdErrorPFInteg = 1.0d-3
! Threshold for checking accuracy of calculation of
! integrated Planc function by using a pre-calculated
! table.
integer:: i
integer:: j
integer:: l
integer:: n
namelist /radiation_RL78_nml/ VMRCO2, DelTimeCalcTransValue, DelTimeCalcTransUnit
VMRCO2 = 382.0d-6
DelTimeCalcTransValue = 3.0
DelTimeCalcTransUnit = 'hrs'
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)
rewind( unit_nml )
read( unit_nml, nml = radiation_RL78_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
end if
! Handle interval time
!
call DCDiffTimeCreate( IntTimeSave, DelTimeCalcTransValue, DelTimeCalcTransUnit ) ! (in)
!!$ do n = 1, nbmax
!!$ ! unit conversion from (cm-1) to (m-1)
!!$ a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
!!$ a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
!!$ ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$ a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
!!$ end do
do n = 1, 12-1
! unit conversion from (cm-1) to (m-1)
a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
! unit conversion from (g-1 cm2) to (kg-1 m2)
a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
end do
! O3 unit conversion
do n = 12, 21
! unit conversion from (cm-1) to (m-1)
a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
! unit conversion from "wrong (g-1 cm2)" to (kg-1 m2)
! The values written by Roewe and Liou (1978) seems to be wrong.
! In order to convert cm-1 atm-1 to g-1 cm2, the value has to be MULTIPLIED
! by ~400.
! The factor ~400 is
! 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0 * 1.0d4 * 1.0d-3.
! However, it seems that the value was DIVIDED by ~400 to convert the value
! in cm-1 atm-1 to g-1 cm2 by Roewe and Liou (1978).
! So, first, the value written in Roewe and Liou (1978) is multiplied by ~400
! to obtain the correct value in cm-1 atm-1.
! Then the value is multiplied by ~40 to convert unit from cm-1 atm-1 to m2 kg-1.
! The factor ~40 is
! 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0.
a_BandParam(3,n) = a_BandParam(3,n) * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0 * 1.0d4 * 1.0d-3 * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0
end do
do n = 21+1, nbmax
! unit conversion from (cm-1) to (m-1)
a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
! unit conversion from (g-1 cm2) to (kg-1 m2)
a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
end do
!!$ do n = 1, nbH2Omax
!!$ ! unit conversion from (cm-1) to (m-1)
!!$ a_H2OBandParam(1,n) = a_H2OBandParam(1,n) * 1.0d2
!!$ a_H2OBandParam(2,n) = a_H2OBandParam(2,n) * 1.0d2
!!$ ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$ a_H2OBandParam(3,n) = a_H2OBandParam(3,n) * 1.0d-4 * 1.0d3
!!$
!!$ ! unit conversion from (cm-1) to (m-1)
!!$ a_H2OContParam(1,n) = a_H2OContParam(1,n) * 1.0d2
!!$ a_H2OContParam(2,n) = a_H2OContParam(2,n) * 1.0d2
!!$ ! unit conversion from (g-1 cm2 atm-1) to (kg-1 m2 Pa-1)
!!$ a_H2OContParam(3,n) = a_H2OContParam(3,n) * 1.0d-4 * 1.0d3 * 1.0d-5
!!$ end do
! unit conversion from (cm-1) to (m-1)
a_CO2BandParam(1) = a_CO2BandParam(1) * 1.0d2
a_CO2BandParam(2) = a_CO2BandParam(2) * 1.0d2
! unit conversion from (g-1 cm2) to (kg-1 m2)
a_CO2BandParam(3) = a_CO2BandParam(3) * 1.0d-4 * 1.0d3
!!$ do n = 1, nbO3max
!!$ ! unit conversion from (cm-1) to (m-1)
!!$ a_O3BandParam(1,n) = a_O3BandParam(1,n) * 1.0d2
!!$ a_O3BandParam(2,n) = a_O3BandParam(2,n) * 1.0d2
!!$ ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$ a_O3BandParam(3,n) = a_O3BandParam(3,n) * 1.0d-4 * 1.0d3
!!$ end do
!!$ do n = 1, 3
!!$ ! unit conversion from (cm-1) to (m-1)
!!$ a_TMPFORO3BandParam(1,n) = a_TMPFORO3BandParam(1,n) * 1.0d2
!!$ a_TMPFORO3BandParam(2,n) = a_TMPFORO3BandParam(2,n) * 1.0d2
!!$ ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$ a_TMPFORO3BandParam(3,n) = a_TMPFORO3BandParam(3,n) * 1.0d-4 * 1.0d3
!!$ end do
! Values by Miyoshi and Morita (19??)
!
H2OScaleIndex = 1.0d0
CO2ScaleIndex = 0.9d0
O3ScaleIndex = 0.4d0
NGaussQuad = 5
MeanMolWeight = 28.0d-3
H2OMolWeight = 18.0d-3
CO2MolWeight = 44.0d-3
allocate( xyrra_TransSave (0:imax-1, 1:jmax, 0:kmax, 0:kmax, 1:nbmax) )
!----------------------------------------------------
! Check accuracy of integration of Planc function by using a pre-calculated table.
!
! This routine is called once here, to initialize a pre-calculated table.
n = 1
xy_TempTMP = 300.0d0
call CalcIntegratedPFWithTable2D( n, imax, jmax, xy_TempTMP, xy_PFTable, .false. )
do n = 1, nbmax
FlagCheckLoopExit = .false.
l = 1
do
do j = 1, jmax
do i = 0, imax-1
xy_TempTMP(i,j) = a_TableTemp(1) + ( a_TableTemp(2) - a_TableTemp(1) ) * 0.5d0 + ( a_TableTemp(2) - a_TableTemp(1) ) * ( imax * jmax * ( l - 1 ) + imax * ( j - 1 ) + i )
end do
end do
do j = 1, jmax
do i = 0, imax-1
if ( xy_TempTMP(i,j) > a_TableTemp(ntmax) ) then
xy_TempTMP(i,j) = a_TableTemp(ntmax)
FlagCheckLoopExit = .true.
end if
end do
end do
call CalcIntegratedPF2D( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, imax, jmax, xy_TempTMP, xy_PF )
call Integ_DPFDT_GQ_Array2D( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, xy_TempTMP, xy_DPFDT )
call CalcIntegratedPFWithTable2D( n, imax, jmax, xy_TempTMP, xy_PFTable, .false. )
call CalcIntegratedPFWithTable2D( n, imax, jmax, xy_TempTMP, xy_DPFDTTable, .true. )
do j = 1, jmax
do i = 0, imax-1
ErrorPFInteg = abs( xy_PF (i,j) - xy_PFTable (i,j) ) / xy_PF (i,j)
if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
call MessageNotify( 'E', module_name, 'Error of integrated PF, %f, is greater than threshold, %f, in band %d.', d = (/ ErrorPFInteg, ThresholdErrorPFInteg /), i = (/n/) )
end if
ErrorPFInteg = abs( xy_DPFDT(i,j) - xy_DPFDTTable(i,j) ) / xy_DPFDT(i,j)
if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
call MessageNotify( 'E', module_name, 'Error of integrated DPFDT, %f, is greater than threshold, %f, in band %d', d = (/ ErrorPFInteg, ThresholdErrorPFInteg /), i = (/n/) )
end if
end do
end do
if ( FlagCheckLoopExit ) exit
l = l + 1
end do
end do
!----------------------------------------------------
do n = 1, nRefAtm
RefAtmPro(2,n) = RefAtmPro(2,n) * 1.0d2
!!$ RefAtmPro(5,n) = RefAtmPro(5,n) * 1.0d-3 / ( 48.0d0 * 1.66d-27 ) &
!!$ & / ( RefAtmPro(2,n) / ( 1.38d-23 * RefAtmPro(3,n) ) )
RefAtmPro(5,n) = RefAtmPro(5,n) * 1.0d-3 / ( RefAtmPro(2,n) / ( 8.314d0 / 48.0d-3 * RefAtmPro(3,n) ) )
end do
sw_fs = .false.
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' DelTimeCalcTrans = %f [%c]', d = (/ DelTimeCalcTransValue /), c1 = trim( DelTimeCalcTransUnit ) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
end subroutine RadiationRoeweLiouInit
| Subroutine : | |
| index : | integer , intent(in ) |
| sw_log : | logical , intent(in ) |
| xyz_Press(0:imax-1,1:jmax,1:kmax) : | real(DP), intent(in ) |
| xyz_Array(0:imax-1,1:jmax,1:kmax) : | real(DP), intent(out) |
subroutine SetRefAtmPro( index, sw_log, xyz_Press, xyz_Array )
use dc_types
use gridset, only: imax, jmax, kmax
integer , intent(in ) :: index
logical , intent(in ) :: sw_log
real(DP), intent(in ) :: xyz_Press(0:imax-1,1:jmax,1:kmax)
real(DP), intent(out) :: xyz_Array(0:imax-1,1:jmax,1:kmax)
!
! local variables
!
real(DP):: z_RefAtmProArray(1:nRefAtm)
integer :: k, kk
if ( sw_log ) then
do k = 1, nRefAtm
z_RefAtmProArray(k) = log( RefAtmPro(index,k) + 1.0d-100 )
end do
else
do k = 1, nRefAtm
z_RefAtmProArray(k) = RefAtmPro(index,k)
end do
end if
do k = 1, kmax
if( xyz_Press(0,1,k) <= RefAtmPro(2,nRefAtm) ) then
xyz_Array(0,1,k) = z_RefAtmProArray(k)
else
search_loop : do kk = 2, nRefAtm
if( RefAtmPro(2,kk) < xyz_Press(0,1,k) ) exit search_loop
end do search_loop
if( kk > nRefAtm ) stop 'Unexpected error in setting temperature profile'
xyz_Array(0,1,k) = ( z_RefAtmProArray( kk ) - z_RefAtmProArray( kk-1 ) ) / ( log( RefAtmPro(2,kk) / RefAtmPro(2,kk-1) ) ) * ( log( xyz_Press(0,1,k) / RefAtmPro(2,kk-1) ) ) + z_RefAtmProArray( kk-1 )
end if
end do
if ( sw_log ) then
do k = 1, kmax
xyz_Array(0,1,k) = exp( xyz_Array(0,1,k) )
end do
else
do k = 1, kmax
xyz_Array(0,1,k) = xyz_Array(0,1,k)
end do
end if
do k = 1, kmax
xyz_Array(:,:,k) = xyz_Array(0,1,k)
end do
end subroutine SetRefAtmPro
| Subroutine : | |
| x1 : | real(DP), intent(in ) |
| x2 : | real(DP), intent(in ) |
| n : | integer , intent(in ) |
| x(n) : | real(DP), intent(out) |
| w(n) : | real(DP), intent(out) |
subroutine gauleg( x1, x2, n, x, w )
!
! Physical constants settings
!
use constants, only: PI ! $ \pi $ .
! Circular constant
real(DP), intent(in ) :: x1,x2
integer , intent(in ) :: n
real(DP), intent(out) :: x(n),w(n)
real(DP):: eps
real(DP):: z1, z, xm, xl, pp, p3, p2, p1
integer :: i, j, m
!
! changed at 2005/09/14
!
!!$ eps = 1.0e-11
eps = 1.0e-15
m = ( n+1 ) / 2
xm = ( x2+x1 ) / 2.0d0
xl = ( x2-x1 ) / 2.0d0
do i = 1, m
z = cos( pi*dble(i-0.25)/(dble(n+0.5)) )
100 continue
p1 = 1.0d0
p2 = 0.0d0
!!$ do j = 1, n
do j = 1, n
p3 = p2
p2 = p1
p1 = ( (2.0d0*dble(j)-1.0d0)*z*p2-(dble(j)-1.0d0)*p3 ) / dble(j)
enddo
pp = dble(n) * (z*p1-p2) / (z*z-1.0d0)
z1 = z
z = z1 - p1 / pp
!
! changed at 2005/09/14
!
!!$ if( abs( z-z1 ) .gt. eps ) go to 100
!!$ if( abs( z-z1 ) / z1 .gt. eps ) go to 100
!
! changed before 2008/08/03
!
if( abs( z-z1 ) / abs( z1 + 1.0d-200 ) .gt. eps ) go to 100
x( i ) = xm-xl*z
x( n+1-i ) = xm+xl*z
w( i ) = 2.0d0*xl/((1.0d0-z*z)*pp*pp)
w( n+1-i ) = w(i)
end do
end subroutine gauleg
| Constant : | |||
| module_name = ‘radiation_RL78‘ : | character(*), parameter
|
| Function : | |
| pf : | real(DP) |
| x : | real(DP), intent(in) |
| t : | real(DP), intent(in) |
real(DP) function pf(x,t)
real(DP), intent(in) :: x, t
real(dp), parameter :: c = 2.99792458d8 , planc = 6.6260755d-34 , boltz = 1.380658d-23
pf=2.0d0*planc*c*c*x*x*x /(exp(planc*c*(x+1.0d-10)/(boltz*t))-1.0d0)
return
end function pf
| Constant : | |||
| version = ’$Name: dcpam5-20100224 $’ // ’$Id: radiation_RL78.f90,v 1.4 2010-02-24 08:20:40 yot Exp $’ : | character(*), parameter
|