Class | rad_CK1991 |
In: |
radiation/rad_CK1991.F90
|
Note that Japanese and English are described in parallel.
長波放射モデル.
This is a model of long wave radiation.
Chou, M.-D., and L. Kouvaris, Calculations of transmittion functions in the infrared CO2 and O3 bands, J. Geophys. Res., 96, 9003-9012, 1991.
!$ ! 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) |
!$ ! NAMELIST#radiation_DennouAGCM_nml
Subroutine : | |
xyz_DelAbsMass(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 ) |
Spec : | character(len=*), intent(in ) |
xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) : | real(DP) , intent(out) |
subroutine RadCK1991CalcTrans( xyz_DelAbsMass, xyz_Press, xyz_Temp, Spec, xyrr_Trans ) ! USE statements ! real(DP) , intent(in ) :: xyz_DelAbsMass(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) character(len=*), intent(in ) :: Spec real(DP) , intent(out) :: xyrr_Trans (0:imax-1, 1:jmax, 0:kmax, 0:kmax) ! 初期化確認 ! Initialization check ! if ( .not. rad_ck1991_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if if ( Spec == 'CO2' ) then call RadCK1991CalcTransCore( NTabPress, NCO2TabAbsAmt, a_TabLog10Press, a_CO2TabLog10AbsAmt, aa_CO2TabAlpha, aa_CO2TabBeta, aa_CO2TabAbs, xyz_DelAbsMass, xyz_Press, xyz_Temp, xyrr_Trans ) else if ( Spec == 'O3' ) then call RadCK1991CalcTransCore( NTabPress, NO3TabAbsAmt, a_TabLog10Press, a_O3TabLog10AbsAmt, aa_O3TabAlpha, aa_O3TabBeta, aa_O3TabAbs, xyz_DelAbsMass, xyz_Press, xyz_Temp, xyrr_Trans ) else call MessageNotify( 'E', module_name, 'Specified composition, %c, is inappropriate', c1 = trim(Spec) ) end if end subroutine RadCK1991CalcTrans
Subroutine : |
subroutine RadCK1991Init ! ファイル入出力補助 ! File I/O support ! use dc_iounit, only: FileOpen !!$ ! NAMELIST ファイル入力に関するユーティリティ !!$ ! Utilities for NAMELIST file input !!$ ! !!$ use namelist_util, only: NmlutilMsg, NmlutilAryValid !!$ integer:: unit_nml ! NAMELIST ファイルオープン用装置番号. !!$ ! Unit number for NAMELIST file open !!$ integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT. !!$ ! IOSTAT of NAMELIST read !!$ namelist /rad_RL78_nml/ & !!$ & VMRCO2, & !!$ & DelTimeCalcTransValue, & !!$ & DelTimeCalcTransUnit, & !!$ & flag_save_time if ( rad_CK1991_inited ) return !!$ !!$ VMRCO2 = 382.0d-6 !!$ !!$ DelTimeCalcTransValue = 3.0 !!$ DelTimeCalcTransUnit = 'hrs' !!$ flag_save_time = .false. !!$ !!$ !!$ ! NAMELIST is input !!$ ! !!$ if ( trim(namelist_filename) /= '' ) then !!$ call FileOpen( unit_nml, & ! (out) !!$ & namelist_filename, mode = 'r' ) ! (in) !!$ !!$ rewind( unit_nml ) !!$ read( unit_nml, & ! (in) !!$ & nml = rad_RL78_nml, & ! (out) !!$ & iostat = iostat_nml ) ! (out) !!$ close( unit_nml ) !!$ !!$ call NmlutilMsg( iostat_nml, module_name ) ! (in) !!$ end if ! Convert unit of pressure from mbar to Pa ! a_TabLog10Press = 10.0d0**a_TabLog10Press a_TabLog10Press = a_TabLog10Press * 1.0d2 a_TabLog10Press = log10( a_TabLog10Press ) ! Convert unit of absorber amount from (atm cm)_{STP} to kg m-2 ! To convert from (atm cm)_{STP} to kg m-2, the value is divided by ! 1.0d2 / 101325.0d0 * 8.31432d0 / ( 44.0d-3 ) * 273.15d0, and ! 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 273.15d0, ! for CO2 and O3, respectively. ! MEMO: In a calculation below, GasRUniv variable in constants module can be used. ! But, I do not use it, since unit of value is just converted by ! multiplying a factor. Of course, non-use of GasRUniv must not cause ! significant effect on result. ! a_CO2TabLog10AbsAmt = 10.0d0**a_CO2TabLog10AbsAmt a_CO2TabLog10AbsAmt = a_CO2TabLog10AbsAmt / ( 1.0d2 / 101325.0d0 * 8.31432d0 / ( 44.0d-3 ) * 273.15d0 ) a_CO2TabLog10AbsAmt = log10( a_CO2TabLog10AbsAmt ) ! a_O3TabLog10AbsAmt = 10.0d0**a_O3TabLog10AbsAmt a_O3TabLog10AbsAmt = a_O3TabLog10AbsAmt / ( 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 273.15d0 ) a_O3TabLog10AbsAmt = log10( a_O3TabLog10AbsAmt ) ! ! Convert values absorptance from -1e3 * log10(A) to log10(A) ! aa_CO2TabAbs = 10.0d0**( aa_CO2TabAbs / ( -1.0d3 ) ) aa_CO2TabAbs = log10( aa_CO2TabAbs ) ! aa_O3TabAbs = 10.0d0**( aa_O3TabAbs / ( -1.0d3 ) ) aa_O3TabAbs = log10( aa_O3TabAbs ) ! ! Convert values of alpha absorptance from 1e4 * alpha to alpha ! aa_CO2TabAlpha = aa_CO2TabAlpha / ( 1.0d4 ) ! aa_O3TabAlpha = aa_O3TabAlpha / ( 1.0d4 ) ! ! Convert values of beta absorptance from 1e6 * beta to beta ! aa_CO2TabBeta = aa_CO2TabBeta / ( 1.0d6 ) ! aa_O3TabBeta = aa_O3TabBeta / ( 1.0d6 ) ! 印字 ; 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) ) rad_ck1991_inited = .true. end subroutine RadCK1991Init
Subroutine : | |
NTabPress : | integer , intent(in ) |
NTabAbsAmt : | integer , intent(in ) |
a_TabLog10Press(NTabPress ) : | real(DP), intent(in ) |
a_TabLog10AbsAmt(NTabAbsAmt) : | real(DP), intent(in ) |
aa_Tab(NTabAbsAmt, NTabPress) : | real(DP), intent(in ) |
xy_IndexPress(0:imax-1, 1:jmax) : | integer , intent(in ) |
xy_IndexAbsAmt(0:imax-1, 1:jmax) : | integer , intent(in ) |
xy_Log10Press(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
xy_Log10AbsAmt(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
xy_Array(0:imax-1, 1:jmax) : | real(DP), intent(out) |
subroutine RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_Tab, xy_IndexPress, xy_IndexAbsAmt, xy_Log10Press, xy_Log10AbsAmt, xy_Array ) ! USE statements ! integer , intent(in ) :: NTabPress integer , intent(in ) :: NTabAbsAmt real(DP), intent(in ) :: a_TabLog10Press (NTabPress ) real(DP), intent(in ) :: a_TabLog10AbsAmt(NTabAbsAmt) real(DP), intent(in ) :: aa_Tab (NTabAbsAmt, NTabPress) integer , intent(in ) :: xy_IndexPress (0:imax-1, 1:jmax) integer , intent(in ) :: xy_IndexAbsAmt (0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_Log10Press (0:imax-1, 1:jmax) real(DP), intent(in ) :: xy_Log10AbsAmt (0:imax-1, 1:jmax) real(DP), intent(out) :: xy_Array (0:imax-1, 1:jmax) ! ! Work variables ! real(DP) :: val1 real(DP) :: val2 integer :: ip1 integer :: ip2 integer :: iw1 integer :: iw2 integer :: i integer :: j ! 初期化確認 ! Initialization check ! if ( .not. rad_ck1991_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if do j = 1, jmax do i = 0, imax-1 ip1 = xy_IndexPress (i,j) - 1 ip2 = xy_IndexPress (i,j) iw1 = xy_IndexAbsAmt(i,j) - 1 iw2 = xy_IndexAbsAmt(i,j) val1 = ( aa_Tab (iw2,ip1) - aa_Tab (iw1,ip1) ) / ( a_TabLog10AbsAmt(iw2) - a_TabLog10AbsAmt(iw1) ) * ( xy_Log10AbsAmt (i,j) - a_TabLog10AbsAmt(iw1) ) + aa_Tab (iw1,ip1) val2 = ( aa_Tab (iw2,ip2) - aa_Tab (iw1,ip2) ) / ( a_TabLog10AbsAmt(iw2) - a_TabLog10AbsAmt(iw1) ) * ( xy_Log10AbsAmt (i,j) - a_TabLog10AbsAmt(iw1) ) + aa_Tab (iw1,ip2) xy_Array(i,j) = ( val2 - val1 ) / ( a_TabLog10Press(ip2) - a_TabLog10Press(ip1) ) * ( xy_Log10Press (i,j) - a_TabLog10Press(ip1) ) + val1 end do end do end subroutine RadCK1991Interpolate
Subroutine : | |
NTabPress : | integer , intent(in ) |
NTabAbsAmt : | integer , intent(in ) |
a_TabLog10Press(NTabPress ) : | real(DP), intent(in ) |
a_TabLog10AbsAmt(NTabAbsAmt) : | real(DP), intent(in ) |
aa_TabAlpha(NTabAbsAmt, NTabPress) : | real(DP), intent(in ) |
aa_TabBeta(NTabAbsAmt, NTabPress) : | real(DP), intent(in ) |
aa_TabAbs(NTabAbsAmt, NTabPress) : | real(DP), intent(in ) |
xyz_DelAbsMass(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 ) |
xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) : | real(DP), intent(out) |
subroutine RadCK1991CalcTransCore( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabAlpha, aa_TabBeta, aa_TabAbs, xyz_DelAbsMass, xyz_Press, xyz_Temp, xyrr_Trans ) ! USE statements ! integer , intent(in ) :: NTabPress integer , intent(in ) :: NTabAbsAmt real(DP), intent(in ) :: a_TabLog10Press (NTabPress ) real(DP), intent(in ) :: a_TabLog10AbsAmt(NTabAbsAmt) real(DP), intent(in ) :: aa_TabAlpha (NTabAbsAmt, NTabPress) real(DP), intent(in ) :: aa_TabBeta (NTabAbsAmt, NTabPress) real(DP), intent(in ) :: aa_TabAbs (NTabAbsAmt, NTabPress) real(DP), intent(in ) :: xyz_DelAbsMass (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(out) :: xyrr_Trans (0:imax-1, 1:jmax, 0:kmax, 0:kmax) ! ! Work variables ! real(DP), parameter :: RefTemp = 250.0_DP real(DP) :: xy_EffTemp (0:imax-1, 1:jmax) real(DP) :: xy_EffPress (0:imax-1, 1:jmax) real(DP) :: xy_AbsAmt (0:imax-1, 1:jmax) real(DP) :: xy_Log10EffPress(0:imax-1, 1:jmax) real(DP) :: xy_Log10AbsAmt (0:imax-1, 1:jmax) integer :: xy_IndexPress (0:imax-1, 1:jmax) integer :: xy_IndexAbsAmt (0:imax-1, 1:jmax) real(DP) :: xy_Alpha (0:imax-1, 1:jmax) real(DP) :: xy_Beta (0:imax-1, 1:jmax) real(DP) :: xy_AbsAtRefTemp (0:imax-1, 1:jmax) real(DP) :: xy_Abs (0:imax-1, 1:jmax) integer :: i integer :: j integer :: k integer :: kk integer :: l ! 初期化確認 ! Initialization check ! if ( .not. rad_ck1991_inited ) then call MessageNotify( 'E', module_name, 'This module has not been initialized.' ) end if do k = 0, kmax do kk = k+1, kmax xy_EffTemp = 0.0_DP xy_EffPress = 0.0_DP xy_AbsAmt = 0.0_DP do l = k+1, kk xy_EffTemp = xy_EffTemp + xyz_Temp (:,:,l) * xyz_DelAbsMass(:,:,l) xy_EffPress = xy_EffPress + xyz_Press(:,:,l) * xyz_DelAbsMass(:,:,l) xy_AbsAmt = xy_AbsAmt + xyz_DelAbsMass(:,:,l) end do xy_EffTemp = xy_EffTemp / ( xy_AbsAmt + 1.0d-100 ) xy_EffPress = xy_EffPress / ( xy_AbsAmt + 1.0d-100 ) xy_Log10EffPress = log10( xy_EffPress + 1.0d-100 ) xy_Log10AbsAmt = log10( xy_AbsAmt + 1.0d-100 ) do j = 1, jmax do i = 0, imax-1 if ( xy_Log10EffPress(i,j) < a_TabLog10Press (1) ) then !!$ call MessageNotify( 'E', module_name, & !!$ & 'at k = %d and kk = %d, Log10EffPress(%d,%d) = %f < %f. too small', & !!$ & i = (/k, kk, i, j/), & !!$ & d = (/xy_Log10EffPress(i,j), a_TabLog10Press(1)/) ) xy_IndexPress(i,j) = 2 else #ifdef VECTOR xy_IndexPress(i,j) = 1 + 1 Loop_Search_Press : do l = 1+1, NTabPress-1 if ( a_TabLog10Press (l) < xy_Log10EffPress(i,j) ) then xy_IndexPress(i,j) = l + 1 end if end do Loop_Search_Press #else Loop_Search_Press : do l = 1+1, NTabPress if ( a_TabLog10Press (l) > xy_Log10EffPress(i,j) ) then exit Loop_Search_Press end if end do Loop_Search_Press if ( l > NTabPress ) then l = NTabPress !!$ call MessageNotify( 'E', module_name, & !!$ & 'at k = %d and kk = %d, Log10EffPress(%d,%d) = %f > %f. too large', & !!$ & i = (/k, kk, i, j/), & !!$ & d = (/xy_Log10EffPress(i,j), a_TabLog10Press(NTabPress)/) ) end if xy_IndexPress(i,j) = l #endif end if if ( xy_Log10AbsAmt(i,j) < a_TabLog10AbsAmt(1) ) then !!$ call MessageNotify( 'E', module_name, & !!$ & 'at k = %d and kk = %d, Log10AbsAmt(%d,%d) = %f < %f. too small', & !!$ & i = (/k, kk, i, j/), & !!$ & d = (/xy_Log10AbsAmt(i,j), a_TabLog10AbsAmt(1)/) ) xy_IndexAbsAmt(i,j) = 2 else #ifdef VECTOR xy_IndexAbsAmt(i,j) = 1 + 1 Loop_Search_AbsAmt : do l = 1+1, NTabAbsAmt-1 if ( a_TabLog10AbsAmt(l) < xy_Log10AbsAmt(i,j) ) then xy_IndexAbsAmt(i,j) = l + 1 end if end do Loop_Search_AbsAmt #else Loop_Search_AbsAmt : do l = 1+1, NTabAbsAmt if ( a_TabLog10AbsAmt(l) > xy_Log10AbsAmt(i,j) ) then exit Loop_Search_AbsAmt end if end do Loop_Search_AbsAmt if ( l > NTabAbsAmt ) then call MessageNotify( 'E', module_name, 'at k = %d and kk = %d, Log10AbsAmt(%d,%d) = %f > %f. too large', i = (/k, kk, i, j/), d = (/xy_Log10AbsAmt(i,j), a_TabLog10AbsAmt(NTabAbsAmt)/) ) end if xy_IndexAbsAmt(i,j) = l #endif end if end do end do call RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabAlpha, xy_IndexPress, xy_IndexAbsAmt, xy_Log10EffPress, xy_Log10AbsAmt, xy_Alpha ) call RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabBeta, xy_IndexPress, xy_IndexAbsAmt, xy_Log10EffPress, xy_Log10AbsAmt, xy_Beta ) call RadCK1991Interpolate( NTabPress, NTabAbsAmt, a_TabLog10Press, a_TabLog10AbsAmt, aa_TabAbs, xy_IndexPress, xy_IndexAbsAmt, xy_Log10EffPress, xy_Log10AbsAmt, xy_AbsAtRefTemp ) xy_AbsAtRefTemp = 10.0d0**xy_AbsAtRefTemp xy_Abs = xy_AbsAtRefTemp * ( 1.0_DP + xy_Alpha * ( xy_EffTemp - RefTemp ) + xy_Beta * ( xy_EffTemp - RefTemp )**2 ) xyrr_Trans(:,:,k,kk) = 1.0_DP - xy_Abs end do end do ! ! correction ! do k = 0, kmax do kk = k+2, kmax do j = 1, jmax do i = 0, imax-1 if ( xyrr_Trans(i,j,k,kk) > xyrr_Trans(i,j,k,kk-1) ) then xyrr_Trans(i,j,k,kk) = xyrr_Trans(i,j,k,kk-1) end if end do end do end do end do do k = 0, kmax do kk = 0, k-1 xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,kk,k) end do kk = k xyrr_Trans(:,:,k,kk) = 1.0_DP end do end subroutine RadCK1991CalcTransCore
Variable : | |||
rad_ck1991_inited = .false. : | logical, save
|
Constant : | |||
version = ’$Name: dcpam5-20120921 $’ // ’$Id: rad_CK1991.F90,v 1.1 2011-11-30 03:46:07 yot Exp $’ : | character(*), parameter
|