!= プランク関数の計算 ! != Calculate Planck function ! ! Authors:: Yoshiyuki O. Takahashi ! Version:: $Id: planck_func.f90,v 1.6 2014/05/07 09:39:21 murashin Exp $ ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved. ! License:: See COPYRIGHT[link:../../../COPYRIGHT] ! module planck_func ! != プランク関数の計算 ! != Calculate Planck function ! ! Note that Japanese and English are described in parallel. ! ! ! ! ! !== Procedures List ! !!$ ! 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 ! !!$ ! NAMELIST#radiation_DennouAGCM_nml ! ! モジュール引用 ; USE statements ! ! 種別型パラメタ ! Kind type parameter ! use dc_types, only: DP, & ! 倍精度実数型. Double precision. & STRING, & ! 文字列. Strings. & TOKEN ! キーワード. Keywords. ! 宣言文 ; Declaration statements ! implicit none private ! 公開手続き ! Public procedure ! public :: aaa_PF public :: PF public :: DPFDT public :: Integ_PF_GQ_Array3D public :: Integ_PF_GQ_Array2D public :: Integ_DPFDT_GQ_Array2D public :: Integ_DPFDT_GQ_Array3D ! 公開変数 ! Public variables ! !!$ logical, save :: planck_func_inited = .false. ! 初期設定フラグ. ! Initialization flag ! 非公開変数 ! Private variables ! real(DP), parameter :: & & SOL = 2.99792458e8_DP , & & Planc = 6.6260755e-34_DP , & & Boltz = 1.380658e-23_DP character(*), parameter:: module_name = 'planck_func' ! モジュールの名称. ! Module name character(*), parameter:: version = & & '$Name: $' // & & '$Id: planck_func.f90,v 1.6 2014/05/07 09:39:21 murashin Exp $' ! モジュールのバージョン ! Module version contains !-------------------------------------------------------------------------------------- function aaa_PF( & & is, ie, js, je, ks, ke, & & WN, aaa_Temp & & ) & result( aaa_Res ) ! ! 温度, 比湿, 気圧から, 放射フラックスを計算します. ! ! Calculate radiation flux from temperature, specific humidity, and ! air pressure. ! ! モジュール引用 ; USE statements ! ! 宣言文 ; Declaration statements ! integer , intent(in) :: is integer , intent(in) :: ie integer , intent(in) :: js integer , intent(in) :: je integer , intent(in) :: ks integer , intent(in) :: ke real(DP), intent(in) :: WN real(DP), intent(in) :: aaa_Temp(is:ie, js:je, ks:ke) real(DP) :: aaa_Res (is:ie, js:je, ks:ke) ! 作業変数 ! Work variables ! ! 実行文 ; Executable statement ! aaa_Res = 2.0_DP * Planc * SOL * SOL * WN * WN * WN & / ( exp( Planc * SOL * ( WN+1.0e-10_DP ) / ( Boltz * aaa_Temp ) ) - 1.0_DP ) end function aaa_PF !-------------------------------------------------------------------------------------- function PF( WN, Temp ) result( Res ) ! ! 温度, 比湿, 気圧から, 放射フラックスを計算します. ! ! Calculate radiation flux from temperature, specific humidity, and ! air pressure. ! ! モジュール引用 ; USE statements ! ! 宣言文 ; Declaration statements ! real(DP), intent(in) :: WN real(DP), intent(in) :: Temp real(DP) :: Res ! 作業変数 ! Work variables ! real(DP) :: aaa_Temp(1,1,1) real(DP) :: aaa_Res (1,1,1) ! 実行文 ; Executable statement ! aaa_Temp(1,1,1) = Temp aaa_Res = & & aaa_PF( & & 1, 1, 1, 1, 1, 1, & & WN, aaa_Temp & & ) Res = aaa_Res(1,1,1) end function PF !-------------------------------------------------------------------------------------- subroutine Integ_PF_GQ_Array3D( & & wn1, wn2, num, & & is, ie, js, je, ks, ke, & & aaa_Temp, & & aaa_PFInted & & ) ! ガウス重み, 分点の計算 ! Calculate Gauss node and Gaussian weight ! use gauss_quad, only : GauLeg real(DP), intent(in ) :: wn1,wn2 integer , intent(in ) :: num integer , intent(in ) :: is, ie integer , intent(in ) :: js, je integer , intent(in ) :: ks, ke real(DP), intent(in ) :: aaa_Temp (is:ie, js:je, ks:ke) real(DP), intent(out) :: aaa_PFInted(is:ie, js:je, ks:ke) ! ! local variables ! real(DP):: x( num ), w( num ) integer :: l call GauLeg( wn1, wn2, num, x, w ) aaa_PFInted(:,:,:) = 0.0_DP do l = 1, num aaa_PFInted(:,:,:) = aaa_PFInted(:,:,:) & & + aaa_PF( & & is, ie, js, je, ks, ke, & & x(l), aaa_Temp & & ) & & * w( l ) end do end subroutine Integ_PF_GQ_Array3D !-------------------------------------------------------------------------------------- subroutine Integ_PF_GQ_Array2D( & & wn1, wn2, num, & & is, ie, js, je, & & xy_Temp, & & xy_PFInted & & ) real(DP), intent(in ) :: wn1,wn2 integer , intent(in ) :: num integer , intent(in ) :: is integer , intent(in ) :: ie integer , intent(in ) :: js integer , intent(in ) :: je real(DP), intent(in ) :: xy_Temp (is:ie, js:je) real(DP), intent(out) :: xy_PFInted(is:ie, js:je) ! ! local variables ! real(DP) :: xyz_Temp3d (is:ie, js:je, 1:1) real(DP) :: xyz_PFInted3d(is:ie, js:je, 1:1) xyz_Temp3d(:,:,1) = xy_Temp(:,:) call Integ_PF_GQ_Array3D( & & wn1, wn2, num, & & is, ie, js, je, 1, 1, & & xyz_Temp3d, & & xyz_PFInted3d & & ) xy_PFInted(:,:) = xyz_PFInted3d(:,:,1) end subroutine Integ_PF_GQ_Array2D !-------------------------------------------------------------------------------------- function DPFDT( & & WN, & ! (in ) & Temp & ! (in ) & ) & & result( Res ) ! USE statements ! ! 宣言文 ; Declaration statements ! real(DP), intent(in ) :: WN real(DP), intent(in ) :: Temp real(DP) :: Res ! 作業変数 ! Work variables ! real(DP) :: aaa_Temp(1,1,1) real(DP) :: aaa_Res (1,1,1) aaa_Temp(1,1,1) = Temp aaa_Res = aaa_DPFDT( & & 1, 1, 1, 1, 1, 1, & ! (in ) & WN, & ! (in ) & aaa_Temp & ! (in ) & ) Res = aaa_Res(1,1,1) end function DPFDT !-------------------------------------------------------------------------------------- function aaa_DPFDT( & & is, ie, js, je, ks, ke, & ! (in ) & WN, & ! (in ) & aaa_Temp & ! (in ) & ) & & result( aaa_Res ) ! USE statements ! integer , intent(in ) :: is integer , intent(in ) :: ie integer , intent(in ) :: js integer , intent(in ) :: je integer , intent(in ) :: ks integer , intent(in ) :: ke real(DP), intent(in ) :: WN real(DP), intent(in ) :: aaa_Temp(is:ie, js:je, ks:ke) real(DP) :: aaa_Res (is:ie, js:je, ks:ke) real(DP) :: aaa_ExpTerm(is:ie, js:je, ks:ke) real(DP) :: aaa_PF (is:ie, js:je, ks:ke) aaa_ExpTerm = exp( Planc * SOL * ( WN + 1.0e-10_DP ) / ( Boltz * aaa_Temp ) ) aaa_PF = 2.0_DP * Planc * SOL * SOL * WN * WN * WN & / ( aaa_ExpTerm - 1.0_DP ) aaa_Res = & & 1.0_DP / ( 2.0_DP * SOL * WN * WN * Boltz ) & & * ( aaa_PF / aaa_Temp )**2 & & * aaa_ExpTerm end function aaa_DPFDT !-------------------------------------------------------------------------------------- subroutine Integ_DPFDT_GQ_Array3D( & & WN1, WN2, Num, & ! (in ) & is, ie, js, je, ks, ke, aaa_Temp, & ! (in ) & aaa_DPFDTInted & ! (out) & ) ! USE statements ! ! ガウス重み, 分点の計算 ! Calculate Gauss node and Gaussian weight ! use gauss_quad, only : GauLeg real(DP), intent(in ) :: WN1 real(DP), intent(in ) :: WN2 integer , intent(in ) :: Num integer , intent(in ) :: is integer , intent(in ) :: ie integer , intent(in ) :: js integer , intent(in ) :: je integer , intent(in ) :: ks integer , intent(in ) :: ke real(DP), intent(in ) :: aaa_Temp (is:ie, js:je, ks:ke) real(DP), intent(out) :: aaa_DPFDTInted(is:ie, js:je, ks:ke) ! ! local variables ! real(DP):: GP( Num ) real(DP):: GW( Num ) integer :: l call GauLeg( WN1, WN2, Num, GP, GW ) aaa_DPFDTInted = 0.0_DP do l = 1, num aaa_DPFDTInted = aaa_DPFDTInted & & + aaa_DPFDT( & & is, ie, js, je, ks, ke, & ! (in ) & GP(l), & ! (in ) & aaa_Temp & ! (in ) & ) & & * GW(l) end do end subroutine Integ_DPFDT_GQ_Array3D !-------------------------------------------------------------------------------------- subroutine Integ_DPFDT_GQ_Array2D( & & WN1, WN2, Num, & ! (in ) & is, ie, js, je, aa_Temp, & ! (in ) & aa_DPFDTInted & ! (out) & ) ! USE statements ! real(DP), intent(in ) :: WN1 real(DP), intent(in ) :: WN2 integer , intent(in ) :: Num integer , intent(in ) :: is integer , intent(in ) :: ie integer , intent(in ) :: js integer , intent(in ) :: je real(DP), intent(in ) :: aa_Temp (is:ie, js:je) real(DP), intent(out) :: aa_DPFDTInted(is:ie, js:je) ! ! local variables ! real(DP) :: aaa_Temp (is:ie, js:je, 1:1) real(DP) :: aaa_DPFDTInted(is:ie, js:je, 1:1) aaa_Temp(:,:,1) = aa_Temp call Integ_DPFDT_GQ_Array3D( & & WN1, WN2, Num, & ! (in ) & is, ie, js, je, 1, 1, aaa_Temp, & ! (in ) & aaa_DPFDTInted & ! (out) & ) aa_DPFDTInted = aaa_DPFDTInted(:,:,1) end subroutine Integ_DPFDT_GQ_Array2D !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------- subroutine PlanckFuncPrepPFTable( & & WNs, WNe, NGaussQuad, TableTempMin, TableTempMax, ntmax, & & a_TableTemp, a_TableIPF, a_TableIDPFDT & & ) ! メッセージ出力 ! Message output ! use dc_message, only: MessageNotify ! ガウス重み, 分点の計算 ! Calculate Gauss node and Gaussian weight ! use gauss_quad, only : GauLeg real(DP), intent(in ) :: WNs real(DP), intent(in ) :: WNe integer , intent(in ) :: NGaussQuad real(DP), intent(in ) :: TableTempMin real(DP), intent(in ) :: TableTempMax integer , intent(in ) :: ntmax real(DP), intent(out) :: a_TableTemp (1:ntmax) real(DP), intent(out) :: a_TableIPF (1:ntmax) real(DP), intent(out) :: a_TableIDPFDT(1:ntmax) ! Local variables ! real(DP) :: TableTempIncrement integer :: nn real(DP), allocatable :: aa_TempTMP (:,:) real(DP), allocatable :: aa_PF (:,:) real(DP), allocatable :: aa_DPFDT (:,:) real(DP), allocatable :: aa_PFTable (:,:) real(DP), allocatable :: aa_DPFDTTable(:,:) 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 :: a_GQP(:) real(DP) , allocatable :: a_GQW(:) integer:: i integer:: j integer:: l integer:: m ! Preparation of tables for calculation of Plank function ! TableTempIncrement = ( TableTempMax - TableTempMin ) / ( ntmax - 1 ) do m = 1, ntmax a_TableTemp(m) = TableTempMin + TableTempIncrement * ( m - 1 ) end do a_TableIPF = 0.0_DP a_TableIDPFDT = 0.0_DP allocate( a_GQP(1:NGaussQuad) ) allocate( a_GQW(1:NGaussQuad) ) call GauLeg( & & WNs, WNe, NGaussQuad, & ! (in ) & a_GQP, a_GQW & ! (out) & ) do m = 1, ntmax do l = 1, NGaussQuad a_TableIPF (m) = a_TableIPF (m) + PF ( a_GQP(l), a_TableTemp(m) ) * a_GQW(l) a_TableIDPFDT(m) = a_TableIDPFDT(m) + DPFDT( a_GQP(l), a_TableTemp(m) ) * a_GQW(l) end do end do deallocate( a_GQP ) deallocate( a_GQW ) !---------------------------------------------------- ! Check accuracy of integration of Planc function by using a pre-calculated table. ! nn = ntmax-1 allocate( aa_TempTMP (1:nn, 1:1) ) allocate( aa_PF (1:nn, 1:1) ) allocate( aa_DPFDT (1:nn, 1:1) ) allocate( aa_PFTable (1:nn, 1:1) ) allocate( aa_DPFDTTable(1:nn, 1:1) ) j = 1 do i = 1, nn aa_TempTMP(i,j) = & & a_TableTemp(1) & & + ( a_TableTemp(2) - a_TableTemp(1) ) * 0.5_DP & & + ( a_TableTemp(2) - a_TableTemp(1) ) * ( i - 1 ) end do call Integ_PF_GQ_Array2D( & & WNs, WNe, NGaussQuad, & & 1, nn, 1, 1, aa_TempTMP, & & aa_PF & & ) call Integ_DPFDT_GQ_Array2D( & & WNs, WNe, NGaussQuad, & ! (in ) & 1, nn, 1, 1, aa_TempTMP, & ! (in ) & aa_DPFDT & ! (out) & ) call CalcIntegratedPFWithTable2D( & & ntmax, a_TableTemp, a_TableIPF, & & 1, nn, 1, 1, aa_TempTMP, & & aa_PFTable, & & .false. & & ) call CalcIntegratedPFWithTable2D( & & ntmax, a_TableTemp, a_TableIDPFDT, & & 1, nn, 1, 1, aa_TempTMP, & & aa_DPFDTTable, & & .true. & & ) do i = 1, nn ErrorPFInteg = abs( aa_PF (i,j) - aa_PFTable (i,j) ) / aa_PF (i,j) if ( ErrorPFInteg > ThresholdErrorPFInteg ) then call MessageNotify( 'E', module_name, & & 'Error of integrated PF, %f, is greater than threshold, %f.', & & d = (/ ErrorPFInteg, ThresholdErrorPFInteg /) ) end if ErrorPFInteg = abs( aa_DPFDT(i,j) - aa_DPFDTTable(i,j) ) / aa_DPFDT(i,j) if ( ErrorPFInteg > ThresholdErrorPFInteg ) then call MessageNotify( 'E', module_name, & & 'Error of integrated DPFDT, %f, is greater than threshold, %f.', & & d = (/ ErrorPFInteg, ThresholdErrorPFInteg /) ) end if end do deallocate( aa_TempTMP ) deallocate( aa_PF ) deallocate( aa_DPFDT ) deallocate( aa_PFTable ) deallocate( aa_DPFDTTable ) end subroutine PlanckFuncPrepPFTable !-------------------------------------------------------------------------------------- subroutine CalcIntegratedPFWithTable2D( & & ntmax, a_TableTemp, a_TableIPF, & & is, ie, js, je, xy_Temp, & & xy_IntegPF, & & flag_DPFDT & & ) ! USE statements ! integer , intent(in ) :: ntmax real(DP), intent(in ) :: a_TableTemp(1:ntmax) real(DP), intent(in ) :: a_TableIPF (1:ntmax) integer , intent(in ) :: is integer , intent(in ) :: ie integer , intent(in ) :: js integer , intent(in ) :: je real(DP), intent(in ) :: xy_Temp (is:ie, js:je) real(DP), intent(out) :: xy_IntegPF (is:ie, js:je) logical , intent(in ), optional :: flag_DPFDT ! ! local variables ! real(DP) :: xyz_Temp (is:ie, js:je, 1) real(DP) :: xyz_IntegPF(is:ie, js:je, 1) xyz_Temp(:,:,1) = xy_Temp call CalcIntegratedPFWithTable3D( & & ntmax, a_TableTemp, a_TableIPF, & & is, ie, js, je, 1, 1, xyz_Temp, & & xyz_IntegPF, & & flag_DPFDT & & ) xy_IntegPF = xyz_IntegPF(:,:,1) end subroutine CalcIntegratedPFWithTable2D !-------------------------------------------------------------------------------------- subroutine CalcIntegratedPFWithTable3D( & & ntmax, a_TableTemp, a_TableIPF, & & is, ie, js, je, ks, ke, xyz_Temp, & & xyz_IntegPF, & & flag_DPFDT & & ) ! USE statements ! ! メッセージ出力 ! Message output ! use dc_message, only: MessageNotify integer , intent(in ) :: ntmax real(DP), intent(in ) :: a_TableTemp(1:ntmax) real(DP), intent(in ) :: a_TableIPF (1:ntmax) integer , intent(in ) :: is integer , intent(in ) :: ie integer , intent(in ) :: js integer , intent(in ) :: je integer , intent(in ) :: ks integer , intent(in ) :: ke real(DP), intent(in ) :: xyz_Temp (is:ie, js:je, ks:ke) real(DP), intent(out) :: xyz_IntegPF(is:ie, js:je, ks:ke) logical , intent(in ), optional :: flag_DPFDT ! ! local variables ! real(DP) :: TableTempMin real(DP) :: TableTempMax real(DP) :: TableTempIncrement logical :: local_flag_DPFDT integer :: xyz_TempIndex(is:ie, js:je, ks:ke) integer :: i integer :: j integer :: k integer :: m TableTempMin = a_TableTemp(1) TableTempMax = a_TableTemp(ntmax) TableTempIncrement = ( TableTempMax - TableTempMin ) / ( ntmax - 1 ) do k = ks, ke do j = js, je do i = is, ie 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 = ks, ke do j = js, je do i = is, ie 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) = & & a_TableIPF(m-1) & & * ( 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 ) ) ) & & + a_TableIPF(m ) & & * ( 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 ) ) ) & & + a_TableIPF(m+1) & & * ( 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 = ks, ke do j = js, je do i = is, ie 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_TableIPF( m-1, iband ) xyz_IntegPF(i,j,k) = & & a_TableIPF(m-1) & & * ( 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 ) ) ) & & + a_TableIPF(m ) & & * ( 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 ) ) ) & & + a_TableIPF(m+1) & & * ( 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 !-------------------------------------------------------------------------------------- end module planck_func