| Class | rad_rte_two_stream_app |
| In: |
radiation/rad_rte_two_stream_app.f90
|
Note that Japanese and English are described in parallel.
Solve radiative transfer equation in two stream approximation. Analytic solution is used to calculate radiative flux in a homogeneous atmosphere in which the single scattering albedo and the asymmetry factor are constant. Radiative transfer equation is solved numerically with the method by Toon et al. (1989) to calculate radiative flux in an inhomogeneous atmosphere.
Toon, O. B., C. P. McKay, and A. P. Ackerman, Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres, J. Geophys. Res., 94, 16287-16301, 1989.
| !$ ! 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 : | |
| xy_SurfAlbedo(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| SolarFluxAtTOA : | real(DP), intent(in ) |
| xy_InAngle(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| SSA : | real(DP), intent(in ) |
| AF : | real(DP), intent(in ) |
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax ) : | real(DP), intent(in ) |
| xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax ) : | real(DP), intent(out) |
| FlagSemiInfAtm : | logical , intent(in ), optional |
| FlagSL09 : | logical , intent(in ), optional |
subroutine OLD_RadRTETwoStreamAppHomogAtm( xy_SurfAlbedo, SolarFluxAtTOA, xy_InAngle, SSA, AF, xyr_OptDep, xyr_RadSFlux, FlagSemiInfAtm, FlagSL09 )
! Calculate radiative flux in a homogeneous scattering and absorbing atmosphere.
! Analytical solution is used for calculation of radiative flux.
! Radiative flux in a semi-infinite atmosphere is calculated if FlagSemiInfAtm
! is .true.. If FlagSemiInfAtm is not given or is .false., radiative flux in a finite
! atmosphere (bounded by the surface) is calculated.
!
! If FlagSL09 is .true., short wave radiative flux is calculated with the method by
! Schneider and Liu (2009).
!
! See Meador and Weaver (19??), Toon et al. (1989), Liou (200?), and so on for
! details of radiative transfer equation in this system.
real(DP), intent(in ) :: xy_SurfAlbedo(0:imax-1, 1:jmax)
real(DP), intent(in ) :: SolarFluxAtTOA
real(DP), intent(in ) :: xy_InAngle (0:imax-1, 1:jmax)
real(DP), intent(in ) :: SSA
real(DP), intent(in ) :: AF
real(DP), intent(in ) :: xyr_OptDep (0:imax-1, 1:jmax, 0:kmax )
real(DP), intent(out) :: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax )
logical , intent(in ), optional :: FlagSemiInfAtm
logical , intent(in ), optional :: FlagSL09
!
! cosz : cosine of solar zenith angle
! cosz2 : cosz squared
!
real(DP) :: xy_cosSZA ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInv ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )
real(DP) :: SSAAdj
real(DP) :: AFAdj
real(DP) :: xyr_OptDepAdj(0:imax-1, 1:jmax, 0:kmax )
real(DP) :: Lambda
real(DP) :: LSigma
real(DP) :: Gam1
real(DP) :: Gam2
real(DP) :: xy_Gam3 (0:imax-1, 1:jmax)
real(DP) :: xy_Gam4 (0:imax-1, 1:jmax)
real(DP) :: xyr_Trans (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_TMPVal(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_CUp (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_CDo (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xy_k1 (0:imax-1, 1:jmax)
real(DP) :: xy_k2 (0:imax-1, 1:jmax)
real(DP) :: xyr_ExpLamOptDep(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xy_DWRadSFluxAtTOA(0:imax-1, 1:jmax)
logical :: FlagSemiInfAtmLV
logical :: FlagSL09LV
integer :: i
integer :: j
integer :: k
! Check flags
!
FlagSL09LV = .false.
if ( present( FlagSL09 ) ) then
if ( FlagSL09 ) then
FlagSL09LV = .true.
end if
end if
FlagSemiInfAtmLV = .false.
if ( present( FlagSemiInfAtm ) ) then
if ( FlagSemiInfAtm ) then
FlagSemiInfAtmLV = .true.
end if
end if
if ( FlagSL09LV .and. ( .not. FlagSemiInfAtmLV ) ) then
call MessageNotify( 'E', module_name, 'FlagSemiInfAtm has to be .true. when FlagSL09 is .true.' )
end if
if ( FlagSL09LV ) then
SSAAdj = SSA
AFAdj = AF
xyr_OptDepAdj = xyr_OptDep
else
! Delta-function adjustment
!
SSAAdj = ( 1.0_DP - AF**2 ) * SSA / ( 1.0_DP - SSA * AF**2 )
AFAdj = AF / ( 1.0_DP + AF )
xyr_OptDepAdj = ( 1.0_DP - SSA * AF**2 ) * xyr_OptDep
end if
do j = 1, jmax
do i = 0, imax-1
if ( xy_InAngle(i,j) > 0.0_DP ) then
xy_cosSZA (i,j) = 1.0_DP / xy_InAngle(i,j)
xy_cosSZAInv (i,j) = xy_InAngle(i,j)
xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
else
xy_cosSZA (i,j) = 0.0_DP
xy_cosSZAInv (i,j) = 0.0_DP
xy_cosSZAInvsq(i,j) = 0.0_DP
end if
end do
end do
if ( FlagSL09LV ) then
! Coefficients for Hemispheric mean approximation
!
Gam1 = 2.0_DP - SSAAdj * ( 1.0_DP + AFAdj )
Gam2 = SSAAdj * ( 1.0_DP - AFAdj )
xy_Gam3 = 1.0d100
xy_Gam4 = 1.0d100
else
! Coefficients for Eddington approximation
!
Gam1 = ( 7.0_DP - SSAAdj * ( 4.0_DP + 3.0_DP * AFAdj ) ) / 4.0_DP
Gam2 = -( 1.0_DP - SSAAdj * ( 4.0_DP - 3.0_DP * AFAdj ) ) / 4.0_DP
xy_Gam3 = ( 2.0_DP - 3.0_DP * AFAdj * xy_cosSZA ) / 4.0_DP
xy_Gam4 = 1.0_DP - xy_Gam3
end if
Lambda = sqrt( Gam1**2 - Gam2**2 )
LSigma = Gam2 / ( Gam1 + Lambda )
do k = 0, kmax
xyr_Trans(:,:,k) = exp( - xyr_OptDepAdj(:,:,k) * xy_cosSZAInv )
end do
if ( FlagSL09LV ) then
xyr_CUp = 0.0_DP
xyr_CDo = 0.0_DP
xy_DWRadSFluxAtTOA = SolarFluxAtTOA * xy_CosSZA
else
do k = 0, kmax
xyr_TMPVal(:,:,k) = SSAAdj * SolarFluxAtTOA * xyr_Trans(:,:,k) / ( Lambda**2 - xy_cosSZAInvsq )
xyr_CUp(:,:,k) = xyr_TMPVal(:,:,k) * ( ( Gam1 - xy_cosSZAInv ) * xy_Gam3 + Gam2 * xy_Gam4 )
xyr_CDo(:,:,k) = xyr_TMPVal(:,:,k) * ( ( Gam1 + xy_cosSZAInv ) * xy_Gam4 + Gam2 * xy_Gam3 )
xy_DWRadSFluxAtTOA = 0.0_DP
end do
end if
! A variable used in the following calculation
!
xyr_ExpLamOptDep = exp( Lambda * xyr_OptDepAdj )
if ( FlagSemiInfAtmLV ) then
xy_k1 = 0.0_DP
xy_k2 = xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax)
else
xy_k2 = ( ( xy_SurfAlbedo * LSigma - 1.0_DP ) * ( xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax) ) * xyr_ExpLamOptDep(:,:,0) + LSigma * ( - xyr_CUp(:,:,0) + xy_SurfAlbedo * ( xyr_CDo(:,:,0) + SolarFluxAtTOA * xy_CosSZA * xyr_Trans(:,:,0) ) ) ) / ( ( xy_SurfAlbedo * LSigma - 1.0_DP ) * xyr_ExpLamOptDep(:,:,0) - ( xy_SurfAlbedo - LSigma ) * LSigma / xyr_ExpLamOptDep(:,:,0) )
xy_k1 = ( xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax) - xy_k2 ) / LSigma
end if
! Calculate radiative flux
!
do k = 0, kmax
xyr_RadSFlux(:,:,k) = ( 1.0_DP - LSigma ) * ( xy_k1 * xyr_ExpLamOptDep(:,:,k) - xy_k2 / xyr_ExpLamOptDep(:,:,k) ) + xyr_CUp(:,:,k) - xyr_CDo(:,:,k)
end do
if ( .not. FlagSL09LV ) then
!
! Add direct solar insolation
!
do k = 0, kmax
xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) - SolarFluxAtTOA * xyr_Trans(:,:,k) * xy_cosSZA
end do
end if
do k = 0, kmax
do j = 1, jmax
do i = 0, imax-1
if( xy_cosSZA(i,j) <= 0.0_DP ) then
xyr_RadSFlux(i,j,k) = 0.0_DP
end if
end do
end do
end do
end subroutine OLD_RadRTETwoStreamAppHomogAtm
| Subroutine : | |
| xy_SurfAlbedo(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| SolarFluxAtTOA : | real(DP), intent(in ) |
| xy_InAngle(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| SSA : | real(DP), intent(in ) |
| AF : | real(DP), intent(in ) |
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ) |
| xyr_RadSUwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) |
| xyr_RadSDwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) |
| FlagSemiInfAtm : | logical , intent(in ), optional |
| FlagSL09 : | logical , intent(in ), optional |
subroutine RadRTETwoStreamAppHomogAtm( xy_SurfAlbedo, SolarFluxAtTOA, xy_InAngle, SSA, AF, xyr_OptDep, xyr_RadSUwFlux, xyr_RadSDwFlux, FlagSemiInfAtm, FlagSL09 )
! Calculate radiative flux in a homogeneous scattering and absorbing atmosphere.
! Analytical solution is used for calculation of radiative flux.
! Radiative flux in a semi-infinite atmosphere is calculated if FlagSemiInfAtm
! is .true.. If FlagSemiInfAtm is not given or is .false., radiative flux in a finite
! atmosphere (bounded by the surface) is calculated.
!
! If FlagSL09 is .true., short wave radiative flux is calculated with the method by
! Schneider and Liu (2009).
!
! See Meador and Weaver (19??), Toon et al. (1989), Liou (200?), and so on for
! details of radiative transfer equation in this system.
real(DP), intent(in ) :: xy_SurfAlbedo(0:imax-1, 1:jmax)
real(DP), intent(in ) :: SolarFluxAtTOA
real(DP), intent(in ) :: xy_InAngle (0:imax-1, 1:jmax)
real(DP), intent(in ) :: SSA
real(DP), intent(in ) :: AF
real(DP), intent(in ) :: xyr_OptDep (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadSUwFlux(0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadSDwFlux(0:imax-1, 1:jmax, 0:kmax)
logical , intent(in ), optional :: FlagSemiInfAtm
logical , intent(in ), optional :: FlagSL09
!
! cosz : cosine of solar zenith angle
! cosz2 : cosz squared
!
real(DP) :: xy_cosSZA ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInv ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )
real(DP) :: SSAAdj
real(DP) :: AFAdj
real(DP) :: xyr_OptDepAdj(0:imax-1, 1:jmax, 0:kmax )
real(DP) :: Lambda
real(DP) :: LSigma
real(DP) :: Gam1
real(DP) :: Gam2
real(DP) :: xy_Gam3 (0:imax-1, 1:jmax)
real(DP) :: xy_Gam4 (0:imax-1, 1:jmax)
real(DP) :: xyr_Trans (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_TMPVal(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_CUp (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xyr_CDo (0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xy_k1 (0:imax-1, 1:jmax)
real(DP) :: xy_k2 (0:imax-1, 1:jmax)
real(DP) :: xyr_ExpLamOptDep(0:imax-1, 1:jmax, 0:kmax)
real(DP) :: xy_DWRadSFluxAtTOA(0:imax-1, 1:jmax)
logical :: FlagSemiInfAtmLV
logical :: FlagSL09LV
integer :: i
integer :: j
integer :: k
! 初期化確認
! Initialization check
!
if ( .not. rad_rte_two_stream_app_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! Check flags
!
FlagSL09LV = .false.
if ( present( FlagSL09 ) ) then
if ( FlagSL09 ) then
FlagSL09LV = .true.
end if
end if
FlagSemiInfAtmLV = .false.
if ( present( FlagSemiInfAtm ) ) then
if ( FlagSemiInfAtm ) then
FlagSemiInfAtmLV = .true.
end if
end if
if ( FlagSL09LV .and. ( .not. FlagSemiInfAtmLV ) ) then
call MessageNotify( 'E', module_name, 'FlagSemiInfAtm has to be .true. when FlagSL09 is .true.' )
end if
if ( FlagSL09LV ) then
SSAAdj = SSA
AFAdj = AF
xyr_OptDepAdj = xyr_OptDep
else
! Delta-function adjustment
!
SSAAdj = ( 1.0_DP - AF**2 ) * SSA / ( 1.0_DP - SSA * AF**2 )
AFAdj = AF / ( 1.0_DP + AF )
xyr_OptDepAdj = ( 1.0_DP - SSA * AF**2 ) * xyr_OptDep
end if
do j = 1, jmax
do i = 0, imax-1
if ( xy_InAngle(i,j) > 0.0_DP ) then
xy_cosSZA (i,j) = 1.0_DP / xy_InAngle(i,j)
xy_cosSZAInv (i,j) = xy_InAngle(i,j)
xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
else
xy_cosSZA (i,j) = 0.0_DP
xy_cosSZAInv (i,j) = 0.0_DP
xy_cosSZAInvsq(i,j) = 0.0_DP
end if
end do
end do
if ( FlagSL09LV ) then
! Coefficients for Hemispheric mean approximation
!
Gam1 = 2.0_DP - SSAAdj * ( 1.0_DP + AFAdj )
Gam2 = SSAAdj * ( 1.0_DP - AFAdj )
xy_Gam3 = 1.0d100
xy_Gam4 = 1.0d100
else
! Coefficients for Eddington approximation
!
Gam1 = ( 7.0_DP - SSAAdj * ( 4.0_DP + 3.0_DP * AFAdj ) ) / 4.0_DP
Gam2 = -( 1.0_DP - SSAAdj * ( 4.0_DP - 3.0_DP * AFAdj ) ) / 4.0_DP
xy_Gam3 = ( 2.0_DP - 3.0_DP * AFAdj * xy_cosSZA ) / 4.0_DP
xy_Gam4 = 1.0_DP - xy_Gam3
end if
Lambda = sqrt( Gam1**2 - Gam2**2 )
LSigma = Gam2 / ( Gam1 + Lambda )
do k = 0, kmax
xyr_Trans(:,:,k) = exp( - xyr_OptDepAdj(:,:,k) * xy_cosSZAInv )
end do
if ( FlagSL09LV ) then
xyr_CUp = 0.0_DP
xyr_CDo = 0.0_DP
xy_DWRadSFluxAtTOA = SolarFluxAtTOA * xy_CosSZA
else
do k = 0, kmax
xyr_TMPVal(:,:,k) = SSAAdj * SolarFluxAtTOA * xyr_Trans(:,:,k) / ( Lambda**2 - xy_cosSZAInvsq )
xyr_CUp(:,:,k) = xyr_TMPVal(:,:,k) * ( ( Gam1 - xy_cosSZAInv ) * xy_Gam3 + Gam2 * xy_Gam4 )
xyr_CDo(:,:,k) = xyr_TMPVal(:,:,k) * ( ( Gam1 + xy_cosSZAInv ) * xy_Gam4 + Gam2 * xy_Gam3 )
xy_DWRadSFluxAtTOA = 0.0_DP
end do
end if
! A variable used in the following calculation
!
xyr_ExpLamOptDep = exp( Lambda * xyr_OptDepAdj )
if ( FlagSemiInfAtmLV ) then
xy_k1 = 0.0_DP
xy_k2 = xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax)
else
xy_k2 = ( ( xy_SurfAlbedo * LSigma - 1.0_DP ) * ( xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax) ) * xyr_ExpLamOptDep(:,:,0) + LSigma * ( - xyr_CUp(:,:,0) + xy_SurfAlbedo * ( xyr_CDo(:,:,0) + SolarFluxAtTOA * xy_CosSZA * xyr_Trans(:,:,0) ) ) ) / ( ( xy_SurfAlbedo * LSigma - 1.0_DP ) * xyr_ExpLamOptDep(:,:,0) - ( xy_SurfAlbedo - LSigma ) * LSigma / xyr_ExpLamOptDep(:,:,0) )
xy_k1 = ( xy_DWRadSFluxAtTOA - xyr_CDo(:,:,kmax) - xy_k2 ) / LSigma
end if
! Calculate radiative flux
!
!!$ do k = 0, kmax
!!$ xyr_RadSFlux(:,:,k) = &
!!$ & ( 1.0_DP - LSigma ) &
!!$ & * ( xy_k1 * xyr_ExpLamOptDep(:,:,k) - xy_k2 / xyr_ExpLamOptDep(:,:,k) ) &
!!$ & + xyr_CUp(:,:,k) - xyr_CDo(:,:,k)
!!$ end do
do k = 0, kmax
xyr_RadSUwFlux(:,:,k) = xy_k1 * xyr_ExpLamOptDep(:,:,k) + LSigma * xy_k2 / xyr_ExpLamOptDep(:,:,k) + xyr_CUp(:,:,k)
xyr_RadSDwFlux(:,:,k) = LSigma * xy_k1 * xyr_ExpLamOptDep(:,:,k) + xy_k2 / xyr_ExpLamOptDep(:,:,k) + xyr_CDo(:,:,k)
end do
if ( .not. FlagSL09LV ) then
!
! Add direct solar insolation
!
!!$ do k = 0, kmax
!!$ xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) &
!!$ & - SolarFluxAtTOA * xyr_Trans(:,:,k) * xy_cosSZA
!!$ end do
do k = 0, kmax
xyr_RadSDwFlux(:,:,k) = xyr_RadSDwFlux(:,:,k) + SolarFluxAtTOA * xyr_Trans(:,:,k) * xy_cosSZA
end do
end if
do k = 0, kmax
do j = 1, jmax
do i = 0, imax-1
if( xy_cosSZA(i,j) <= 0.0_DP ) then
!!$ xyr_RadSFlux(i,j,k) = 0.0_DP
xyr_RadSUwFlux(i,j,k) = 0.0_DP
xyr_RadSDwFlux(i,j,k) = 0.0_DP
end if
end do
end do
end do
end subroutine RadRTETwoStreamAppHomogAtm
| Subroutine : |
subroutine RadRTETwoStreamAppInit
!!$ ! ファイル入出力補助
!!$ ! File I/O support
!!$ !
!!$ use dc_iounit, only: FileOpen
!!$
!!$ ! NAMELIST ファイル入力に関するユーティリティ
!!$ ! Utilities for NAMELIST file input
!!$ !
!!$ use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid
! ガウス重み, 分点の計算
! Calculate Gauss node and Gaussian weight
!
use gauss_quad, only : GauLeg
! 宣言文 ; Declaration statements
!
!!$ integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
!!$ ! Unit number for NAMELIST file open
!!$ integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
!!$ ! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
!!$ namelist /rad_rte_two_stream_app_nml/ &
!!$ & NGaussQuad
!
! デフォルト値については初期化手続 "rad_rte_two_stream_app#RadRTETwoStreamAppInit"
! のソースコードを参照のこと.
!
! Refer to source codes in the initialization procedure
! "rad_rte_two_stream_app#RadRTETwoStreamAppInit" for the default values.
!
if ( rad_rte_two_stream_app_inited ) return
! デフォルト値の設定
! Default values settings
!
!!$ NGaussQuad = 8
! NAMELIST の読み込み
! 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_rte_two_stream_app_nml, & ! (out)
!!$ & iostat = iostat_nml ) ! (out)
!!$ close( unit_nml )
!!$
!!$ call NmlutilMsg( iostat_nml, module_name ) ! (in)
!!$ end if
!!$ allocate( a_GQP(1:NGaussQuad) )
!!$ allocate( a_GQW(1:NGaussQuad) )
call GauLeg( 0.0_DP, 1.0_DP, NGaussQuad, a_GQP, a_GQW )
! Initialization of modules used in this module
!
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
!!$ call MessageNotify( 'M', module_name, 'NGaussQuad = %d', i = (/ NGaussQuad /) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
rad_rte_two_stream_app_inited = .true.
end subroutine RadRTETwoStreamAppInit
| Subroutine : | |
| xyz_SSA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyz_AF(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) |
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ) |
| xy_SurfAlbedo(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| xyr_PFInted(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ) |
| xy_SurfPFInted(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| xy_SurfDPFDTInted(0:imax-1, 1:jmax) : | real(DP), intent(in ) |
| xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) |
| xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) |
| xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(out) |
| xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(out) |
subroutine RadRTETwoStreamAppLW( xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted, xyr_RadUwFlux, xyr_RadDwFlux, xyra_DelRadUwFlux, xyra_DelRadDwFlux )
! USE statements
!
real(DP), intent(in ) :: xyz_SSA (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_AF (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyr_OptDep (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ) :: xy_SurfAlbedo (0:imax-1, 1:jmax)
real(DP), intent(in ) :: xyr_PFInted (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ) :: xy_SurfPFInted (0:imax-1, 1:jmax)
real(DP), intent(in ) :: xy_SurfDPFDTInted(0:imax-1, 1:jmax)
real(DP), intent(out) :: xyr_RadUwFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadDwFlux (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP), intent(out) :: xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
! Local variables
!
! 初期化確認
! Initialization check
!
if ( .not. rad_rte_two_stream_app_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
call RadRTETwoStreamAppWrapper( IDSchemeLongWave, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, xyr_PFInted = xyr_PFInted, xy_SurfPFInted = xy_SurfPFInted, xy_SurfDPFDTInted = xy_SurfDPFDTInted, xyra_DelRadUwFlux = xyra_DelRadUwFlux, xyra_DelRadDwFlux = xyra_DelRadDwFlux )
end subroutine RadRTETwoStreamAppLW
| Subroutine : | |||
| xyz_SSA( 0:imax-1, 1:jmax, 1:kmax ) : | real(DP), intent(in ) | ||
| xyz_AF( 0:imax-1, 1:jmax, 1:kmax ) : | real(DP), intent(in ) | ||
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ) | ||
| xy_SurfAlbedo( 0:imax-1, 1:jmax ) : | real(DP), intent(in ) | ||
| SolarFluxTOA : | real(DP), intent(in ) | ||
| xy_InAngle(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) | ||
| xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) |
subroutine RadRTETwoStreamAppSW( xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, SolarFluxTOA, xy_InAngle, xyr_RadUwFlux, xyr_RadDwFlux )
! USE statements
!
real(DP), intent(in ) :: xyz_SSA ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyz_AF ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyr_OptDep (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
real(DP), intent(in ) :: SolarFluxTOA
real(DP), intent(in ) :: xy_InAngle (0:imax-1, 1:jmax)
! sec (入射角).
! sec (angle of incidence)
real(DP), intent(out) :: xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax)
! Local variables
!
! 初期化確認
! Initialization check
!
if ( .not. rad_rte_two_stream_app_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
call RadRTETwoStreamAppWrapper( IDSchemeShortWave, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, SolarFluxTOA = SolarFluxTOA, xy_InAngle = xy_InAngle )
end subroutine RadRTETwoStreamAppSW
| Subroutine : | |||
| xyz_SSA( 0:imax-1, 1:jmax, 1:kmax ) : | real(DP), intent(in ) | ||
| xyz_AF( 0:imax-1, 1:jmax, 1:kmax ) : | real(DP), intent(in ) | ||
| SolarFluxTOA : | real(DP), intent(in ) | ||
| xy_SurfAlbedo( 0:imax-1, 1:jmax ) : | real(DP), intent(in ) | ||
| IDScheme : | integer , intent(in ) | ||
| xy_InAngle(0:imax-1, 1:jmax) : | real(DP), intent(in )
| ||
| xyr_OptDep( 0:imax-1, 1:jmax, 0:kmax ) : | real(DP), intent(in ) | ||
| xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) | ||
| xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) | ||
| js : | integer , intent(in ) | ||
| je : | integer , intent(in ) |
subroutine OLD_RadRTETwoStreamAppCore( xyz_SSA, xyz_AF, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDep, xyr_RadUwFlux, xyr_RadDwFlux, js, je )
! USE statements
!
!!$ use pf_module , only : pfint_gq_array
real(DP), intent(in ) :: xyz_SSA ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyz_AF ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: SolarFluxTOA
real(DP), intent(in ) :: xy_SurfAlbedo( 0:imax-1, 1:jmax )
integer , intent(in ) :: IDScheme
real(DP), intent(in ) :: xy_InAngle (0:imax-1, 1:jmax)
! sec (入射角).
! sec (angle of incidence)
real(DP), intent(in ) :: xyr_OptDep ( 0:imax-1, 1:jmax, 0:kmax )
real(DP), intent(out) :: xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax)
integer , intent(in ) :: js
integer , intent(in ) :: je
!!$ real(DP), intent(in ) :: gt ( 0:imax-1, 1:jmax, 1:kmax )
!!$ real(DP), intent(in ) :: gts ( 0:imax-1, 1:jmax )
!!$ real(DP), intent(in ) :: gph ( 0:imax-1, 1:jmax, 0:kmax )
!!$ real(DP), intent(in ) :: emis ( 0:imax-1, 1:jmax )
!!$ real(DP), intent(in ) :: wn1, wn2
!!$ integer , intent(in ) :: divnum
!!$ real(DP) , intent(out) :: &
!!$ gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$ gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )
!
! ssa2 : Delta-Function Adjusted Single-Scattering Albedo
! af2 : Delta-Function Adjusted Asymmetry Factor
! opdep2 : Delta-Function Adjusted Optical Depth
!
real(DP) :: xyz_SSA2 ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_AF2 ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyr_OpDep2( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: xyr_Trans2( 0:imax-1, 1:jmax, 0:kmax )
!
! gam? : Coefficients of Generalized Radiative Transfer Equation
!
real(DP) :: xyz_Gam1( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_Gam2( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_Gam3( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_Gam4( 0:imax-1, 1:jmax, 1:kmax )
!
! cosz : cosine of solar zenith angle
! cosz2 : cosz squared
!
real(DP) :: xy_cosSZA ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInv ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )
!
! temporary variables
!
real(DP) :: xyz_DTau( 0:imax-1, 1:jmax, 1:kmax )
integer :: i, j, k, l
integer :: ms, me
real(DP) :: xyz_Lambda ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_LSigma ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyaz_smalle( 0:imax-1, 1:jmax, 4, 1:kmax )
!
! cupb : upward C at bottom of layer
! cupt : upward C at top of layer
! cdob : downward C at bottom of layer
! cdot : downward C at top of layer
!
real(DP) :: xyz_cupb( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_cupt( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_cdob( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_cdot( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: aa_TridiagMtx1( 1:imax*jmax, 1:kmax*2 )
real(DP) :: aa_TridiagMtx2( 1:imax*jmax, 1:kmax*2 )
real(DP) :: aa_TridiagMtx3( 1:imax*jmax, 1:kmax*2 )
real(DP) :: aa_Vec ( 1:imax*jmax, 1:kmax*2 )
real(DP) :: xy_TMPVal( 0:imax-1, 1:jmax )
!!$ real(DP) :: mu1
!!$ real(DP) :: gth( im, jm, km+1 ), pfinth( im, jm, km+1 )
!!$ real(DP) :: b0( ijs:ije, 1, km ), b1( ijs:ije, 1, km )
!!$ real(DP) :: gemis
!!$ real(DP) :: inteup( 0:imax-1, 1:jmax, 0:kmax )
!!$ real(DP) :: intedo( 0:imax-1, 1:jmax, 0:kmax )
!!$ real(DP) :: tmpg, tmph, tmpj, tmpk, alp1, alp2, sig1, sig2
!!$ real(DP) :: tmpcos
!!$! integer(i4b), parameter :: quadn = 3
!!$ integer(i4b), parameter :: quadn = 5
!!$! integer(i4b), parameter :: quadn = 8
!!$ real(DP) :: qucos ( quadn ), qudcos( quadn )
! Variables for debug
!
!!$ real(DP) :: xyr_UwFluxDebug( 0:imax-1, 1:jmax, 0:kmax )
!!$ real(DP) :: xyr_DwFluxDebug( 0:imax-1, 1:jmax, 0:kmax )
!!$ call gauleg( 0.0d0, 1.0d0, quadn, qucos, qudcos )
do k = 1, kmax
do j = js, je
do i = 0, imax-1
if ( xyz_ssa(i,j,k) >= 1.0d0 ) then
call MessageNotify( 'E', module_name, 'Single Scattering Albedo = %f', d = (/ xyz_ssa(i,j,k) /) )
end if
end do
end do
end do
if( IDScheme .eq. IDSchemeShortWave ) then
do j = js, je
do i = 0, imax-1
if ( xy_InAngle(i,j) > 0.0d0 ) then
xy_cosSZA (i,j) = 1.0d0 / xy_InAngle(i,j)
xy_cosSZAInv (i,j) = xy_InAngle(i,j)
xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
else
xy_cosSZA (i,j) = 0.0d0
xy_cosSZAInv (i,j) = 0.0d0
xy_cosSZAInvsq(i,j) = 0.0d0
end if
end do
end do
!!$ gth(:,:,:) = 1.0d100
!!$ else
!!$
!!$ stop 'sw != 1 is not supported.'
!!$
!!$ do ij = ijs, ije
!!$ cosz ( ij, 1 ) = 1.0d100
!!$ cosz2( ij, 1 ) = 1.0d100
!!$ end do
!!$
!!$ do ij = ijs, ije
!!$ gth( ij, 1, 1 ) = gt( ij, 1, 1 )
!!$ end do
!!$ do k = 2, km+1-1
!!$ do ij = ijs, ije
!!$ gth( ij, 1, k ) = ( gt( ij, 1, k-1 ) + gt( ij, 1, k ) ) * 0.5d0
!!$ end do
!!$ end do
!!$ do ij = ijs, ije
!!$ gth( ij, 1, km+1 ) = gts( ij, 1 )
!!$ end do
!!$! call pfint_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$! ijs, ije )
!!$ call pfint_gq_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$ ijs, ije )
else
stop 'Unexpected sw number in twostreamapp.f90'
end if
if( IDScheme .eq. IDSchemeShortWave ) then
!
! Delta-Function Adjustment
!
do k = 1, kmax
do j = js, je
xyz_AF2 (:,j,k) = xyz_AF(:,j,k) / ( 1.0d0 + xyz_AF(:,j,k) )
xyz_SSA2(:,j,k) = ( 1.0d0 - xyz_AF(:,j,k)**2 ) * xyz_SSA(:,j,k) / ( 1.0d0 - xyz_SSA(:,j,k) * xyz_AF(:,j,k)**2 )
end do
end do
!!$ opdep2(:,:,:) = ( 1.0d0 - xyz_ssa * xyz_af**2 ) * opdep(:,:,:)
do k = 0, kmax
do j = js, je
xyr_OpDep2(:,j,k) = 0.0d0
end do
end do
do k = kmax-1, 0, -1
do j = js, je
xyr_OpDep2(:,j,k) = xyr_OpDep2(:,j,k+1) + ( 1.0d0 - xyz_SSA(:,j,k+1) * xyz_AF(:,j,k+1)**2 ) * ( xyr_OptDep(:,j,k) - xyr_OptDep(:,j,k+1) )
end do
end do
do k = 0, kmax
do j = js, je
xyr_Trans2(:,j,k) = exp( -xyr_OpDep2(:,j,k) * xy_cosSZAInv(:,j) )
end do
end do
!
! Eddington approximation
!
do k = 1, kmax
do j = js, je
xyz_Gam1(:,j,k) = ( 7.0d0 - xyz_SSA2(:,j,k) * ( 4.0d0 + 3.0d0 * xyz_AF2(:,j,k) ) ) / 4.0d0
xyz_Gam2(:,j,k) = -( 1.0d0 - xyz_SSA2(:,j,k) * ( 4.0d0 - 3.0d0 * xyz_AF2(:,j,k) ) ) / 4.0d0
xyz_Gam3(:,j,k) = ( 2.0d0 - 3.0d0 * xyz_AF2(:,j,k) * xy_cosSZA(:,j) ) / 4.0d0
xyz_Gam4(:,j,k) = 1.0d0 - xyz_Gam3(:,j,k)
end do
end do
!!$ else if( sw .eq. 2 ) then
!!$
!!$ !
!!$ ! In infrared, delta-function adjustment is not done.
!!$ !
!!$ af2 = af
!!$ ssa2 = ssa
!!$
!!$ do k = 1, km+1
!!$ do ij = ijs, ije
!!$ opdep2( ij, 1, k ) = opdep( ij, 1, k )
!!$ end do
!!$ end do
!!$
!!$ do k = 1, km+1
!!$ do ij = ijs, ije
!!$ exp_opdep2( ij, 1, k ) = 1.0d100
!!$ end do
!!$ end do
!!$
!!$ !
!!$ ! Hemispheric mean approximation
!!$ !
!!$ do k = 1, km
!!$ do ij = ijs, ije
!!$ gam1( ij, 1, k ) = 2.0d0 - ssa2 * ( 1.0d0 + af2 )
!!$ gam2( ij, 1, k ) = ssa2 * ( 1.0d0 - af2 )
!!$ gam3( ij, 1, k ) = 1.0d100
!!$ gam4( ij, 1, k ) = 1.0d100
!!$ end do
!!$ end do
else
stop 'Unexpected sw number in twostreamapp.f90'
end if
do k = 1, kmax
do j = js, je
xyz_DTau(:,j,k) = xyr_OpDep2(:,j,k-1) - xyr_OpDep2(:,j,k)
end do
end do
!
! Avoiding singularity when dtau equal to zero
!
do k = 1, kmax
do j = js, je
do i = 0, imax-1
if( ( IDScheme .eq. IDSchemeLongWave ) .and. ( xyz_DTau(i,j,k) .lt. 1.0d-10 ) ) then
xyz_DTau(i,j,k) = 0.0d0
end if
end do
end do
end do
do k = 1, kmax
do j = js, je
! In very small parameter space of single scattering albedo and asymmetry
! factor close to asymmetry factor of 1, xyz_Gam1**2 - xyz_Gam2**2 becomes
! negative value. (yot, 2011/11/20)
!!$ xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) )
!!$ xyz_Lambda(:,j,k) = sqrt( max( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k), 0.0_DP ) )
xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) + 1.0d-10 )
xyz_LSigma(:,j,k) = xyz_Gam2(:,j,k) / ( xyz_Gam1(:,j,k) + xyz_Lambda(:,j,k) )
end do
end do
do k = 1, kmax
do j = js, je
xy_TMPVal(:,j) = exp( - xyz_Lambda(:,j,k) * xyz_DTau(:,j,k) )
xyaz_smalle(:,j,1,k) = xyz_LSigma(:,j,k) * xy_TMPVal(:,j) + 1.0_DP
xyaz_smalle(:,j,2,k) = xyz_LSigma(:,j,k) * xy_TMPVal(:,j) - 1.0_DP
xyaz_smalle(:,j,3,k) = xy_TMPVal(:,j) + xyz_LSigma(:,j,k)
xyaz_smalle(:,j,4,k) = xy_TMPVal(:,j) - xyz_LSigma(:,j,k)
end do
end do
if( IDScheme .eq. IDSchemeShortWave ) then
do k = 1, kmax
do j = js, je
! Lines below will be deleted in future (yot, 2010/09/14).
!!$ xyz_cupb(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1) &
!!$ & * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
!!$ xyz_cupt(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k ) &
!!$ & * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
!!$ xyz_cdob(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1) &
!!$ & * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
!!$ xyz_cdot(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k ) &
!!$ & * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
xy_TMPVal(:,j) = xyz_SSA2(:,j,k) * SolarFluxTOA * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
xyz_cupb(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k-1)
xyz_cupt(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k )
xy_TMPVal(:,j) = xyz_SSA2(:,j,k) * SolarFluxTOA * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
xyz_cdob(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k-1)
xyz_cdot(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k )
end do
end do
!!$ else if( sw .eq. 2 ) then
!!$
!!$ do k = 1, km
!!$ do ij = ijs, ije
!!$ !
!!$ ! Notice!
!!$ ! Avoiding singularity when dtau equal to zero.
!!$ ! dtau occationally has much smaller value.
!!$ ! When this occurs, b1 cannot be calculated correctly.
!!$ !
!!$ if( dtau( ij, 1, k ) .ne. 0.0d0 ) then
!!$ b0( ij, 1, k ) = pfinth( ij, 1, k )
!!$ b1( ij, 1, k ) = ( pfinth( ij, 1, k+1 ) - pfinth( ij, 1, k ) ) / dtau( ij, 1, k )
!!$ else
!!$ b0( ij, 1, k ) = 0.0d0
!!$ b1( ij, 1, k ) = 0.0d0
!!$ end if
!!$
!!$ end do
!!$ end do
!!$
!!$ do k = 1, km
!!$ do ij = ijs, ije
!!$ mu1 = ( 1.0d0 - ssa2 ) / ( gam1( ij, 1, k ) - gam2( ij, 1, k ) )
!!$
!!$ cupb( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$ + b1( ij, 1, k ) &
!!$ * ( dtau( ij, 1, k ) + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$ cupt( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$ + b1( ij, 1, k ) &
!!$ * ( 0.0d0 + 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$ cdob( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$ + b1( ij, 1, k ) &
!!$ * ( dtau( ij, 1, k ) - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$ cdot( ij, 1, k ) = pix2 * mu1 * ( b0( ij, 1, k ) &
!!$ + b1( ij, 1, k ) &
!!$ * ( 0.0d0 - 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) ) )
!!$ end do
!!$ end do
!!$
else
stop 'Unexpected sw number in twostreamapp.f90'
end if
k = 1
l = 1
do j = js, je
do i = 0, imax-1
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = 0.0_DP
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,3,k)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = - ( xyaz_smalle(i,j,2,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,4,k) )
end do
end do
if( IDScheme .eq. IDSchemeShortWave ) then
do j = js, je
do i = 0, imax-1
aa_Vec( i+imax*(j-1)+1, l ) = - xyz_cupb(i,j,k) + xy_SurfAlbedo(i,j) * xyz_cdob(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_Trans2(i,j,0) * xy_cosSZA(i,j)
end do
end do
!!$ else if( sw .eq. 2 ) then
!!$ do ij = ijs, ije
!!$! gemis = 1.0d0
!!$ gemis = emis( ij, 1 )
!!$ ea( (ij-ijs+1), l ) &
!!$ = -cupb( ij, 1, km ) + ( 1.0d0 - gemis ) * cdob( ij, 1, km ) &
!!$ + gemis * pi * pfinth( ij, 1, km+1 )
!!$ end do
else
stop 'Unexpected sw number in twostreamapp.f90'
end if
do k = 1, kmax-1
do j = js, je
do i = 0, imax-1
l = 2 * k ! equation number
!
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,2,k+1) - xyaz_smalle(i,j,3,k ) * xyaz_smalle(i,j,4,k+1)
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k ) * xyaz_smalle(i,j,2,k+1) - xyaz_smalle(i,j,4,k ) * xyaz_smalle(i,j,4,k+1)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k+1) * xyaz_smalle(i,j,4,k+1) - xyaz_smalle(i,j,2,k+1) * xyaz_smalle(i,j,3,k+1)
aa_Vec ( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k+1) * ( - xyz_cdot(i,j,k ) + xyz_cdob(i,j,k+1) ) - xyaz_smalle(i,j,4,k+1) * ( - xyz_cupt(i,j,k ) + xyz_cupb(i,j,k+1) )
l = 2 * k + 1 ! equation number
!
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k ) * xyaz_smalle(i,j,3,k ) - xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,4,k )
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,1,k+1) - xyaz_smalle(i,j,3,k ) * xyaz_smalle(i,j,3,k+1)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,3,k ) * xyaz_smalle(i,j,4,k+1) - xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,2,k+1)
aa_Vec ( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,3,k ) * ( - xyz_cdot(i,j,k ) + xyz_cdob(i,j,k+1) ) - xyaz_smalle(i,j,1,k ) * ( - xyz_cupt(i,j,k ) + xyz_cupb(i,j,k+1) )
end do
end do
end do
k = kmax
l = 2 * kmax
do j = js, je
do i = 0, imax-1
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k)
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = 0.0d0
aa_Vec ( i+imax*(j-1)+1, l ) = -xyz_cdot(i,j,k) + 0.0d0
end do
end do
ms = 0 + imax*(js-1)+1
me = imax-1 + imax*(je-1)+1
call tridiag( imax*jmax, 2*kmax, aa_TridiagMtx1, aa_TridiagMtx2, aa_TridiagMtx3, aa_Vec, ms, me )
if( IDScheme .eq. IDSchemeShortWave ) then
k = 1
do j = js, je
do i = 0, imax-1
xyr_RadUwFlux(i,j,k-1) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) + xyz_cupb(i,j,k)
xyr_RadDwFlux(i,j,k-1) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) + xyz_cdob(i,j,k)
end do
end do
do k = 1, kmax
do j = js, je
do i = 0, imax-1
xyr_RadUwFlux(i,j,k) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) + xyz_cupt(i,j,k)
xyr_RadDwFlux(i,j,k) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) + xyz_cdot(i,j,k)
end do
end do
end do
! Code for debug
!
!!$ do k = 1, kmax
!!$ do j = js, je
!!$ do i = 0, imax-1
!!$ xyr_UwFluxDebug(i,j,k-1) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
!!$ & - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) &
!!$ & + xyz_cupb(i,j,k)
!!$ xyr_DwFluxDebug(i,j,k-1) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
!!$ & - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) &
!!$ & + xyz_cdob(i,j,k)
!!$ end do
!!$ end do
!!$ end do
!!$
!!$ k = kmax
!!$ do j = js, je
!!$ do i = 0, imax-1
!!$ xyr_UwFluxDebug(i,j,k) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
!!$ & + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) &
!!$ & + xyz_cupt(i,j,k)
!!$ xyr_DwFluxDebug(i,j,k) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
!!$ & + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) &
!!$ & + xyz_cdot(i,j,k)
!!$ end do
!!$ end do
!!$
!!$
!!$ i = imax/2
!!$ j = jmax/2
!!$ do k = kmax, 0, -1
!!$ write( 6, * ) k, xyr_UwFlux(i,j,k), xyr_UwFluxDebug(i,j,k), xyr_UwFlux(i,j,k) - xyr_UwFluxDebug(i,j,k)
!!$ end do
!!$ do k = kmax, 0, -1
!!$ write( 6, * ) k, xyr_DwFlux(i,j,k), xyr_DwFluxDebug(i,j,k), xyr_DwFlux(i,j,k) - xyr_DwFluxDebug(i,j,k)
!!$ end do
!!$ k = 0
!!$ write( 6, * ) k, xyr_UwFlux(i,j,k), xy_SurfAlbedo(i,j) * xyr_DwFlux(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_Trans2(i,j,0) * xy_cosSZA(i,j), &
!!$ xyr_UwFlux(i,j,k) - ( xy_SurfAlbedo(i,j) * xyr_DwFlux(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_Trans2(i,j,0) * xy_cosSZA(i,j) )
!!$ stop
!
! Addition of Direct Solar Insident
!
do k = 0, kmax
do j = js, je
xyr_RadDwFlux(:,j,k) = xyr_RadDwFlux(:,j,k) + SolarFluxTOA * xyr_Trans2(:,j,k) * xy_cosSZA(:,j)
end do
end do
!!$ else if( sw .eq. 2 ) then
!!$
!!$ ! Source function technique described by Toon et al. [1989]
!!$ ! is used to calculated infrared heating rate.
!!$ !
!!$ do k = 1, km+1
!!$ do ij = ijs, ije
!!$ upfl( ij, 1, k ) = 0.0d0
!!$ dofl( ij, 1, k ) = 0.0d0
!!$ end do
!!$ end do
!!$
!!$ do l = 1, quadn
!!$ tmpcos = qucos( l )
!!$
!!$ k = 1
!!$ do ij = ijs, ije
!!$ intedo( ij, 1, k ) = 0.0d0
!!$ end do
!!$ do k = 1, km
!!$ do ij = ijs, ije
!!$ tmpj = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$ * ( lambda( ij, 1, k ) + 2.0d0 )
!!$ tmpk = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) &
!!$ * ( 2.0d0 - lambda( ij, 1, k ) )
!!$ sig1 = pix2 * ( b0( ij, 1, k ) &
!!$ - b1( ij, 1, k ) &
!!$ * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$ sig2 = pix2 * b1( ij, 1, k )
!!$ intedo( ij, 1, k+1 ) = intedo( ij, 1, k ) &
!!$ * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$ + tmpj / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$ * ( 1.0d0 &
!!$ - exp( -dtau( ij, 1, k ) &
!!$ * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$ + tmpk / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$ * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$ -exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$ + sig1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$ + sig2 * ( tmpcos * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$ + dtau( ij, 1, k ) - tmpcos )
!!$ end do
!!$ end do
!!$
!!$ k = km+1
!!$ do ij = ijs, ije
!!$! gemis = 1.0d0
!!$ gemis = emis( ij, 1 )
!!$ inteup( ij, 1, k ) = ( 1.0d0 - gemis ) * intedo( ij, 1, km+1 ) &
!!$ + gemis * pix2 * pfinth( ij, 1, km+1 )
!!$ end do
!!$ do k = km, 1, -1
!!$ do ij = ijs, ije
!!$ tmpg = ( ea( (ij-ijs+1), k+k-1 ) + ea( (ij-ijs+1), k+k ) ) &
!!$ * ( 2.0d0 - lambda( ij, 1, k ) )
!!$ tmph = ( ea( (ij-ijs+1), k+k-1 ) - ea( (ij-ijs+1), k+k ) ) * lsigma( ij, 1, k ) &
!!$ * ( lambda( ij, 1, k ) + 2.0d0 )
!!$ alp1 = pix2 * ( b0( ij, 1, k ) &
!!$ + b1( ij, 1, k ) &
!!$ * ( 1.0d0 / ( gam1( ij, 1, k ) + gam2( ij, 1, k ) ) - 0.5d0 ) )
!!$ alp2 = pix2 * b1( ij, 1, k )
!!$ inteup( ij, 1, k ) = inteup( ij, 1, k+1 ) &
!!$ * exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$ + tmpg / ( lambda( ij, 1, k ) * tmpcos - 1.0d0 ) &
!!$ * ( exp( -dtau( ij, 1, k ) / tmpcos ) &
!!$ - exp( -dtau( ij, 1, k ) * lambda( ij, 1, k ) ) ) &
!!$ + tmph / ( lambda( ij, 1, k ) * tmpcos + 1.0d0 ) &
!!$ * ( 1.0d0 &
!!$ -exp( -dtau( ij, 1, k ) &
!!$ * ( lambda( ij, 1, k ) + 1.0d0 / tmpcos ) ) ) &
!!$ + alp1 * ( 1.0d0 - exp( -dtau( ij, 1, k ) / tmpcos ) ) &
!!$ + alp2 * ( tmpcos &
!!$ - ( dtau( ij, 1, k ) + tmpcos ) &
!!$ * exp( -dtau( ij, 1, k ) / tmpcos ) )
!!$ end do
!!$ end do
!!$
!!$ do k = 1, km+1
!!$ do ij = ijs, ije
!!$ upfl( ij, 1, k ) = upfl( ij, 1, k ) &
!!$ + inteup( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$ dofl( ij, 1, k ) = dofl( ij, 1, k ) &
!!$ + intedo( ij, 1, k ) * qucos( l ) * qudcos( l )
!!$ end do
!!$ end do
!!$ end do
!!$
else
stop 'Unexpected sw number in twostreamapp.f90'
end if
!!$ do ij = ijs, ije
!!$ gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$ end do
!!$ if( sw == 1 ) then
!!$ ! visible
!!$ do ij = ijs, ije
!!$ gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$ end do
!!$ else
!!$ ! ir
!!$ do ij = ijs, ije
!!$ gdf( ij, 1 ) = dofl( ij, 1, km+1 )
!!$ end do
!!$ end if
!!$ do ij = ijs, ije
!!$ gdf( ij, 1 ) = upfl( ij, 1, km+1 ) - dofl( ij, 1, km+1 )
!!$ end do
!!$ do k = 1, km
!!$ do ij = ijs, ije
!!$ q( ij, 1, k ) &
!!$ = ( ( upfl( ij, 1, k+1 ) - dofl( ij, 1, k+1 ) ) &
!!$ - ( upfl( ij, 1, k ) - dofl( ij, 1, k ) ) ) &
!!$ / ( gph ( ij, 1, k+1 ) - gph ( ij, 1, k ) ) &
!!$ * grav / cp
!!$ end do
!!$ end do
!!$ if( sw == 1 ) then
!!$ write( 6, * ) 'vis flux=', &
!!$! -gdf((ijs+ije)/2,1), &
!!$! -gdf((ijs+ije)/2,1) * ( 1.0d0 - albedo((ijs+ije)/2,1) ), &
!!$ -gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1), &
!!$! upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1),
!!$ upfl((ijs+ije)/2,1,km+1), &
!!$ upfl((ijs+ije)/2,1,km+1) - gdf((ijs+ije)/2,1) * albedo((ijs+ije)/2,1)
!!$ else
!!$ write( 6, * ) 'ir flux=', &
!!$ -gdf((ijs+ije)/2,1), &
!!$ upfl((ijs+ije)/2,1,km+1)-dofl((ijs+ije)/2,1,km+1), &
!!$ ' sig * Ts^4 = ', sboltz * gts((ijs+ije)/2,1)**4
!!$ end if
if( IDScheme .eq. IDSchemeShortWave ) then
do k = 0, kmax
do j = js, je
do i = 0, imax-1
if( xy_cosSZA(i,j) <= 0.0d0 ) then
xyr_RadUwFlux(i,j,k) = 0.0d0
xyr_RadDwFlux(i,j,k) = 0.0d0
end if
end do
end do
end do
end if
!!$ !
!!$ ! output variables
!!$ !
!!$ do ij = ijs, ije
!!$ goru( ij, 1 ) = upfl( ij, 1, 1 )
!!$ gord( ij, 1 ) = dofl( ij, 1, 1 )
!!$ gsru( ij, 1 ) = upfl( ij, 1, km+1 )
!!$ gsrd( ij, 1 ) = dofl( ij, 1, km+1 )
!!$ end do
!!$ if( sw .eq. 1 ) then
!!$ do ij = ijs, ije
!!$ if( cosz( ij, 1 ) .le. 0.0d0 ) then
!!$ goru( ij, 1 ) = 0.0d0
!!$ gord( ij, 1 ) = 0.0d0
!!$ gsru( ij, 1 ) = 0.0d0
!!$ gsrd( ij, 1 ) = 0.0d0
!!$ end if
!!$ end do
!!$ end if
!!$ do ij = ijs, ije
!!$ gor ( ij, 1 ) = goru( ij, 1 ) - gord( ij, 1 )
!!$ gsr ( ij, 1 ) = gsru( ij, 1 ) - gsrd( ij, 1 )
!!$ end do
end subroutine OLD_RadRTETwoStreamAppCore
| Subroutine : | |||
| IDScheme : | integer , intent(in ) | ||
| xyz_SSA(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) | ||
| xyz_AF(0:imax-1, 1:jmax, 1:kmax) : | real(DP), intent(in ) | ||
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ) | ||
| xy_SurfAlbedo(0:imax-1, 1:jmax ) : | real(DP), intent(in ) | ||
| xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) | ||
| xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) | ||
| js : | integer , intent(in ) | ||
| je : | integer , intent(in ) | ||
| SolarFluxTOA : | real(DP), intent(in ), optional | ||
| xy_InAngle(0:imax-1, 1:jmax) : | real(DP), intent(in ), optional
| ||
| xyr_PFInted(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ), optional | ||
| xy_SurfPFInted(0:imax-1, 1:jmax) : | real(DP), intent(in ), optional | ||
| xy_SurfDPFDTInted(0:imax-1, 1:jmax) : | real(DP), intent(in ), optional | ||
| xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(out), optional | ||
| xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(out), optional |
subroutine RadRTETwoStreamAppCore( IDScheme, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, js, je, SolarFluxTOA, xy_InAngle, xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted, xyra_DelRadUwFlux, xyra_DelRadDwFlux )
! USE statements
!
!!$ use pf_module , only : pfint_gq_array
integer , intent(in ) :: IDScheme
real(DP), intent(in ) :: xyz_SSA (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyz_AF (0:imax-1, 1:jmax, 1:kmax)
real(DP), intent(in ) :: xyr_OptDep (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ) :: xy_SurfAlbedo(0:imax-1, 1:jmax )
real(DP), intent(out) :: xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax)
integer , intent(in ) :: js
integer , intent(in ) :: je
real(DP), intent(in ), optional :: SolarFluxTOA
real(DP), intent(in ), optional :: xy_InAngle (0:imax-1, 1:jmax)
! sec (入射角).
! sec (angle of incidence)
real(DP), intent(in ), optional :: xyr_PFInted (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ), optional :: xy_SurfPFInted (0:imax-1, 1:jmax)
real(DP), intent(in ), optional :: xy_SurfDPFDTInted(0:imax-1, 1:jmax)
real(DP), intent(out), optional :: xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP), intent(out), optional :: xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
!!$ real(DP), intent(in ) :: gt ( 0:imax-1, 1:jmax, 1:kmax )
!!$ real(DP), intent(in ) :: gts ( 0:imax-1, 1:jmax )
!!$ real(DP), intent(in ) :: gph ( 0:imax-1, 1:jmax, 0:kmax )
!!$ real(DP), intent(in ) :: emis ( 0:imax-1, 1:jmax )
!!$ real(DP), intent(in ) :: wn1, wn2
!!$ integer , intent(in ) :: divnum
!!$ real(DP) , intent(out) :: &
!!$ gor ( im, jm ), goru ( im, jm ), gord ( im, jm ), &
!!$ gsr ( im, jm ), gsru ( im, jm ), gsrd ( im, jm )
!
! ssa2 : Delta-Function Adjusted Single-Scattering Albedo
! af2 : Delta-Function Adjusted Asymmetry Factor
! opdep2 : Delta-Function Adjusted Optical Depth
!
real(DP) :: xyz_SSA2 ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_AF2 ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyr_OpDep2( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: xyr_Trans2( 0:imax-1, 1:jmax, 0:kmax )
!
! gam? : Coefficients of Generalized Radiative Transfer Equation
!
real(DP) :: xyz_Gam1( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_Gam2( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_Gam3( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_Gam4( 0:imax-1, 1:jmax, 1:kmax )
!
! cosz : cosine of solar zenith angle
! cosz2 : cosz squared
!
real(DP) :: xy_cosSZA ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInv ( 0:imax-1, 1:jmax )
real(DP) :: xy_cosSZAInvsq( 0:imax-1, 1:jmax )
!
! temporary variables
!
real(DP) :: xyz_DelTau( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_Lambda ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_LGamma ( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyaz_smalle( 0:imax-1, 1:jmax, 4, 1:kmax )
real(DP) :: xy_SurfSrc( 0:imax-1, 1:jmax )
!
! CUpB : upward C at bottom of layer
! CUpT : upward C at top of layer
! CDoB : downward C at bottom of layer
! CDoT : downward C at top of layer
!
real(DP) :: xyz_CUpB( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_CUpT( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_CDoB( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_CDoT( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: aa_TridiagMtx1( 1:imax*jmax, 1:kmax*2 )
real(DP) :: aa_TridiagMtx2( 1:imax*jmax, 1:kmax*2 )
real(DP) :: aa_TridiagMtx3( 1:imax*jmax, 1:kmax*2 )
real(DP) :: aa_Vec ( 1:imax*jmax, 1:kmax*2 )
real(DP) :: xy_TMPVal( 0:imax-1, 1:jmax )
real(DP) :: xyz_B0( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: xyz_B1( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: Mu
real(DP) :: Mu1
!!$ real(DP) :: gth( im, jm, km+1 ), pfinth( im, jm, km+1 )
!!$ real(DP) :: b0( ijs:ije, 1, km ), b1( ijs:ije, 1, km )
!!$ real(DP) :: gemis
real(DP) :: xyr_IUw( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: xyr_IDw( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: FactG
real(DP) :: FactH
real(DP) :: FactJ
real(DP) :: FactK
real(DP) :: Alp1
real(DP) :: Alp2
real(DP) :: Sig1
real(DP) :: Sig2
integer :: i, j, k, l
integer :: ms, me
! Variables for debug
!
!!$ real(DP) :: xyr_UwFluxDebug( 0:imax-1, 1:jmax, 0:kmax )
!!$ real(DP) :: xyr_DwFluxDebug( 0:imax-1, 1:jmax, 0:kmax )
! 初期化確認
! Initialization check
!
if ( .not. rad_rte_two_stream_app_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!
! Arguments are checked.
!
select case ( IDScheme )
case ( IDSchemeShortWave )
if ( .not. present( SolarFluxTOA ) ) call MessageNotify( 'E', module_name, 'SolarFluxTOA has to be present.' )
if ( .not. present( xy_InAngle ) ) call MessageNotify( 'E', module_name, 'xy_InAngle has to be present.' )
if ( present( xyr_PFInted ) ) call MessageNotify( 'E', module_name, 'xyr_PFInted need not be present.' )
if ( present( xy_SurfPFInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfPFInted need not be present.' )
if ( present( xy_SurfDPFDTInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted need not be present.' )
if ( present( xyra_DelRadUwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux need not be present.' )
if ( present( xyra_DelRadDwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux need not be present.' )
case ( IDSchemeLongWave )
if ( present( SolarFluxTOA ) ) call MessageNotify( 'E', module_name, 'SolarFluxTOA need not be present.' )
if ( present( xy_InAngle ) ) call MessageNotify( 'E', module_name, 'xy_InAngle need not be present.' )
if ( .not. present( xyr_PFInted ) ) call MessageNotify( 'E', module_name, 'xyr_PFInted has to be present.' )
if ( .not. present( xy_SurfPFInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfPFInted has to be present.' )
if ( .not. present( xy_SurfDPFDTInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted has to be present.' )
if ( .not. present( xyra_DelRadUwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux has to be present.' )
if ( .not. present( xyra_DelRadDwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux has to be present.' )
case default
call MessageNotify( 'E', module_name, 'Unexpected IDScheme, %d', i = (/ IDScheme /) )
end select
do k = 1, kmax
do j = js, je
do i = 0, imax-1
if ( xyz_ssa(i,j,k) >= 1.0d0 ) then
call MessageNotify( 'E', module_name, 'Single Scattering Albedo = %f', d = (/ xyz_ssa(i,j,k) /) )
end if
end do
end do
end do
select case ( IDScheme )
case ( IDSchemeShortWave )
do j = js, je
do i = 0, imax-1
if ( xy_InAngle(i,j) > 0.0d0 ) then
xy_cosSZA (i,j) = 1.0d0 / xy_InAngle(i,j)
xy_cosSZAInv (i,j) = xy_InAngle(i,j)
xy_cosSZAInvsq(i,j) = xy_cosSZAInv(i,j)**2
else
xy_cosSZA (i,j) = 0.0d0
xy_cosSZAInv (i,j) = 0.0d0
xy_cosSZAInvsq(i,j) = 0.0d0
end if
end do
end do
case ( IDSchemeLongWave )
do j = js, je
xy_cosSZA (:,j) = 1.0d100
xy_cosSZAInv (:,j) = 1.0d100
xy_cosSZAInvsq(:,j) = 1.0d100
end do
!!$ gth(:,:,:) = 1.0d100
!!$ else
!!$
!!$ stop 'sw != 1 is not supported.'
!!$
!!$ do ij = ijs, ije
!!$ cosz ( ij, 1 ) = 1.0d100
!!$ cosz2( ij, 1 ) = 1.0d100
!!$ end do
!!$
!!$ do ij = ijs, ije
!!$ gth( ij, 1, 1 ) = gt( ij, 1, 1 )
!!$ end do
!!$ do k = 2, km+1-1
!!$ do ij = ijs, ije
!!$ gth( ij, 1, k ) = ( gt( ij, 1, k-1 ) + gt( ij, 1, k ) ) * 0.5d0
!!$ end do
!!$ end do
!!$ do ij = ijs, ije
!!$ gth( ij, 1, km+1 ) = gts( ij, 1 )
!!$ end do
!!$! call pfint_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$! ijs, ije )
!!$ call pfint_gq_array( wn1, wn2, divnum, im, jm, km+1, gth, pfinth, &
!!$ ijs, ije )
case default
call MessageNotify( 'E', module_name, 'Unexpected IDScheme, %d', i = (/ IDScheme /) )
end select
!
! Delta-Function Adjustment
!
do k = 1, kmax
do j = js, je
xyz_AF2 (:,j,k) = xyz_AF(:,j,k) / ( 1.0d0 + xyz_AF(:,j,k) )
xyz_SSA2(:,j,k) = ( 1.0d0 - xyz_AF(:,j,k)**2 ) * xyz_SSA(:,j,k) / ( 1.0d0 - xyz_SSA(:,j,k) * xyz_AF(:,j,k)**2 )
end do
end do
!
do k = 0, kmax
do j = js, je
xyr_OpDep2(:,j,k) = 0.0d0
end do
end do
do k = kmax-1, 0, -1
do j = js, je
xyr_OpDep2(:,j,k) = xyr_OpDep2(:,j,k+1) + ( 1.0d0 - xyz_SSA(:,j,k+1) * xyz_AF(:,j,k+1)**2 ) * ( xyr_OptDep(:,j,k) - xyr_OptDep(:,j,k+1) )
end do
end do
select case ( IDScheme )
case ( IDSchemeShortWave )
do k = 0, kmax
do j = js, je
xyr_Trans2(:,j,k) = exp( -xyr_OpDep2(:,j,k) * xy_cosSZAInv(:,j) )
end do
end do
!
! Eddington approximation
!
do k = 1, kmax
do j = js, je
xyz_Gam1(:,j,k) = ( 7.0d0 - xyz_SSA2(:,j,k) * ( 4.0d0 + 3.0d0 * xyz_AF2(:,j,k) ) ) / 4.0d0
xyz_Gam2(:,j,k) = -( 1.0d0 - xyz_SSA2(:,j,k) * ( 4.0d0 - 3.0d0 * xyz_AF2(:,j,k) ) ) / 4.0d0
xyz_Gam3(:,j,k) = ( 2.0d0 - 3.0d0 * xyz_AF2(:,j,k) * xy_cosSZA(:,j) ) / 4.0d0
xyz_Gam4(:,j,k) = 1.0d0 - xyz_Gam3(:,j,k)
end do
end do
case ( IDSchemeLongWave )
!
! Treatment if delta-function adjustment is not performed.
!
!!$ do k = 1, kmax
!!$ do j = js, je
!!$ xyz_AF2 (:,j,k) = xyz_AF (:,j,k)
!!$ xyz_SSA2(:,j,k) = xyz_SSA(:,j,k)
!!$ end do
!!$ end do
!!$ do k = 0, kmax
!!$ do j = js, je
!!$ xyr_OpDep2(:,j,k) = xyr_OptDep(:,j,k)
!!$ end do
!!$ end do
do k = 0, kmax
do j = js, je
xyr_Trans2(:,j,k) = 1.0d100
end do
end do
!
! Hemispheric mean approximation
!
do k = 1, kmax
do j = js, je
xyz_Gam1(:,j,k) = 2.0_DP - xyz_SSA2(:,j,k) * ( 1.0_DP + xyz_AF2(:,j,k) )
xyz_Gam2(:,j,k) = xyz_SSA2(:,j,k) * ( 1.0_DP - xyz_AF2(:,j,k) )
xyz_Gam3(:,j,k) = 1.0d100
xyz_Gam4(:,j,k) = 1.0d100
end do
end do
case default
call MessageNotify( 'E', module_name, 'Unexpected IDScheme, %d', i = (/ IDScheme /) )
end select
do k = 1, kmax
do j = js, je
xyz_DelTau(:,j,k) = xyr_OpDep2(:,j,k-1) - xyr_OpDep2(:,j,k)
end do
end do
select case ( IDScheme )
case ( IDSchemeLongWave )
!
! Avoiding singularity when dtau equal to zero
!
do k = 1, kmax
do j = js, je
do i = 0, imax-1
if( xyz_DelTau(i,j,k) < 1.0d-10 ) then
xyz_DelTau(i,j,k) = 0.0d0
end if
end do
end do
end do
end select
do k = 1, kmax
do j = js, je
! In very small parameter space of single scattering albedo and asymmetry
! factor close to asymmetry factor of 1, xyz_Gam1**2 - xyz_Gam2**2 becomes
! negative value. (yot, 2011/11/20)
!!$ xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) )
!!$ xyz_Lambda(:,j,k) = sqrt( max( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k), 0.0_DP ) )
xyz_Lambda(:,j,k) = sqrt( xyz_Gam1(:,j,k) * xyz_Gam1(:,j,k) - xyz_Gam2(:,j,k) * xyz_Gam2(:,j,k) + 1.0d-10 )
xyz_LGamma(:,j,k) = xyz_Gam2(:,j,k) / ( xyz_Gam1(:,j,k) + xyz_Lambda(:,j,k) )
end do
end do
do k = 1, kmax
do j = js, je
xy_TMPVal(:,j) = exp( - xyz_Lambda(:,j,k) * xyz_DelTau(:,j,k) )
xyaz_smalle(:,j,1,k) = xyz_LGamma(:,j,k) * xy_TMPVal(:,j) + 1.0_DP
xyaz_smalle(:,j,2,k) = xyz_LGamma(:,j,k) * xy_TMPVal(:,j) - 1.0_DP
xyaz_smalle(:,j,3,k) = xy_TMPVal(:,j) + xyz_LGamma(:,j,k)
xyaz_smalle(:,j,4,k) = xy_TMPVal(:,j) - xyz_LGamma(:,j,k)
end do
end do
select case ( IDScheme )
case ( IDSchemeShortWave )
do k = 1, kmax
do j = js, je
! Lines below will be deleted in future (yot, 2010/09/14).
!!$ xyz_CUpB(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1) &
!!$ & * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
!!$ xyz_CUpT(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k ) &
!!$ & * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
!!$ xyz_CDoB(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k-1) &
!!$ & * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
!!$ xyz_CDoT(:,j,k) = &
!!$ & xyz_SSA2(:,j,k) * SolarFluxTOA * xyr_Trans2(:,j,k ) &
!!$ & * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) &
!!$ & + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) &
!!$ & / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
xy_TMPVal(:,j) = xyz_SSA2(:,j,k) * SolarFluxTOA * ( ( xyz_Gam1(:,j,k) - xy_cosSZAInv(:,j) ) * xyz_Gam3(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam4(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
xyz_CUpB(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k-1)
xyz_CUpT(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k )
xy_TMPVal(:,j) = xyz_SSA2(:,j,k) * SolarFluxTOA * ( ( xyz_Gam1(:,j,k) + xy_cosSZAInv(:,j) ) * xyz_Gam4(:,j,k) + xyz_Gam2(:,j,k) * xyz_Gam3(:,j,k) ) / ( xyz_Lambda(:,j,k) * xyz_Lambda(:,j,k) - xy_cosSZAInvsq(:,j) )
xyz_CDoB(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k-1)
xyz_CDoT(:,j,k) = xy_TMPVal(:,j) * xyr_Trans2(:,j,k )
end do
end do
case ( IDSchemeLongWave )
do k = 1, kmax
do j = js, je
do i = 0, imax-1
!
! Notice!
! Avoiding singularity when dtau equal to zero.
! dtau occationally has much smaller value.
! When this occurs, b1 cannot be calculated correctly.
!
if( xyz_DelTau(i,j,k) /= 0.0_DP ) then
xyz_B0(i,j,k) = xyr_PFInted(i,j,k)
xyz_B1(i,j,k) = ( xyr_PFInted(i,j,k-1) - xyr_PFInted(i,j,k) ) / xyz_DelTau(i,j,k)
else
xyz_B0(i,j,k) = 0.0_DP
xyz_B1(i,j,k) = 0.0_DP
end if
end do
end do
end do
do k = 1, kmax
do j = js, je
do i = 0, imax-1
Mu1 = ( 1.0_DP - xyz_SSA2(i,j,k) ) / ( xyz_Gam1(i,j,k) - xyz_Gam2(i,j,k) )
!!$ xyz_CUpB(i,j,k) = 2.0_DP * PI * Mu1 &
xyz_CUpB(i,j,k) = 2.0_DP * Mu1 * ( xyz_B0(i,j,k) + xyz_B1(i,j,k) * ( xyz_DelTau(i,j,k) + 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
!!$ xyz_CUpT(i,j,k) = 2.0_DP * PI * Mu1 &
xyz_CUpT(i,j,k) = 2.0_DP * Mu1 * ( xyz_B0(i,j,k) + xyz_B1(i,j,k) * ( 0.0_DP + 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
!!$ xyz_CDoB(i,j,k) = 2.0_DP * PI * Mu1 &
xyz_CDoB(i,j,k) = 2.0_DP * Mu1 * ( xyz_B0(i,j,k) + xyz_B1(i,j,k) * ( xyz_DelTau(i,j,k) - 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
!!$ xyz_CDoT(i,j,k) = 2.0_DP * PI * Mu1 &
xyz_CDoT(i,j,k) = 2.0_DP * Mu1 * ( xyz_B0(i,j,k) + xyz_B1(i,j,k) * ( 0.0_DP - 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) ) )
end do
end do
end do
case default
call MessageNotify( 'E', module_name, 'Unexpected IDScheme, %d', i = (/ IDScheme /) )
end select
select case( IDScheme )
case ( IDSchemeShortWave )
do j = js, je
do i = 0, imax-1
xy_SurfSrc(:,j) = xy_SurfAlbedo(:,j) * SolarFluxTOA * xyr_Trans2(:,j,0) * xy_cosSZA(:,j)
end do
end do
case ( IDSchemeLongWave )
do j = js, je
do i = 0, imax-1
xy_SurfSrc(:,j) = xy_SurfPFInted(:,j)
end do
end do
!!$ else if( sw .eq. 2 ) then
!!$ do ij = ijs, ije
!!$! gemis = 1.0d0
!!$ gemis = emis( ij, 1 )
!!$ ea( (ij-ijs+1), l ) &
!!$ = -CUpB( ij, 1, km ) + ( 1.0d0 - gemis ) * CDoB( ij, 1, km ) &
!!$ + gemis * pi * pfinth( ij, 1, km+1 )
!!$ end do
case default
call MessageNotify( 'E', module_name, 'Unexpected IDScheme, %d', i = (/ IDScheme /) )
end select
k = 1
l = 1
do j = js, je
do i = 0, imax-1
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = 0.0_DP
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,3,k)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = - ( xyaz_smalle(i,j,2,k) - xy_SurfAlbedo(i,j) * xyaz_smalle(i,j,4,k) )
end do
end do
do j = js, je
do i = 0, imax-1
aa_Vec( i+imax*(j-1)+1, l ) = - xyz_CUpB(i,j,k) + xy_SurfAlbedo(i,j) * xyz_CDoB(i,j,k) + xy_SurfSrc(i,j)
end do
end do
do k = 1, kmax-1
do j = js, je
do i = 0, imax-1
l = 2 * k ! equation number
!
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,2,k+1) - xyaz_smalle(i,j,3,k ) * xyaz_smalle(i,j,4,k+1)
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k ) * xyaz_smalle(i,j,2,k+1) - xyaz_smalle(i,j,4,k ) * xyaz_smalle(i,j,4,k+1)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k+1) * xyaz_smalle(i,j,4,k+1) - xyaz_smalle(i,j,2,k+1) * xyaz_smalle(i,j,3,k+1)
aa_Vec ( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k+1) * ( - xyz_CDoT(i,j,k ) + xyz_CDoB(i,j,k+1) ) - xyaz_smalle(i,j,4,k+1) * ( - xyz_CUpT(i,j,k ) + xyz_CUpB(i,j,k+1) )
l = 2 * k + 1 ! equation number
!
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k ) * xyaz_smalle(i,j,3,k ) - xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,4,k )
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,1,k+1) - xyaz_smalle(i,j,3,k ) * xyaz_smalle(i,j,3,k+1)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,3,k ) * xyaz_smalle(i,j,4,k+1) - xyaz_smalle(i,j,1,k ) * xyaz_smalle(i,j,2,k+1)
aa_Vec ( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,3,k ) * ( - xyz_CDoT(i,j,k ) + xyz_CDoB(i,j,k+1) ) - xyaz_smalle(i,j,1,k ) * ( - xyz_CUpT(i,j,k ) + xyz_CUpB(i,j,k+1) )
end do
end do
end do
k = kmax
l = 2 * kmax
do j = js, je
do i = 0, imax-1
aa_TridiagMtx1( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,1,k)
aa_TridiagMtx2( i+imax*(j-1)+1, l ) = xyaz_smalle(i,j,2,k)
aa_TridiagMtx3( i+imax*(j-1)+1, l ) = 0.0d0
aa_Vec ( i+imax*(j-1)+1, l ) = -xyz_CDoT(i,j,k) + 0.0d0
end do
end do
ms = 0 + imax*(js-1)+1
me = imax-1 + imax*(je-1)+1
call tridiag( imax*jmax, 2*kmax, aa_TridiagMtx1, aa_TridiagMtx2, aa_TridiagMtx3, aa_Vec, ms, me )
select case ( IDScheme )
case ( IDSchemeShortWave )
k = 1
do j = js, je
do i = 0, imax-1
xyr_RadUwFlux(i,j,k-1) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) + xyz_CUpB(i,j,k)
xyr_RadDwFlux(i,j,k-1) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) + xyz_CDoB(i,j,k)
end do
end do
do k = 1, kmax
do j = js, je
do i = 0, imax-1
xyr_RadUwFlux(i,j,k) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) + xyz_CUpT(i,j,k)
xyr_RadDwFlux(i,j,k) = aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) + xyz_CDoT(i,j,k)
end do
end do
end do
! Code for debug
!
!!$ do k = 1, kmax
!!$ do j = js, je
!!$ do i = 0, imax-1
!!$ xyr_UwFluxDebug(i,j,k-1) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
!!$ & - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) &
!!$ & + xyz_CUpB(i,j,k)
!!$ xyr_DwFluxDebug(i,j,k-1) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
!!$ & - aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) &
!!$ & + xyz_CDoB(i,j,k)
!!$ end do
!!$ end do
!!$ end do
!!$
!!$ k = kmax
!!$ do j = js, je
!!$ do i = 0, imax-1
!!$ xyr_UwFluxDebug(i,j,k) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,3,k) &
!!$ & + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,4,k) &
!!$ & + xyz_CUpT(i,j,k)
!!$ xyr_DwFluxDebug(i,j,k) = &
!!$ & aa_Vec( i+imax*(j-1)+1, 2*k-1) * xyaz_smalle(i,j,1,k) &
!!$ & + aa_Vec( i+imax*(j-1)+1, 2*k ) * xyaz_smalle(i,j,2,k) &
!!$ & + xyz_CDoT(i,j,k)
!!$ end do
!!$ end do
!!$
!!$
!!$ i = imax/2
!!$ j = jmax/2
!!$ do k = kmax, 0, -1
!!$ write( 6, * ) k, xyr_UwFlux(i,j,k), xyr_UwFluxDebug(i,j,k), xyr_UwFlux(i,j,k) - xyr_UwFluxDebug(i,j,k)
!!$ end do
!!$ do k = kmax, 0, -1
!!$ write( 6, * ) k, xyr_DwFlux(i,j,k), xyr_DwFluxDebug(i,j,k), xyr_DwFlux(i,j,k) - xyr_DwFluxDebug(i,j,k)
!!$ end do
!!$ k = 0
!!$ write( 6, * ) k, xyr_UwFlux(i,j,k), xy_SurfAlbedo(i,j) * xyr_DwFlux(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_Trans2(i,j,0) * xy_cosSZA(i,j), &
!!$ xyr_UwFlux(i,j,k) - ( xy_SurfAlbedo(i,j) * xyr_DwFlux(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * xyr_Trans2(i,j,0) * xy_cosSZA(i,j) )
!!$ stop
!
! Addition of Direct Solar Insident
!
do k = 0, kmax
do j = js, je
xyr_RadDwFlux(:,j,k) = xyr_RadDwFlux(:,j,k) + SolarFluxTOA * xyr_Trans2(:,j,k) * xy_cosSZA(:,j)
end do
end do
!
! Set zero at nightside
!
do k = 0, kmax
do j = js, je
do i = 0, imax-1
if( xy_cosSZA(i,j) <= 0.0d0 ) then
xyr_RadUwFlux(i,j,k) = 0.0d0
xyr_RadDwFlux(i,j,k) = 0.0d0
end if
end do
end do
end do
case ( IDSchemeLongWave )
!!$ else if( sw .eq. 2 ) then
! Source function technique described by Toon et al. [1989]
! is used to calculated infrared heating rate.
!
do k = 0, kmax
do j = js, je
xyr_RadUwFlux(:,j,k) = 0.0_DP
xyr_RadDwFlux(:,j,k) = 0.0_DP
end do
end do
do k = 0, kmax
do j = js, je
xyra_DelRadUwFlux(:,j,k,0) = 0.0_DP
xyra_DelRadUwFlux(:,j,k,1) = 0.0_DP
xyra_DelRadDwFlux(:,j,k,0) = 0.0_DP
xyra_DelRadDwFlux(:,j,k,1) = 0.0_DP
end do
end do
do l = 1, NGaussQuad
Mu = a_GQP( l )
k = kmax
do j = js, je
xyr_IDw(:,j,k) = 0.0_DP
end do
do k = kmax, 1, -1
do j = js, je
do i = 0, imax-1
Mu1 = ( 1.0_DP - xyz_SSA2(i,j,k) ) / ( xyz_Gam1(i,j,k) - xyz_Gam2(i,j,k) )
FactJ = ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) + aa_Vec( i+imax*(j-1)+1, 2*k ) ) * xyz_LGamma(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0_DP / Mu1 )
FactK = ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) - aa_Vec( i+imax*(j-1)+1, 2*k ) ) * ( 1.0_DP / Mu1 - xyz_Lambda(i,j,k) )
Sig1 = 2.0_DP * ( xyz_B0(i,j,k) - xyz_B1(i,j,k) * ( 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) - Mu1 ) )
Sig2 = 2.0_DP * xyz_B1(i,j,k)
xyr_IDw(i,j,k-1) = xyr_IDw(i,j,k) * exp( - xyz_DelTau(i,j,k) / Mu ) + FactJ / ( xyz_Lambda(i,j,k) * Mu + 1.0_DP ) * ( 1.0_DP - exp( - xyz_DelTau(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0d0 / Mu ) ) ) + FactK / ( xyz_Lambda(i,j,k) * Mu - 1.0_DP ) * ( exp( - xyz_DelTau(i,j,k) / Mu ) - exp( - xyz_DelTau(i,j,k) * xyz_Lambda(i,j,k) ) ) + Sig1 * ( 1.0_DP - exp( - xyz_DelTau(i,j,k) / Mu ) ) + Sig2 * ( Mu * exp( - xyz_DelTau(i,j,k) / Mu ) + xyz_DelTau(i,j,k) - Mu )
end do
end do
end do
k = 0
do j = js, je
! gemis = 1.0d0
!!$ gemis = emis( ij, 1 )
!!$ inteup( ij, 1, k ) = ( 1.0d0 - gemis ) * intedo( ij, 1, km+1 ) &
!!$ + gemis * pix2 * pfinth( ij, 1, km+1 )
xyr_IUw(:,j,k) = xy_SurfAlbedo(:,j) * xyr_IDw(:,j,0) + 2.0_DP * xy_SurfPFInted(:,j)
end do
do k = 1, kmax
do j = js, je
do i = 0, imax-1
FactG = ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) + aa_Vec( i+imax*(j-1)+1, 2*k ) ) * ( 1.0_DP / Mu1 - xyz_Lambda(i,j,k) )
FactH = ( aa_Vec( i+imax*(j-1)+1, 2*k-1 ) - aa_Vec( i+imax*(j-1)+1, 2*k ) ) * xyz_LGamma(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0_DP / Mu1 )
Alp1 = 2.0_DP * ( xyz_B0(i,j,k) + xyz_B1(i,j,k) * ( 1.0_DP / ( xyz_Gam1(i,j,k) + xyz_Gam2(i,j,k) ) - Mu1 ) )
Alp2 = 2.0_DP * xyz_B1(i,j,k)
xyr_IUw(i,j,k) = xyr_IUw(i,j,k-1) * exp( - xyz_DelTau(i,j,k) / Mu ) + FactG / ( xyz_Lambda(i,j,k) * Mu - 1.0_DP ) * ( exp( - xyz_DelTau(i,j,k) / Mu ) - exp( - xyz_DelTau(i,j,k) * xyz_Lambda(i,j,k) ) ) + FactH / ( xyz_Lambda(i,j,k) * Mu + 1.0_DP ) * ( 1.0_DP - exp( - xyz_DelTau(i,j,k) * ( xyz_Lambda(i,j,k) + 1.0d0 / Mu ) ) ) + Alp1 * ( 1.0_DP - exp( - xyz_DelTau(i,j,k) / Mu ) ) + Alp2 * ( Mu - ( xyz_DelTau(i,j,k) + Mu ) * exp( - xyz_Deltau(i,j,k) / Mu ) )
end do
end do
end do
do k = 0, kmax
do j = js, je
xyr_RadUwFlux(:,j,k) = xyr_RadUwFlux(:,j,k) + Mu * xyr_IUw(:,j,k) * a_GQW( l )
xyr_RadDwFlux(:,j,k) = xyr_RadDwFlux(:,j,k) + Mu * xyr_IDw(:,j,k) * a_GQW( l )
end do
end do
do k = 0, kmax
do j = js, je
xyra_DelRadUwFlux(:,j,k,0) = xyra_DelRadUwFlux(:,j,k,0) + Mu * 2.0_DP * xy_SurfDPFDTInted(:,j) * exp( - ( xyr_OpDep2(:,j,0) - xyr_OpDep2(:,j,k) ) / Mu ) * a_GQW( l )
end do
end do
end do ! do l = 1, NGaussQuad
case default
call MessageNotify( 'E', module_name, 'Unexpected IDScheme, %d', i = (/ IDScheme /) )
end select
end subroutine RadRTETwoStreamAppCore
| Subroutine : | |||
| IDScheme : | integer , intent(in ) | ||
| xyz_SSA( 0:imax-1, 1:jmax, 1:kmax ) : | real(DP), intent(in ) | ||
| xyz_AF( 0:imax-1, 1:jmax, 1:kmax ) : | real(DP), intent(in ) | ||
| xyr_OptDep(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ) | ||
| xy_SurfAlbedo( 0:imax-1, 1:jmax ) : | real(DP), intent(in ) | ||
| xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) | ||
| xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(out) | ||
| SolarFluxTOA : | real(DP), intent(in ), optional | ||
| xy_InAngle(0:imax-1, 1:jmax) : | real(DP), intent(in ), optional
| ||
| xyr_PFInted(0:imax-1, 1:jmax, 0:kmax) : | real(DP), intent(in ), optional | ||
| xy_SurfPFInted(0:imax-1, 1:jmax) : | real(DP), intent(in ), optional | ||
| xy_SurfDPFDTInted(0:imax-1, 1:jmax) : | real(DP), intent(in ), optional | ||
| xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(out), optional | ||
| xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) : | real(DP), intent(out), optional |
subroutine RadRTETwoStreamAppWrapper( IDScheme, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, SolarFluxTOA, xy_InAngle, xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted, xyra_DelRadUwFlux, xyra_DelRadDwFlux )
! USE statements
!
! OpenMP
!
!$ use omp_lib
integer , intent(in ) :: IDScheme
real(DP), intent(in ) :: xyz_SSA ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyz_AF ( 0:imax-1, 1:jmax, 1:kmax )
real(DP), intent(in ) :: xyr_OptDep (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ) :: xy_SurfAlbedo ( 0:imax-1, 1:jmax )
real(DP), intent(out) :: xyr_RadUwFlux(0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(out) :: xyr_RadDwFlux(0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ), optional :: SolarFluxTOA
real(DP), intent(in ), optional :: xy_InAngle (0:imax-1, 1:jmax)
! sec (入射角).
! sec (angle of incidence)
real(DP), intent(in ), optional :: xyr_PFInted (0:imax-1, 1:jmax, 0:kmax)
real(DP), intent(in ), optional :: xy_SurfPFInted (0:imax-1, 1:jmax)
real(DP), intent(in ), optional :: xy_SurfDPFDTInted(0:imax-1, 1:jmax)
real(DP), intent(out), optional :: xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
real(DP), intent(out), optional :: xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
! Local variables
!
integer :: js
integer :: je
integer :: nthreads
integer, allocatable :: a_js(:)
integer, allocatable :: a_je(:)
integer :: n
! 初期化確認
! Initialization check
!
if ( .not. rad_rte_two_stream_app_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
!
! Arguments are checked.
!
select case ( IDScheme )
case ( IDSchemeShortWave )
if ( .not. present( SolarFluxTOA ) ) call MessageNotify( 'E', module_name, 'SolarFluxTOA has to be present.' )
if ( .not. present( xy_InAngle ) ) call MessageNotify( 'E', module_name, 'xy_InAngle has to be present.' )
if ( present( xyr_PFInted ) ) call MessageNotify( 'E', module_name, 'xyr_PFInted need not be present.' )
if ( present( xy_SurfPFInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfPFInted need not be present.' )
if ( present( xy_SurfDPFDTInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted need not be present.' )
if ( present( xyra_DelRadUwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux need not be present.' )
if ( present( xyra_DelRadDwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux need not be present.' )
case ( IDSchemeLongWave )
if ( present( SolarFluxTOA ) ) call MessageNotify( 'E', module_name, 'SolarFluxTOA need not be present.' )
if ( present( xy_InAngle ) ) call MessageNotify( 'E', module_name, 'xy_InAngle need not be present.' )
if ( .not. present( xyr_PFInted ) ) call MessageNotify( 'E', module_name, 'xyr_PFInted has to be present.' )
if ( .not. present( xy_SurfPFInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfPFInted has to be present.' )
if ( .not. present( xy_SurfDPFDTInted ) ) call MessageNotify( 'E', module_name, 'xy_SurfDPFDTInted has to be present.' )
if ( .not. present( xyra_DelRadUwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadUwFlux has to be present.' )
if ( .not. present( xyra_DelRadDwFlux ) ) call MessageNotify( 'E', module_name, 'xyra_DelRadLwFlux has to be present.' )
case default
call MessageNotify( 'E', module_name, 'Unexpected IDScheme, %d', i = (/ IDScheme /) )
end select
nthreads = 1
!$ nthreads = omp_get_max_threads()
!!$ !$ write( 6, * ) "Number of processors : ", omp_get_num_procs()
!!$ !$ write( 6, * ) "Number of threads : ", nthreads
allocate( a_js(0:nthreads-1) )
allocate( a_je(0:nthreads-1) )
do n = 0, nthreads-1
if ( n == 0 ) then
a_js(n) = 1
else
a_js(n) = a_je(n-1) + 1
end if
a_je(n) = a_js(n ) + jmax / nthreads - 1
if ( n + 1 <= mod( jmax, nthreads ) ) then
a_je(n) = a_je(n) + 1
end if
end do
!$OMP PARALLEL DEFAULT(PRIVATE) &
!$OMP SHARED(nthreads,a_js,a_je, &
!$OMP xyz_SSA, xyz_AF, &
!$OMP SolarFluxTOA, &
!$OMP xy_SurfAlbedo, &
!$OMP IDScheme, &
!$OMP xy_InAngle, &
!$OMP xyr_OptDep, &
!$OMP xyr_RadUwFlux, &
!$OMP xyr_RadDwFlux, &
!$OMP xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted, &
!$OMP xyra_DelRadUwFlux, xyra_DelRadDwFlux)
!$OMP DO
do n = 0, nthreads-1
js = a_js(n)
je = a_je(n)
if ( js > je ) cycle
!!$ write( 6, * ) n, js, je
call RadRTETwoStreamAppCore( IDScheme, xyz_SSA, xyz_AF, xyr_OptDep, xy_SurfAlbedo, xyr_RadUwFlux, xyr_RadDwFlux, js, je, SolarFluxTOA = SolarFluxTOA, xy_InAngle = xy_InAngle, xyr_PFInted = xyr_PFInted, xy_SurfPFInted = xy_SurfPFInted, xy_SurfDPFDTInted = xy_SurfDPFDTInted, xyra_DelRadUwFlux = xyra_DelRadUwFlux, xyra_DelRadDwFlux = xyra_DelRadDwFlux )
end do
!$OMP END DO
!$OMP END PARALLEL
deallocate( a_js )
deallocate( a_je )
end subroutine RadRTETwoStreamAppWrapper
| Constant : | |||
| module_name = ‘rad_rte_two_stream_app‘ : | character(*), parameter
|
| Variable : | |||
| rad_rte_two_stream_app_inited = .false. : | logical, save
|
| Subroutine : | |
| mm : | integer , intent(in ) |
| jmx : | integer , intent(in ) |
| a( mm, jmx ) : | real(DP), intent(in ) |
| b( mm, jmx ) : | real(DP), intent(in ) |
| c( mm, jmx ) : | real(DP), intent(in ) |
| f( mm, jmx ) : | real(DP), intent(inout) |
| ms : | integer , intent(in ) |
| me : | integer , intent(in ) |
subroutine tridiag( mm, jmx, a, b, c, f, ms, me )
integer , intent(in ) :: mm, jmx
real(DP), intent(in ) :: a( mm, jmx ),b( mm, jmx )
real(DP), intent(in ) :: c( mm, jmx )
real(DP), intent(inout) :: f( mm, jmx )
integer , intent(in ) :: ms, me
! Local variables
!
real(DP) :: q( mm, jmx ), p
integer :: j, m
! Forward elimination sweep
!
do m = ms, me
q( m, 1 ) = - c( m, 1 ) / b( m, 1 )
f( m, 1 ) = f( m, 1 ) / b( m, 1 )
end do
do j = 2, jmx
do m = ms, me
p = 1.0d0 / ( b( m, j ) + a( m, j ) * q( m, j-1 ) )
q( m, j ) = - c( m, j ) * p
f( m, j ) = ( f( m, j ) - a( m, j ) * f( m, j-1 ) ) * p
end do
end do
! Backward pass
!
do j = jmx - 1, 1, -1
do m = ms, me
f( m, j ) = f( m, j ) + q( m, j ) * f( m, j+1 )
end do
end do
end subroutine tridiag
| Subroutine : | |
| jmx : | integer , intent(in ) |
| a(jmx) : | real(DP), intent(in ) |
| b(jmx) : | real(DP), intent(in ) |
| c(jmx) : | real(DP), intent(in ) |
| f(jmx) : | real(DP), intent(inout) |
subroutine tridiag1( jmx, a, b, c, f )
integer , intent(in ) :: jmx
real(DP), intent(in ) :: a(jmx),b(jmx)
real(DP), intent(in ) :: c(jmx)
real(DP), intent(inout) :: f(jmx)
! Local variables
!
real(DP) :: q(jmx), p
integer :: j
! Forward elimination sweep
!
q( 1 ) = - c( 1 ) / b( 1 )
f( 1 ) = f( 1 ) / b( 1 )
do j = 2, jmx
p = 1.0d0 / ( b( j ) + a( j ) * q( j-1 ) )
q( j ) = - c( j ) * p
f( j ) = ( f( j ) - a( j ) * f( j-1 ) ) * p
end do
! Backward pass
!
do j = jmx - 1, 1, -1
f( j ) = f( j ) + q( j ) * f( j+1 )
end do
end subroutine tridiag1
| Constant : | |||
| version = ’$Name: dcpam5-20120813 $’ // ’$Id: rad_rte_two_stream_app.f90,v 1.4 2012-05-08 15:12:12 yot Exp $’ : | character(*), parameter
|