subroutine RadiationTwoStreamApp( QeRatio, xyz_ssa, xyz_af, SolarFluxTOA, xy_SurfAlbedo, IDScheme, xy_InAngle, xyr_OptDepBase, xyr_RadFlux )
! USE statements
!
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
!
! Physical constants settings
!
use constants, only: Grav, CpDry, PI, StB ! $ \sigma_{SB} $ .
! ステファンボルツマン定数.
! Stefan-Boltzmann constant
! 格子点設定
! Grid points settings
!
use gridset, only: imax, jmax, kmax ! 鉛直層数.
! Number of vertical level
!!$ use pf_module , only : pfint_gq_array
real(DP), intent(in ) :: QeRatio
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_OptDepBase( 0:imax-1, 1:jmax, 0:kmax )
real(DP), intent(out) :: xyr_RadFlux ( 0:imax-1, 1:jmax, 0:kmax )
!!$ 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 )
!
! opdep : Optical depth
!
real(DP) :: opdep( 0:imax-1, 1:jmax, 0:kmax )
!
! 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) :: opdep2 ( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: exp_opdep2( 0:imax-1, 1:jmax, 0:kmax )
!
! gam? : Coefficients of Generalized Radiative Transfer Equation
!
real(DP) :: gam1( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: gam2( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: gam3( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: gam4( 0:imax-1, 1:jmax, 1:kmax )
!
! upfl : Upward Flux
! dofl : Downward Flux
!
real(DP) :: xyr_UwFlux( 0:imax-1, 1:jmax, 0:kmax )
real(DP) :: xyr_DwFlux( 0:imax-1, 1:jmax, 0: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) :: dtau( 0:imax-1, 1:jmax, 1:kmax )
integer :: i, j, k, l
!!$ integer :: m
real(DP) :: lambda( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: lsigma( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: 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) :: cupb( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: cupt( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: cdob( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: cdot( 0:imax-1, 1:jmax, 1:kmax )
real(DP) :: aa( imax*jmax, kmax*2 )
real(DP) :: ba( imax*jmax, kmax*2 )
real(DP) :: da( imax*jmax, kmax*2 )
real(DP) :: ea( imax*jmax, kmax*2 )
real(DP) :: 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 )
!!$ call gauleg( 0.0d0, 1.0d0, quadn, qucos, qudcos )
do k = 1, kmax
do j = 1, jmax
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 = 1, jmax
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
!
! calculation of Dust Optical Depth
!
opdep(:,:,:) = xyr_OptDepBase * QeRatio
if( IDScheme .eq. IDSchemeShortWave ) then
!
! Delta-Function Adjustment
!
xyz_af2 = xyz_af / ( 1.0d0 + xyz_af )
xyz_ssa2 = ( 1.0d0 - xyz_af**2 ) * xyz_ssa / ( 1.0d0 - xyz_ssa * xyz_af**2 )
!!$ opdep2(:,:,:) = ( 1.0d0 - xyz_ssa * xyz_af**2 ) * opdep(:,:,:)
opdep2(:,:,:) = 0.0d0
do k = kmax-1, 0, -1
opdep2(:,:,k) = opdep2(:,:,k+1) + ( 1.0d0 - xyz_ssa(:,:,k+1) * xyz_af(:,:,k+1)**2 ) * ( opdep(:,:,k) - opdep(:,:,k+1) )
end do
do k = 0, kmax
exp_opdep2(:,:,k) = exp( -opdep2(:,:,k) * xy_cosSZAInv )
end do
!
! Eddington approximation
!
do k = 1, kmax
gam1(:,:,k) = ( 7.0d0 - xyz_ssa2(:,:,k) * ( 4.0d0 + 3.0d0 * xyz_af2(:,:,k) ) ) / 4.0d0
gam2(:,:,k) = -( 1.0d0 - xyz_ssa2(:,:,k) * ( 4.0d0 - 3.0d0 * xyz_af2(:,:,k) ) ) / 4.0d0
gam3(:,:,k) = ( 2.0d0 - 3.0d0 * xyz_af2(:,:,k) * xy_cosSZA ) / 4.0d0
gam4(:,:,k) = 1.0d0 - gam3(:,:,k)
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
dtau(:,:,k) = opdep2(:,:,k-1) - opdep2(:,:,k)
end do
!
! Avoiding singularity when dtau equal to zero
!
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if( ( IDScheme .eq. IDSchemeLongWave ) .and. ( dtau(i,j,k) .lt. 1.0d-10 ) ) then
dtau(i,j,k) = 0.0d0
end if
end do
end do
end do
lambda(:,:,:) = sqrt( gam1(:,:,:) * gam1(:,:,:) - gam2(:,:,:) * gam2(:,:,:) )
lsigma(:,:,:) = gam2(:,:,:) / ( gam1(:,:,:) + lambda(:,:,:) )
do k = 1, kmax
tmpval(:,:) = exp( - lambda(:,:,k) * dtau(:,:,k) )
smalle(:,:,1,k) = 1.0d0 + lsigma(:,:,k) * tmpval(:,:)
smalle(:,:,2,k) = 1.0d0 - lsigma(:,:,k) * tmpval(:,:)
smalle(:,:,3,k) = lsigma(:,:,k) + tmpval(:,:)
smalle(:,:,4,k) = lsigma(:,:,k) - tmpval(:,:)
end do
if( IDScheme .eq. IDSchemeShortWave ) then
do k = 1, kmax
cupb(:,:,k) = xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k-1) * ( ( gam1(:,:,k) - xy_cosSZAInv ) * gam3(:,:,k) + gam2(:,:,k) * gam4(:,:,k) ) / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
cupt(:,:,k) = xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k ) * ( ( gam1(:,:,k) - xy_cosSZAInv ) * gam3(:,:,k) + gam2(:,:,k) * gam4(:,:,k) ) / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
cdob(:,:,k) = xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k-1) * ( ( gam1(:,:,k) + xy_cosSZAInv ) * gam4(:,:,k) + gam2(:,:,k) * gam3(:,:,k) ) / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
cdot(:,:,k) = xyz_ssa2(:,:,k) * SolarFluxTOA * exp_opdep2(:,:,k ) * ( ( gam1(:,:,k) + xy_cosSZAInv ) * gam4(:,:,k) + gam2(:,:,k) * gam3(:,:,k) ) / ( lambda(:,:,k) * lambda(:,:,k) - xy_cosSZAInvsq )
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 = 1, jmax
do i = 0, imax-1
aa( i+imax*(j-1)+1, l ) = 0.0d0
ba( i+imax*(j-1)+1, l ) = smalle(i,j,2,k) - xy_SurfAlbedo(i,j) * smalle(i,j,4,k)
da( i+imax*(j-1)+1, l ) = smalle(i,j,1,k) - xy_SurfAlbedo(i,j) * smalle(i,j,3,k)
end do
end do
if( IDScheme .eq. IDSchemeShortWave ) then
do j = 1, jmax
do i = 0, imax-1
ea( i+imax*(j-1)+1, l ) = - cupb(i,j,k) + xy_SurfAlbedo(i,j) * cdob(i,j,k) + xy_SurfAlbedo(i,j) * SolarFluxTOA * exp_opdep2(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+1, kmax
!!$
!!$ do j = 1, jmax
!!$ do i = 0, imax-1
!!$
!!$ ! for even l
!!$ !
!!$ l = 2 * (k-1)
!!$ !
!!$ aa( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,3,k) * smalle(i,j,4,k-1) &
!!$ & - smalle(i,j,1,k) * smalle(i,j,2,k-1)
!!$ ba( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,1,k) * smalle(i,j,1,k-1) &
!!$ & - smalle(i,j,3,k) * smalle(i,j,3,k-1)
!!$ da( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,2,k) * smalle(i,j,3,k ) &
!!$ & - smalle(i,j,1,k) * smalle(i,j,4,k )
!!$ ea( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,3,k) * ( cupt(i,j,k-1) - cupb(i,j,k) ) &
!!$ & - smalle(i,j,1,k) * ( cdot(i,j,k-1) - cdob(i,j,k) )
!!$
!!$ ! for odd l
!!$ !
!!$ l = 2 * (k-1) + 1
!!$
!!$ aa( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,1,k-1) * smalle(i,j,4,k-1) &
!!$ & - smalle(i,j,2,k-1) * smalle(i,j,3,k-1)
!!$ ba( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,2,k ) * smalle(i,j,2,k-1) &
!!$ & - smalle(i,j,4,k ) * smalle(i,j,4,k-1)
!!$ da( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,1,k ) * smalle(i,j,2,k-1) &
!!$ & - smalle(i,j,3,k ) * smalle(i,j,4,k-1)
!!$ ea( i+imax*(j-1), l ) = &
!!$ & smalle(i,j,2,k-1) * ( cupt(i,j,k-1) - cupb(i,j,k) ) &
!!$ & - smalle(i,j,4,k-1) * ( cdot(i,j,k-1) - cdob(i,j,k) )
!!$ end do
!!$ end do
!!$ end do
do k = 1, kmax-1
do j = 1, jmax
do i = 0, imax-1
l = 2 * k ! equation number
!
aa( i+imax*(j-1)+1, l ) = smalle(i,j,3,k+1) * smalle(i,j,4,k ) - smalle(i,j,1,k+1) * smalle(i,j,2,k )
ba( i+imax*(j-1)+1, l ) = smalle(i,j,1,k+1) * smalle(i,j,1,k ) - smalle(i,j,3,k+1) * smalle(i,j,3,k )
da( i+imax*(j-1)+1, l ) = smalle(i,j,2,k+1) * smalle(i,j,3,k+1) - smalle(i,j,1,k+1) * smalle(i,j,4,k+1)
ea( i+imax*(j-1)+1, l ) = smalle(i,j,3,k+1) * ( cupt(i,j,k ) - cupb(i,j,k+1) ) - smalle(i,j,1,k+1) * ( cdot(i,j,k ) - cdob(i,j,k+1) )
l = 2 * k + 1 ! equation number
!
aa( i+imax*(j-1)+1, l ) = smalle(i,j,1,k ) * smalle(i,j,4,k ) - smalle(i,j,2,k ) * smalle(i,j,3,k )
ba( i+imax*(j-1)+1, l ) = smalle(i,j,2,k+1) * smalle(i,j,2,k ) - smalle(i,j,4,k+1) * smalle(i,j,4,k )
da( i+imax*(j-1)+1, l ) = smalle(i,j,1,k+1) * smalle(i,j,2,k ) - smalle(i,j,3,k+1) * smalle(i,j,4,k )
ea( i+imax*(j-1)+1, l ) = smalle(i,j,2,k ) * ( cupt(i,j,k ) - cupb(i,j,k+1) ) - smalle(i,j,4,k ) * ( cdot(i,j,k ) - cdob(i,j,k+1) )
end do
end do
end do
k = kmax
l = 2 * kmax
do j = 1, jmax
do i = 0, imax-1
aa( i+imax*(j-1)+1, l ) = -smalle(i,j,2,k)
ba( i+imax*(j-1)+1, l ) = smalle(i,j,1,k)
da( i+imax*(j-1)+1, l ) = 0.0d0
ea( i+imax*(j-1)+1, l ) = -cdot(i,j,k) + 0.0d0
end do
end do
!!$ call tridiag( imax*jmax, 2*kmax, aa, ba, da, ea )
call tridiaginv( imax*jmax, 2*kmax, aa, ba, da, ea )
if( IDScheme .eq. IDSchemeShortWave ) then
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
xyr_UwFlux(i,j,k-1) = ea( i+imax*(j-1)+1, 2*k ) * smalle(i,j,1,k) + ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,2,k) + cupb(i,j,k)
xyr_DwFlux(i,j,k-1) = ea( i+imax*(j-1)+1, 2*k ) * smalle(i,j,3,k) + ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,4,k) + cdob(i,j,k)
end do
end do
end do
k = kmax
do j = 1, jmax
do i = 0, imax-1
xyr_UwFlux(i,j,k) = ea( i+imax*(j-1)+1, 2*k ) * smalle(i,j,3,k) - ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,4,k) + cupt(i,j,k)
xyr_DwFlux(i,j,k) = ea( i+imax*(j-1)+1, 2*k ) * smalle(i,j,1,k) - ea( i+imax*(j-1)+1, 2*k-1) * smalle(i,j,2,k) + cdot(i,j,k)
end do
end do
!
! Addition of Direct Solar Insident
!
do k = 0, kmax
xyr_DwFlux(:,:,k) = xyr_DwFlux(:,:,k) + SolarFluxTOA * exp_opdep2(:,:,k) * xy_cosSZA
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
xyr_RadFlux = xyr_UwFlux - xyr_DwFlux
if( IDScheme .eq. IDSchemeShortWave ) then
do j = 1, jmax
do i = 0, imax-1
if( xy_cosSZA(i,j) <= 0.0d0 ) then
do k = 0, kmax
xyr_RadFlux(i,j,k) = 0.0d0
end do
end if
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 RadiationTwoStreamApp