| Class | rad_Earth_LW_V2_4 |
| In: |
radiation/rad_Earth_LW_V2_4.f90
|
Note that Japanese and English are described in parallel.
長波放射モデル.
This is a model of long wave radiation for the Earth‘s atmospehre. Radiation in the wavenumber range from 0 to 3000 cm-1 is calculated following the scheme by Chou et al. (2001). But absorptions by CFC and weak bands desinated as band 10 in Chou et al. (2001) are neglected.
Chou, M.-D., M. J. Suarez, X.-Z. Liang, and M. M.-H. Yan, A thermal infrared radiation parameterization for atmospheric studies, NASA Technical Report Series on Global Modeling and Data Assimilation, 19, NASA/TM-2001-104606, 2001.
| RadEarthLWV24Flux : | 放射フラックスの計算 |
| ———— : | ———— |
| RadEarthLWV24Flux : | Calculate radiation flux |
| Subroutine : | |
| xyz_DelCO2Mass(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_DelN2OMass(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_DelCH4Mass(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 ) |
| xy_SurfTemp(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| xyz_QCO2(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_QH2OVap(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_QN2O(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_QCH4(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_CloudCover(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_CloudWatREff(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_CloudIceREff(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| 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), intent(out) |
subroutine RadEarthLWV24Flux( xyz_DelCO2Mass, xyz_DelH2OVapMass, xyz_DelH2OLiqMass, xyz_DelH2OSolMass, xyz_DelO3Mass, xyz_DelN2OMass, xyz_DelCH4Mass, xyz_Press, xyz_Temp, xy_SurfTemp, xyz_QCO2, xyz_QH2OVap, xyz_QN2O, xyz_QCH4, xyz_CloudCover, xyz_CloudWatREff, xyz_CloudIceREff, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )
! USE statements
!
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 時刻管理
! Time control
!
use timeset, only: TimeN, TimesetClockStart, TimesetClockStop
!!$ ! Chou et al (1991) による長波放射モデル
!!$ ! Long radiation model described by Chou et al (1991)
!!$ !
!!$ use rad_C1991, only : &
!!$ & RadC1991CalcTransMAH2O
! Chou and Kouvaris (1991) による長波放射モデル
! Long radiation model described by Chou and Kouvaris (1991)
!
use rad_CK1991, only : RadCK1991CalcTrans
! Chou et al. (2001) による長波放射モデル
! Long radiation model described by Chou et al. (2001)
!
use rad_C2001, only : RadC2001CalcTransBand3CO2, RadC2001CalcTransBand3H2O, RadC2001CalcTrans, RadC2001ReduceCloudOptDep, RadC2001CalcCloudOptProp , RadC2001CalcIntegratedPF2D, RadC2001CalcIntegratedPF3D
! 散乱を無視した放射伝達方程式
! Radiative transfer equation without considering scattering
!
use rad_rte_nonscat, only : RadRTENonScat !, RadRTENonScatAnotherForm
! 雲関系ルーチン
! Cloud-related routines
!
use cloud_utils, only : CloudUtilsLocalizeCloud, CloudUtilsCalcOverlapCloudTrans
real(DP), intent(in ):: xyz_DelCO2Mass (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_DelN2OMass (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_DelCH4Mass (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 ):: xy_SurfTemp (0:imax-1, 1:jmax)
real(DP), intent(in ):: xyz_QCO2 (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_QH2OVap (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_QN2O (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_QCH4 (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_CloudCover (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_CloudWatREff (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ):: xyz_CloudIceREff (0:imax-1, 1:jmax, 1: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)
!
! Work variables
!
real(DP) :: xy_SurfEmis (0:imax-1, 1:jmax)
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_CloudWatSSA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudIceSSA (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudWatAF (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xyz_CloudIceAF (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_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) :: xyrr_TransEach (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
real(DP) :: xyr_RadFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyra_DelRadFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP) :: xyz_IntPF (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_SurfIntPF(0:imax-1, 1:jmax)
real(DP) :: xy_IntDPFDT0(0:imax-1, 1:jmax)
real(DP) :: xy_IntDPFDT1(0:imax-1, 1:jmax)
!!$ real(DP) :: xyr_RadFluxMA (0:imax-1, 1:jmax, 0:kmax)
!!$ real(DP) :: xyra_DelRadFluxMA(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP) :: xyz_IntPF2 (0:imax-1, 1:jmax, 1:kmax)
real(DP) :: xy_SurfIntPF2(0:imax-1, 1:jmax)
real(DP) :: xy_IntDPFDT02(0:imax-1, 1:jmax)
real(DP) :: xy_IntDPFDT12(0:imax-1, 1:jmax)
real(DP) :: xyr_RadUwFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_RadDwFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyra_DelRadUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP) :: xyra_DelRadDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
integer :: n
integer :: i, j, k
! 初期化確認
! Initialization check
!
if ( .not. rad_Earth_LW_V2_4_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! 計算時間計測開始
! Start measurement of computation time
!
call TimesetClockStart( module_name )
! Check whether the transmittance is saved or not
!
if ( .not. FlagTransSaved ) then
if ( TimeN - PrevTimeSave >= IntTimeSave ) then
call MessageNotify( 'M', module_name, 'Transmittance is not saved, but criterion for transmittance calculation is met.' )
else
call MessageNotify( 'M', module_name, 'Transmittance is not saved, and criterion for transmittance calculation ' // 'is not met. However, transmittance will be calculated.' )
end if
end if
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.
LOOP_band_trans: do n = 1, nbmax
xyrr_Trans = 1.0_DP
if ( n == nbmax ) then
! Now, nothing is done when n = nbmax.
else if ( n == 3 ) then
! 540-800 cm-1
! Calculation of H2O line and continuum transmittance
call RadC2001CalcTransBand3H2O( xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransEach )
xyrr_Trans = xyrr_Trans * xyrr_TransEach
! Calculation of CO2 transmittance
if ( FlagHighAlt ) then
! Transmittance calculation for middle atmospehre as well as lower atmosphere
call RadCK1991CalcTrans( xyz_DelCO2Mass, xyz_Press, xyz_Temp, 'CO2', xyrr_TransEach )
else
! Transmittance calculation for lower atmoshere
call RadC2001CalcTransBand3CO2( xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2, xyrr_TransEach )
end if
xyrr_Trans = xyrr_Trans * xyrr_TransEach
else
! Calculation of H2O continuum transmittance
call RadC2001CalcTrans( 'H2OCont', n, xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransEach )
xyrr_Trans = xyrr_Trans * xyrr_TransEach
! Calculation of H2O line transmittance
call RadC2001CalcTrans( 'H2OLine', n, xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransEach )
xyrr_Trans = xyrr_Trans * xyrr_TransEach
! Calculation of O3 transmittance
if ( n == 5 ) then
! 980-1100 cm-1
call RadCK1991CalcTrans( xyz_DelO3Mass, xyz_Press, xyz_Temp, 'O3', xyrr_TransEach )
xyrr_Trans = xyrr_Trans * xyrr_TransEach
end if
! Calculation of N2O transmittance
if ( ( n == 6 ) .or. ( n == 7 ) .or. ( n == 10 ) ) then
if ( maxval( xyz_DelN2OMass ) > 0.0_DP ) then
call RadC2001CalcTrans( 'N2OLine', n, xyz_DelN2OMass, xyz_Press, xyz_Temp, xyz_QN2O, xyrr_TransEach )
xyrr_Trans = xyrr_Trans * xyrr_TransEach
end if
end if
! Calculation of CH4 transmittance
if ( ( n == 6 ) .or. ( n == 7 ) ) then
if ( maxval( xyz_DelCH4Mass ) > 0.0_DP ) then
call RadC2001CalcTrans( 'CH4Line', n, xyz_DelCH4Mass, xyz_Press, xyz_Temp, xyz_QCH4, xyrr_TransEach )
xyrr_Trans = xyrr_Trans * xyrr_TransEach
end if
end if
! Calculation of CO2 transmittance in weak bands
if ( ( n == 4 ) .or. ( n == 5 ) .or. ( n == 10 ) ) then
call RadC2001CalcTrans( 'WeakBandCO2Line', n, xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2, xyrr_TransEach )
xyrr_Trans = xyrr_Trans * xyrr_TransEach
end if
end if
xyrra_TransSaved(:,:,:,:,n) = xyrr_Trans
end do LOOP_band_trans
!
! Calculation of transmittance of water vapor by using a method for middle
! atmosphere
!
!!$ call RadC1991CalcTransMAH2O( &
!!$ & xyr_Press, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
!!$ & xyrr_Trans & ! (out)
!!$ & )
xyrr_Trans = -1.0d100
xyrr_TransMASaved = xyrr_Trans
end if
!
! Calculate radiative flux
!
xy_SurfEmis = 1.0_DP
xyr_RadLUwFlux = 0.0_DP
xyr_RadLDwFlux = 0.0_DP
xyra_DelRadLUwFlux = 0.0_DP
xyra_DelRadLDwFlux = 0.0_DP
LOOP_band_RTE : do n = 1, nbmax
xyz_CloudREff = xyz_CloudWatREff
call RadC2001CalcCloudOptProp( 'Liquid', n, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudWatSSA, xyz_CloudWatAF )
xyz_DelCloudWatOptDep = xyz_CloudExtCoef * xyz_DelH2OLiqMass
!
xyz_CloudREff = xyz_CloudIceREff
call RadC2001CalcCloudOptProp( 'Ice', n, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudIceSSA, xyz_CloudIceAF )
xyz_DelCloudIceOptDep = xyz_CloudExtCoef * xyz_DelH2OSolMass
call RadC2001ReduceCloudOptDep( xyz_CloudWatSSA, xyz_CloudWatAF, xyz_DelCloudWatOptDep )
call RadC2001ReduceCloudOptDep( xyz_CloudIceSSA, xyz_CloudIceAF, xyz_DelCloudIceOptDep )
!
call CloudUtilsLocalizeCloud( xyz_CloudCover, xyz_DelCloudWatOptDep )
call CloudUtilsLocalizeCloud( xyz_CloudCover, xyz_DelCloudIceOptDep )
!
xyz_TransCloudOneLayer = exp( - ( xyz_DelCloudWatOptDep + xyz_DelCloudIceOptDep ) * DiffFactor )
!
call CloudUtilsCalcOverlapCloudTrans( xyz_TransCloudOneLayer, xyz_CloudCover, xyrr_TransCloud )
! Now, nothing is done when n = nbmax.
if ( n == nbmax ) cycle
xyrr_Trans = xyrra_TransSaved(:,:,:,:,n)
xyrr_Trans = xyrr_Trans * xyrr_TransCloud
call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_SurfIntPF, 1, jmax )
call CalcIntegratedPFWithTable3D( n, kmax, xyz_Temp, xyz_IntPF, 1, jmax )
call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_IntDPFDT0, 1, jmax, .true. )
call CalcIntegratedPFWithTable2D( n, xyz_Temp(:,:,1), xy_IntDPFDT1, 1, jmax, .true. )
xy_SurfIntPF = xy_SurfEmis * PI * xy_SurfIntPF
xyz_IntPF = PI * xyz_IntPF
xy_IntDPFDT0 = xy_SurfEmis * PI * xy_IntDPFDT0
xy_IntDPFDT1 = PI * xy_IntDPFDT1
! Lines below are under testing.
!
!!$ xy_SurfIntPF2 = xy_SurfIntPF
!!$ xyz_IntPF2 = xyz_IntPF
!!$ xy_IntDPFDT02 = xy_IntDPFDT0
!!$ xy_IntDPFDT12 = xy_IntDPFDT1
!!$
!!$
!!$ call RadC2001CalcIntegratedPF2D( &
!!$ & n, xy_SurfTemp, &
!!$ & xy_SurfIntPF &
!!$ & )
!!$ call RadC2001CalcIntegratedPF3D( &
!!$ & n, kmax, xyz_Temp, &
!!$ & xyz_IntPF &
!!$ & )
!!$
!!$ call RadC2001CalcIntegratedPF2D( &
!!$ & n, xy_SurfTemp, &
!!$ & xy_IntDPFDT0, &
!!$ & .true. &
!!$ & )
!!$ call RadC2001CalcIntegratedPF2D( &
!!$ & n, xyz_Temp(:,:,1), &
!!$ & xy_IntDPFDT1, &
!!$ & .true. &
!!$ & )
!!$
!!$ xy_SurfIntPF = xy_SurfEmis * xy_SurfIntPF
!!$ xyz_IntPF = xyz_IntPF
!!$ xy_IntDPFDT0 = xy_SurfEmis * xy_IntDPFDT0
!!$ xy_IntDPFDT1 = xy_IntDPFDT1
!!$
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ write( 20+n, * ) xy_SurfTemp(i,j), xy_SurfIntPF2(i,j), xy_SurfIntPF(i,j), &
!!$ & xy_IntDPFDT02(i,j), xy_IntDPFDT02(i,j)
!!$ end do
!!$ end do
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ write( 40+n, * ) xyz_Temp(i,j,1), xy_IntDPFDT12(i,j), xy_IntDPFDT12(i,j)
!!$ end do
!!$ end do
!!$ do k = 1, kmax
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$ write( 60+n, * ) xyz_Temp(i,j,k), xyz_IntPF2(i,j,k), xyz_IntPF(i,j,k)
!!$ end do
!!$ end do
!!$ end do
!!$ call flush( 20+n )
!!$ call flush( 40+n )
!!$ call flush( 60+n )
!!$ call OLD_RadRTENonScat( &
!!$ & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0, & ! (in)
!!$ & xyrr_Trans, & ! (in)
!!$ & xyr_RadFlux, xyra_DelRadFlux & ! (out)
!!$ & )
call RadRTENonScat( xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0, xyrr_Trans, xyr_RadUwFlux, xyr_RadDwFlux, xyra_DelRadUwFlux, xyra_DelRadDwFlux )
!!$ i = 0
!!$ j = jmax/2+1
!!$ do k = 0, kmax
!!$ write( 6, * ) k, xyr_RadFlux(i,j,k), &
!!$ & xyr_RadFlux(i,j,k) - ( xyr_RadUwFlux(i,j,k) - xyr_RadDwFlux(i,j,k) )
!!$ end do
!!$ do k = 0, kmax
!!$ write( 6, * ) k, &
!!$ & xyra_DelRadFlux(i,j,k,0), &
!!$ & xyra_DelRadFlux(i,j,k,0) - ( xyra_DelRadUwFlux(i,j,k,0) - xyra_DelRadDwFlux(i,j,k,0) ), &
!!$ & xyra_DelRadFlux(i,j,k,1), &
!!$ & xyra_DelRadFlux(i,j,k,1) - ( xyra_DelRadUwFlux(i,j,k,1) - xyra_DelRadDwFlux(i,j,k,1) )
!!$ end do
!!$ if ( ( n == 1 ) .or. ( n == 2 ) .or. ( n == 9 ) ) then
!!$ !
!!$ ! For bands 0-340, 340-540, 1380-1900
!!$ ! merge with flux calculated with a method for middle atmosphere
!!$ !
!!$
!!$ xyrr_Trans = xyrr_TransMASaved
!!$
!!$ call RadELWV22IntegRTE( &
!!$ & n, & ! (in )
!!$ & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans, & ! (in )
!!$ & xyr_RadFluxMA, xyra_DelRadFluxMA & ! (out)
!!$ & )
!!$
!!$ call RadDcpamELWV22CutMergeFlux( &
!!$ & xyr_Press, & ! (in)
!!$ & xyr_RadFlux, xyra_DelRadFlux, & ! (inout)
!!$ & xyr_RadFluxMA, xyra_DelRadFluxMA & ! (in) optional
!!$ & )
!!$
!!$ else if ( ( n == 4 ) .or. ( n == 6 ) .or. ( n == 8 ) ) then
!!$ !
!!$ ! For bands 800-980, 1100-1380, 1900-3000
!!$ ! flux above a pressure level is modified to be constant
!!$ !
!!$
!!$ call RadDcpamELWV22CutMergeFlux( &
!!$ & xyr_Press, & ! (in)
!!$ & xyr_RadFlux, xyra_DelRadFlux & ! (inout)
!!$ & )
!!$
!!$ end if
xyr_RadLUwFlux = xyr_RadLUwFlux + xyr_RadUwFlux
xyr_RadLDwFlux = xyr_RadLDwFlux + xyr_RadDwFlux
xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadUwFlux
xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadDwFlux
end do LOOP_band_RTE
!!$ i = 0
!!$ j = jmax / 2 + 1
!!$ write( 73, * ) xy_SurfTemp(i,j), 0.0d0, 0.0d0, xyr_Press(i,j,0)
!!$ do k = 1, kmax
!!$ write( 73, * ) xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), xyz_QO3(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_RadLFlux(i,j,k-1) - xyr_RadLFlux(i,j,k) ) &
!!$ & / ( xyr_Press(i,j,k-1) - xyr_Press(i,j,k) ) &
!!$ & / 1004.6 * 9.8, &
!!$ & xyz_Press(i,j,k)
!!$ end do
!!$ call flush( 83 )
!!$
!!$ i = 0
!!$ j = jmax / 2 + 1
!!$ do k = 0, kmax
!!$ write( 93, * ) xyr_RadLFlux(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 RadEarthLWV24Flux
| Subroutine : | |
| FlagSnow : | logical, intent(in) |
This procedure input/output NAMELIST#rad_Earth_LW_V2_4_nml .
subroutine RadEarthLWV24Init( FlagSnow )
! USE statements
!
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
! 暦と日時の取り扱い
! Calendar and Date handler
!
use dc_calendar, only: DCCalConvertByUnit
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid
! Chou and Kouvaris (1991) による長波放射モデル
! Long radiation model described by Chou and Kouvaris (1991)
!
use rad_CK1991, only : RadCK1991Init
! Chou et al. (2001) による長波放射モデル
! Long radiation model described by Chou et al. (2001)
!
use rad_C2001, only : RadC2001Init
! 散乱を無視した放射伝達方程式
! Radiative transfer equation without considering scattering
!
use rad_rte_nonscat, only : RadRTENonScatInit
! 雲関系ルーチン
! Cloud-related routines
!
use cloud_utils, only : CloudUtilsInit
logical, intent(in) :: FlagSnow
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
integer :: n
namelist /rad_Earth_LW_V2_4_nml/ FlagHighAlt, DelTimeCalcTransValue, DelTimeCalcTransUnit, flag_save_time
if ( rad_Earth_LW_V2_4_inited ) return
FlagHighAlt = .false.
DelTimeCalcTransValue = 3.0
DelTimeCalcTransUnit = 'hrs'
flag_save_time = .false.
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)
rewind( unit_nml )
read( unit_nml, nml = rad_Earth_LW_V2_4_nml, iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
end if
! Handle interval time
!
IntTimeSave = DCCalConvertByUnit( DelTimeCalcTransValue, DelTimeCalcTransUnit, 'sec' ) ! (in)
do n = 1, nbmax
! unit conversion from (cm-1) to (m-1)
aa_BandParam(1,n) = aa_BandParam(1,n) * 1.0d2
aa_BandParam(2,n) = aa_BandParam(2,n) * 1.0d2
end do
! allocate a variable for saving transmittance
!
allocate( xyrra_TransSaved (0:imax-1,1:jmax,0:kmax,0:kmax,1:nbmax) )
allocate( xyrr_TransMASaved(0:imax-1,1:jmax,0:kmax,0:kmax) )
call RadEarthLWV24PrepPFTable
! Initialization of modules used in this module
!
! Chou and Kouvaris (1991) による長波放射モデル
! Long radiation model described by Chou and Kouvaris (1991)
!
call RadCK1991Init
! Chou et al. (2001) による長波放射モデル
! Long radiation model described by Chou et al. (2001)
!
call RadC2001Init
! 散乱を無視した放射伝達方程式
! Radiative transfer equation without considering scattering
!
call RadRTENonScatInit
! 雲関系ルーチン
! Cloud-related routines
!
call CloudUtilsInit( FlagSnow )
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' FlagHighAlt = %b', l = (/ FlagHighAlt /) )
call MessageNotify( 'M', module_name, ' DelTimeCalcTrans = %f [%c]', d = (/ DelTimeCalcTransValue /), c1 = trim( DelTimeCalcTransUnit ) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
rad_Earth_LW_V2_4_inited = .true.
end subroutine RadEarthLWV24Init
| Subroutine : | |
| iband : | integer , intent(in ) |
| xy_IntegPF(0:imax-1, 1:jmax) : | real(DP), intent(out) |
| js : | integer , intent(in ) |
| je : | integer , intent(in ) |
| flag_DPFDT : | logical , intent(in ), optional |
subroutine CalcIntegratedPFWithTable2D( iband, xy_Temp, xy_IntegPF, js, je, flag_DPFDT )
! USE statements
!
integer , intent(in ) :: iband
real(DP), intent(in ) :: xy_temp (0:imax-1, 1:jmax)
real(DP), intent(out) :: xy_IntegPF(0:imax-1, 1:jmax)
integer , intent(in ) :: js
integer , intent(in ) :: je
logical , intent(in ), optional :: flag_DPFDT
!
! local variables
!
real(DP) :: xyz_Temp (0:imax-1, 1:jmax, 1)
real(DP) :: xyz_IntegPF(0:imax-1, 1:jmax, 1)
integer :: j
do j = js, je
xyz_Temp(:,j,1) = xy_Temp(:,j)
end do
call CalcIntegratedPFWithTable3D( iband, 1, xyz_temp, xyz_IntegPF, js, je, flag_DPFDT )
do j = js, je
xy_IntegPF(:,j) = xyz_IntegPF(:,j,1)
end do
end subroutine CalcIntegratedPFWithTable2D
| Subroutine : | |
| iband : | integer , intent(in ) |
| km : | integer , intent(in ) |
| xyz_temp(0:imax-1, 1:jmax, 1:km) : | real(DP), intent(in ) |
| xyz_IntegPF(0:imax-1, 1:jmax, 1:km) : | real(DP), intent(out) |
| js : | integer , intent(in ) |
| je : | integer , intent(in ) |
| flag_DPFDT : | logical , intent(in ), optional |
subroutine CalcIntegratedPFWithTable3D( iband, km, xyz_temp, xyz_IntegPF, js, je, flag_DPFDT )
! USE statements
!
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
integer , intent(in ) :: iband
integer , intent(in ) :: km
real(DP), intent(in ) :: xyz_temp (0:imax-1, 1:jmax, 1:km)
real(DP), intent(out) :: xyz_IntegPF(0:imax-1, 1:jmax, 1:km)
logical , intent(in ), optional :: flag_DPFDT
integer , intent(in ) :: js
integer , intent(in ) :: je
!
! local variables
!
logical :: local_flag_DPFDT
integer :: xyz_TempIndex(0:imax-1, 1:jmax, 1:km)
integer :: i
integer :: j
integer :: k
integer :: m
do k = 1, km
do j = js, je
do i = 0, imax-1
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 = js, je
do i = 0, imax-1
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 = js, je
do i = 0, imax-1
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
| Variable : | |||
| IntTimeSave : | real(DP), save
|
| Variable : | |||
| PrevTimeSave : | real(DP), save
|
| Subroutine : |
subroutine RadEarthLWV24PrepPFTable
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! ガウス重み, 分点の計算
! Calculate Gauss node and Gaussian weight
!
use gauss_quad, only : GauLeg
! プランク関数の計算
! Calculate Planck function
!
use planck_func, only : PF, DPFDT, Integ_PF_GQ_Array2D, Integ_DPFDT_GQ_Array2D
integer , parameter :: NGaussQuad = 5
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.
! Variables for preparation for calculation of Plank function
!
real(DP) , allocatable :: GQP(:)
real(DP) , allocatable :: GQW(:)
integer:: i
integer:: j
integer:: l
integer:: m
integer:: n
! Preparation of tables for calculation of Plank function
!
TableTempMin = 50.0d0
TableTempMax = 600.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( aa_BandParam(1,n), aa_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 )
!----------------------------------------------------
! 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, xy_TempTMP, xy_PFTable, 1, jmax, .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 Integ_PF_GQ_Array2D( aa_BandParam(1,n), aa_BandParam(2,n), NGaussQuad, 0, imax-1, 1, jmax, xy_TempTMP, xy_PF )
call Integ_DPFDT_GQ_Array2D( aa_BandParam(1,n), aa_BandParam(2,n), NGaussQuad, 0, imax-1, 1, jmax, xy_TempTMP, xy_DPFDT )
call CalcIntegratedPFWithTable2D( n, xy_TempTMP, xy_PFTable, 1, jmax, .false. )
call CalcIntegratedPFWithTable2D( n, xy_TempTMP, xy_DPFDTTable, 1, jmax, .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
end subroutine RadEarthLWV24PrepPFTable
| Constant : | |||
| module_name = ‘rad_Earth_LW_V2_4‘ : | character(*), parameter
|
| Variable : | |||
| rad_Earth_LW_V2_4_inited = .false. : | logical, save
|
| Constant : | |||
| version = ’$Name: dcpam5-20150129 $’ // ’$Id: rad_Earth_LW_V2_4.f90,v 1.3 2015/01/29 12:06:43 yot Exp $’ : | character(*), parameter
|