Class radiation_RL78
In: radiation/radiation_RL78.f90

Roewe and Liou (1978) による長波放射モデル

Long radiation model described by Roewe and Liou (1978)

Note that Japanese and English are described in parallel.

長波放射モデル.

This is a model of long wave radiation.

References

 Roewe, D., and K.-N. Liou, Influence of cirrus clouds on the infrared cooling
   rate in the troposphere and lower stratosphere, J. Appl. Met., 17, 92-106, 1978.

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

Methods

Included Modules

dc_types dc_date_types gridset constants timeset dc_date read_time_series set_cloud dc_message dc_iounit namelist_util

Public Instance methods

Subroutine :
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QVap(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QO3(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
xyz_Press(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xy_SurfTemp(0:imax-1, 1:jmax) :real(DP), intent(in )
xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) :real(DP), intent(out)

[Source]

  subroutine RadiationRL78Flux( xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, xy_SurfTemp, xyr_RadLFlux, xyra_DelRadLFlux )


    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav, PI      ! $ \pi $ .
                                 ! Circular constant

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimeN, EndTime, TimesetClockStart, TimesetClockStop


    ! 日付および時刻の取り扱い
    ! Date and time handler
    !
    use dc_date, only: operator(-), operator(>=), operator(+), operator(==), toChar

    use read_time_series, only: SetValuesFromTimeSeriesWrapper

    use set_cloud, only : SetCloudLW

    ! OpenMP
    !
    !$ use omp_lib

    real(DP), intent(in ):: xyz_Temp        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QVap        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QO3         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyr_Press       (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ):: xyz_Press       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xy_SurfTemp     (0:imax-1, 1:jmax)
    real(DP), intent(out):: xyr_RadLFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


    ! 
    ! Work variables
    !
    real(DP):: xyz_QCO2              (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_CloudDelOptDep    (0:imax-1, 1:jmax, 1:kmax)

    logical :: flag_dry_atmosphere

    logical :: a_TransFlag(1:nbmax)


    integer :: n

    integer :: nt

    integer :: js
    integer :: je

    integer :: nthreads
    integer, allocatable :: a_js(:)
    integer, allocatable :: a_je(:)



    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    if ( sw_fs ) then
      call RadiationRoeweLiouInit
    end if


    nthreads = 1
    !$ nthreads  = omp_get_max_threads()

    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



!!$    xyz_QCO2(:,:,:) = 382.0d-6 * CO2MolWeight / MeanMolWeight
    xyz_QCO2(:,:,:) = VMRCO2 * CO2MolWeight / MeanMolWeight
!!$    xyz_QO3 (:,:,:) = 0.0d0


!!$    xyz_QVap(:,:,:) = 0.0d0
!!$    xyz_QCO2(:,:,:) = 300.0d-6 * CO2MolWeight / MeanMolWeight
!!$    xyz_QCO2(:,:,:) = 0.0d0
!!$    call SetRefAtmPro(      &
!!$      & 5, .true., xyz_Press, &
!!$      & xyz_QO3                 &
!!$      & )


!!$    if ( .true. ) then
    if ( ( TimeN - PrevTimeSave >= IntTimeSave ) .or. ( .not. FlagTransSaved ) ) then

!!$      write( 6, * ) 'CalcTrans'

      if ( .not. FlagTransSaved ) then
        PrevTimeSave = TimeN
      else
        PrevTimeSave = PrevTimeSave + IntTimeSave
      end if

      FlagTransSaved = .true.


      ! Check for dry atmosphere
      !
      flag_dry_atmosphere = .false.
      if ( flag_save_time ) then
        if ( all( xyz_QVap <= 0.0d0 ) ) then
          flag_dry_atmosphere = .true.
!!$          write( 6, * ) 'Dry atmosphere'
        else
          flag_dry_atmosphere = .false.
        end if
      end if


      call SetCloudLW( xyz_CloudDelOptDep )




      !$OMP PARALLEL DEFAULT(PRIVATE) &
      !$OMP SHARED(nthreads,a_js,a_je, &
      !$OMP        xyz_QCO2, &
      !$OMP        xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, xy_SurfTemp,  &
      !$OMP        xyz_CloudDelOptDep, &
      !$OMP        flag_dry_atmosphere)
      !$OMP DO

      LOOP_thread_trans : do nt = 0, nthreads-1

        js = a_js(nt)
        je = a_je(nt)

        if ( js > je ) cycle

!!$      write( 6, * ) n, js, je


        call RadiationRL78CalcTotTrans( xyz_QCO2, xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, xyz_CloudDelOptDep, flag_dry_atmosphere, js, je )

      end do LOOP_thread_trans

      !$OMP END DO
      !$OMP END PARALLEL


    end if



    a_TransFlag = .false.
    if ( flag_save_time ) then
      do n = 1, nbmax
        if ( all( xyrra_TransSave(:,:,0,kmax,n) == 1.0d0 ) ) then
          a_TransFlag(n) = .true.
        else
          a_TransFlag(n) = .false.
        end if
      end do
    end if


    !$OMP PARALLEL DEFAULT(PRIVATE) &
    !$OMP SHARED(nthreads,a_js,a_je, &
    !$OMP        a_TransFlag,                       &
    !$OMP        xyz_Temp, xy_SurfTemp,             &
    !$OMP        xyr_RadLFlux, xyra_DelRadLFlux)
    !$OMP DO

    LOOP_thread_RTE : do nt = 0, nthreads-1

      js = a_js(nt)
      je = a_je(nt)

      if ( js > je ) cycle


      call RadiationRL78IntegRTE( a_TransFlag, xyz_Temp, xy_SurfTemp, xyr_RadLFlux, xyra_DelRadLFlux, js, je )

    end do LOOP_thread_RTE

    !$OMP END DO
    !$OMP END PARALLEL


    deallocate( a_js )
    deallocate( a_je )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )


  end subroutine RadiationRL78Flux

Private Instance methods

AvogNum
Constant :
AvogNum = 6.0221415d23 :real(DP), parameter
CO2MolWeight
Variable :
CO2MolWeight :real(DP), save
CO2ScaleIndex
Variable :
CO2ScaleIndex :real(DP), save
Subroutine :
wn1 :real(DP), intent(in )
wn2 :real(DP), intent(in )
num :integer , intent(in )
im :integer , intent(in )
jm :integer , intent(in )
km :integer , intent(in )
temp( im, jm, km ) :real(DP), intent(in )
pfinted( im, jm, km ) :real(DP), intent(out)

[Source]

  subroutine CalcIntegratedPF( wn1, wn2, num, im, jm, km, temp, pfinted )


    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: im, jm, km
    real(DP), intent(in ) :: temp   ( im, jm, km )
    real(DP), intent(out) :: pfinted( im, jm, km )


    !
    ! local variables
    !
    real(DP):: x( num ), w( num )
    integer :: i
    integer :: j
    integer :: k
    integer :: l


    call gauleg( wn1, wn2, num, x, w )

    pfinted(:,:,:) = 0.0d0

    do l = 1, num
      do k = 1, km
        do j = 1, jm
          do i = 1, im
            pfinted(i,j,k) = pfinted(i,j,k) + pf( x( l ), temp(i,j,k) ) * w( l )
          end do
        end do
      end do
    end do


  end subroutine CalcIntegratedPF
Subroutine :
wn1 :real(DP), intent(in )
wn2 :real(DP), intent(in )
num :integer , intent(in )
im :integer , intent(in )
jm :integer , intent(in )
temp( im, jm ) :real(DP), intent(in )
pfinted( im, jm ) :real(DP), intent(out)

[Source]

  subroutine CalcIntegratedPF2D( wn1, wn2, num, im, jm, temp, pfinted )


    real(DP), intent(in ) :: wn1,wn2
    integer , intent(in ) :: num
    integer , intent(in ) :: im, jm
    real(DP), intent(in ) :: temp   ( im, jm )
    real(DP), intent(out) :: pfinted( im, jm )


    !
    ! local variables
    !
    real(dp)     :: temp3d( im, jm, 1 ), pfinted3d( im, jm, 1 )


    temp3d(:,:,1) = temp(:,:)
    call CalcIntegratedPF( wn1, wn2, num, im, jm, 1, temp3d, pfinted3d )
    pfinted(:,:) = pfinted3d(:,:,1)


  end subroutine CalcIntegratedPF2D
Subroutine :
iband :integer , intent(in )
xy_IntegPF(0:imax-1, 1:jmax) :real(DP), intent(out)
js :integer , intent(in )
je :integer , intent(in )
flag_DPFDT :logical , intent(in ), optional

[Source]

  subroutine CalcIntegratedPFWithTable2D( iband, xy_Temp, xy_IntegPF, js, je, flag_DPFDT )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax    ! 
                               ! Number of grid points in latitude

    integer , intent(in )           :: iband
    real(DP), intent(in )           :: xy_temp   (0:imax-1, 1:jmax)
    real(DP), intent(out)           :: xy_IntegPF(0:imax-1, 1:jmax)
    integer , intent(in )           :: js
    integer , intent(in )           :: je
    logical , intent(in ), optional :: flag_DPFDT

    !
    ! local variables
    !
    real(DP) :: xyz_Temp   (0:imax-1, 1:jmax, 1)
    real(DP) :: xyz_IntegPF(0:imax-1, 1:jmax, 1)
    integer  :: j


    do j = js, je
      xyz_Temp(:,j,1) = xy_Temp(:,j)
    end do

    call CalcIntegratedPFWithTable3D( iband, 1, xyz_temp, xyz_IntegPF, js, je, flag_DPFDT )

    do j = js, je
      xy_IntegPF(:,j) = xyz_IntegPF(:,j,1)
    end do


  end subroutine CalcIntegratedPFWithTable2D
Subroutine :
iband :integer , intent(in )
km :integer , intent(in )
xyz_temp(0:imax-1, 1:jmax, 1:km) :real(DP), intent(in )
xyz_IntegPF(0:imax-1, 1:jmax, 1:km) :real(DP), intent(out)
js :integer , intent(in )
je :integer , intent(in )
flag_DPFDT :logical , intent(in ), optional

[Source]

  subroutine CalcIntegratedPFWithTable3D( iband, km, xyz_temp, xyz_IntegPF, js, je, flag_DPFDT )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    integer , intent(in ) :: iband
    integer , intent(in ) :: km
    real(DP), intent(in ) :: xyz_temp   (0:imax-1, 1:jmax, 1:km)
    real(DP), intent(out) :: xyz_IntegPF(0:imax-1, 1:jmax, 1:km)
    logical , intent(in ), optional :: flag_DPFDT
    integer , intent(in )           :: js
    integer , intent(in )           :: je

    !
    ! local variables
    !
    logical                     :: local_flag_DPFDT

    integer                     :: xyz_TempIndex(0:imax-1, 1:jmax, 1:km)
    integer                     :: i
    integer                     :: j
    integer                     :: k
    integer                     :: m
    integer                     :: n


    do k = 1, km
      do j = js, je
        do i = 0, imax-1

          if ( ( xyz_Temp(i,j,k) < a_TableTemp(1)     ) .or. ( xyz_Temp(i,j,k) > a_TableTemp(ntmax) ) ) then
            call MessageNotify( 'E', module_name, 'Temperature is not appropriate, Temp(%d,%d,%d) = %f.', i = (/i, j, k/), d = (/xyz_Temp(i,j,k)/) )
          end if

          xyz_TempIndex(i,j,k) = int( ( xyz_Temp(i,j,k) - TableTempMin ) / TableTempIncrement ) + 2

          if ( xyz_TempIndex(i,j,k) == 1 ) then
             xyz_TempIndex(i,j,k) = 2
          else if ( xyz_TempIndex(i,j,k) >= ntmax ) then
             xyz_TempIndex(i,j,k) = ntmax - 1
          end if

!!$          xyz_TempIndex(i,j,k) = ntmax-1
!!$          search_index: do m = 2, ntmax-1
!!$            if ( a_TableTemp(m) >= xyz_Temp(i,j,k) ) then
!!$              xyz_TempIndex(i,j,k) = m
!!$              exit search_index
!!$            end if
!!$          end do search_index

        end do
      end do
    end do


    local_flag_DPFDT = .false.
    if ( present( flag_DPFDT ) ) then
      if ( flag_DPFDT ) then
        local_flag_DPFDT = .true.
      end if
    end if

    if ( .not. local_flag_DPFDT ) then
      do k = 1, km
        do j = js, je
          do i = 0, imax-1
            m = xyz_TempIndex(i,j,k)

!!$            xyz_IntegPF(i,j,k) = &
!!$              &   ( aa_TableIPF( m, iband ) - aa_TableIPF( m-1, iband ) ) &
!!$              & / ( a_TableTemp( m )        - a_TableTemp( m-1 )        ) &
!!$              & * ( xyz_Temp(i,j,k)         - a_TableTemp( m-1 )        ) &
!!$              & +   aa_TableIPF( m-1, iband )

            xyz_IntegPF(i,j,k) = aa_TableIPF(m-1,iband) * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) ) * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m-1 ) - a_TableTemp( m   ) ) * ( a_TableTemp( m-1 ) - a_TableTemp( m+1 ) ) ) + aa_TableIPF(m  ,iband) * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m   ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m   ) - a_TableTemp( m+1 ) ) ) + aa_TableIPF(m+1,iband) * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) ) / ( ( a_TableTemp( m+1 ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m+1 ) - a_TableTemp( m   ) ) )
          end do
        end do
      end do
    else
      do k = 1, km
        do j = js, je
          do i = 0, imax-1
            m = xyz_TempIndex(i,j,k)

!!$            xyz_IntegPF(i,j,k) = &
!!$              &   ( aa_TableIDPFDT( m, iband ) - aa_TableIDPFDT( m-1, iband ) ) &
!!$              & / ( a_TableTemp   ( m )        - a_TableTemp   ( m-1 )        ) &
!!$              & * ( xyz_Temp(i,j,k)         - a_TableTemp( m-1 )        ) &
!!$              & +   aa_TableIDPFDT( m-1, iband )

            xyz_IntegPF(i,j,k) = aa_TableIDPFDT(m-1,iband) * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) ) * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m-1 ) - a_TableTemp( m   ) ) * ( a_TableTemp( m-1 ) - a_TableTemp( m+1 ) ) ) + aa_TableIDPFDT(m  ,iband) * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp   (i,j,k) - a_TableTemp( m+1 ) ) / ( ( a_TableTemp( m   ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m   ) - a_TableTemp( m+1 ) ) ) + aa_TableIDPFDT(m+1,iband) * ( xyz_Temp   (i,j,k) - a_TableTemp( m-1 ) ) * ( xyz_Temp   (i,j,k) - a_TableTemp( m   ) ) / ( ( a_TableTemp( m+1 ) - a_TableTemp( m-1 ) ) * ( a_TableTemp( m+1 ) - a_TableTemp( m   ) ) )
          end do
        end do
      end do
    end if


  end subroutine CalcIntegratedPFWithTable3D
Function :
DPFDT :real(DP)
WN :real(DP), intent(in )
Temp :real(DP), intent(in )

[Source]

  real(DP) function DPFDT( WN, Temp )

    ! USE statements
    !

    real(DP), intent(in ) :: WN
    real(DP), intent(in ) :: Temp


    real(DP), parameter :: c      = 2.99792458d8           , planc  = 6.6260755d-34          , boltz  = 1.380658d-23
    real(DP) :: ExpTerm
    real(DP) :: PF


    ExpTerm = exp( planc * c * ( WN + 1.0d-10 ) / ( boltz * Temp ) )

    PF = 2.0d0 * planc * c * c * WN * WN * WN / ( ExpTerm - 1.0d0 )

    DPFDT = 1.0d0 / ( 2.0d0 * c * WN * WN * boltz ) * ( pf / Temp )**2 * ExpTerm


  end function DPFDT
Subroutine :
WN :real(DP), intent(in )
xy_Temp(0:imax-1, 1:jmax) :real(DP), intent(in )
xy_DPFDT(0:imax-1, 1:jmax) :real(DP), intent(out)

[Source]

  subroutine DPFDT_Array2D( WN, xy_Temp, xy_DPFDT )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    real(DP), intent(in ) :: WN
    real(DP), intent(in ) :: xy_Temp (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xy_DPFDT(0:imax-1, 1:jmax)


    real(DP), parameter :: c      = 2.99792458d8           , planc  = 6.6260755d-34          , boltz  = 1.380658d-23
    real(DP) :: xy_ExpTerm(0:imax-1, 1:jmax)
    real(DP) :: xy_PF     (0:imax-1, 1:jmax)


    xy_ExpTerm(:,:) = exp( planc * c * ( WN + 1.0d-10 ) / ( boltz * xy_Temp ) )

    xy_PF(:,:) = 2.0d0 * planc * c * c * WN * WN * WN / ( xy_ExpTerm(:,:) - 1.0d0 )

    xy_DPFDT(:,:) = 1.0d0 / ( 2.0d0 * c * WN * WN * boltz ) * ( xy_pf(:,:) / xy_Temp(:,:) )**2 * xy_ExpTerm(:,:)


  end subroutine DPFDT_ARRAY2D
DiffFactor
Constant :
DiffFactor = 1.66d0 :real(DP), parameter
FlagTransSaved
Variable :
FlagTransSaved :logical , save
H2OMolWeight
Variable :
H2OMolWeight :real(DP), save
H2OScaleIndex
Variable :
H2OScaleIndex :real(DP), save
IntTimeSave
Variable :
IntTimeSave :type(DC_DIFFTIME), save
: 長波フラックスを計算する時間間隔. Interval time of long wave flux calculation
Subroutine :
WN1 :real(DP), intent(in )
WN2 :real(DP), intent(in )
Num :integer , intent(in )
xy_Temp(0:imax-1, 1:jmax) :real(DP), intent(in )
xy_DPFDTInted(0:imax-1, 1:jmax) :real(DP), intent(out)

[Source]

  subroutine Integ_DPFDT_GQ_Array2D( WN1, WN2, Num, xy_Temp, xy_DPFDTInted )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    real(DP), intent(in ) :: WN1
    real(DP), intent(in ) :: WN2
    integer , intent(in ) :: Num
    real(DP), intent(in ) :: xy_Temp      (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xy_DPFDTInted(0:imax-1, 1:jmax)


    !
    ! local variables
    !
    real(DP):: GP( Num )
    real(DP):: GW( Num )
    real(DP):: xy_DPFDT(0:imax-1, 1:jmax)
    integer :: l


    call gauleg( WN1, WN2, Num, GP, GW )

    xy_DPFDTInted(:,:) = 0.0d0

    do l = 1, num
      call DPFDT_Array2D( GP(l), xy_Temp, xy_DPFDT )

      xy_DPFDTInted(:,:) = xy_DPFDTInted(:,:) + xy_DPFDT(:,:) * GW(l)
    end do


  end subroutine Integ_DPFDT_GQ_Array2D
MeanMolWeight
Variable :
MeanMolWeight :real(DP), save
NGaussQuad
Variable :
NGaussQuad :integer , save
O3ScaleIndex
Variable :
O3ScaleIndex :real(DP), save
PrevTimeSave
Variable :
PrevTimeSave :type(DC_DIFFTIME), save
: 前回長波フラックスを計算した時刻. Time when long wave flux is calculated
Subroutine :
xyz_QCO2(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QVap(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QO3(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyr_Press(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
xyz_Press(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_CloudDelOptDep(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
flag_dry_atmosphere :logical , intent(in )
js :integer , intent(in )
je :integer , intent(in )

[Source]

  subroutine RadiationRL78CalcTotTrans( xyz_QCO2, xyz_Temp, xyz_QVap, xyz_QO3, xyr_Press, xyz_Press, xyz_CloudDelOptDep, flag_dry_atmosphere, js, je )


    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav    ! $ g $ [m s-2].
                                 ! 
                                 ! Gravitational acceleration

    real(DP), intent(in ):: xyz_QCO2          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_Temp          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QVap          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QO3           (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyr_Press         (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ):: xyz_Press         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudDelOptDep(0:imax-1, 1:jmax, 1:kmax)
    logical , intent(in ):: flag_dry_atmosphere
    integer , intent(in ):: js
    integer , intent(in ):: je


    ! 
    ! Work variables
    !
    real(DP):: RefPress

    real(DP):: xyz_H2ODelAbsAmt      (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_H2ODelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_O3DelScaledAbsAmt (0:imax-1, 1:jmax, 1:kmax)

    real(DP):: xyrr_TransCO2         (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    real(DP):: xyz_TransCloudOneLayer(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyrr_TransCloud       (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    real(DP):: xyrr_Trans            (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    integer :: j, k, k1, k2
    integer :: n



    if ( sw_fs ) then
      call RadiationRoeweLiouInit
    end if


    do k = 1, kmax
      do j = js, je
        xyz_H2ODelAbsAmt(:,j,k) = xyz_QVap(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do

    RefPress = 1.0d5

    do k = 1, kmax
      do j = js, je
        xyz_H2ODelScaledAbsAmt(:,j,k) = ( xyz_Press(:,j,k) / RefPress )**H2OScaleIndex * xyz_QVap(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do
    do k = 1, kmax
      do j = js, je
        xyz_CO2DelScaledAbsAmt(:,j,k) = ( xyz_Press(:,j,k) / RefPress )**CO2ScaleIndex * xyz_QCO2(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do
    do k = 1, kmax
      do j = js, je
        xyz_O3DelScaledAbsAmt(:,j,k) = ( xyz_Press(:,j,k) / RefPress )**O3ScaleIndex * xyz_QO3(:,j,k) * ( xyr_Press(:,j,k-1) - xyr_Press(:,j,k) ) / Grav
      end do
    end do



    !-----

    ! Transmission in CO2 15 micron band
    !
    call RadiationRoeweLiouCalcCO2Trans( xyz_CO2DelScaledAbsAmt, xyrr_TransCO2, js, je )


    !-----

    ! Transmission by cloud

!!$      if ( all( xyz_CloudDelOptDep <= 0.0d0 ) ) then
!!$
!!$        xyrr_TransCloud(:,:,:,:) = 1.0d0
!!$
!!$      else

    do k = 1, kmax
      do j = js, je
        xyz_TransCloudOneLayer(:,j,k) = exp( - xyz_CloudDelOptDep(:,j,k) * DiffFactor )
      end do
    end do

    do k1 = 0, kmax
      k2 = k1
      do j = js, je
        xyrr_TransCloud(:,j,k1,k2) = 1.0d0
      end do
      do k2 = k1+1, kmax
        do j = js, je
          xyrr_TransCloud(:,j,k1,k2) = xyrr_TransCloud(:,j,k1,k2-1) * xyz_TransCloudOneLayer(:,j,k2)
        end do
      end do
    end do
    do k1 = 0, kmax
      do k2 = 0, k1-1
        do j = js, je
          xyrr_TransCloud(:,j,k1,k2) = xyrr_TransCloud(:,j,k2,k1)
        end do
      end do
    end do

!!$      end if

    !-----

    do n = 1, nbmax

      if ( n <= 5 ) then
        ! H2O band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans, js, je )
        end if
      else if ( n <= 8 ) then
        ! H2O band & CO2 band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans, js, je )
        end if
        do k2 = 0, kmax
          do k1 = 0, kmax
            do j = js, je
              xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k1,k2) * xyrr_TransCO2(:,j,k1,k2)
            end do
          end do
        end do
      else if ( n == 9 ) then
        ! H2O band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans, js, je )
        end if
      else if ( n == 10 ) then
        ! H2O band & continuum
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans, js, je )
          call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
        end if
      else if ( n == 11 ) then
        ! H2O continuum
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
          call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
        end if
      else if ( n <= 21 ) then
        ! O3 band & H2O continuum
        call RadiationRoeweLiouCalcBandTrans( n, xyz_O3DelScaledAbsAmt, xyrr_Trans, js, je )
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k1,k2)
              end do
            end do
          end do
        else
          call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
        end if
      else if ( n == 22 ) then
        ! H2O continuum
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
          call RadiationRLCalcH2OContTrans( n, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )
        end if
      else if ( n <= 31 ) then
        ! H2O band
        if ( flag_dry_atmosphere ) then
          do k2 = 0, kmax
            do k1 = 0, kmax
              do j = js, je
                xyrr_Trans(:,j,k1,k2) = 1.0d0
              end do
            end do
          end do
        else
          call RadiationRoeweLiouCalcBandTrans( n, xyz_H2ODelScaledAbsAmt, xyrr_Trans, js, je )
        end if
      else if ( n == 32 ) then
        do k2 = 0, kmax
          do k1 = 0, kmax
            do j = js, je
              xyrr_Trans(:,j,k1,k2) = 1.0d0
            end do
          end do
        end do
      else
        write( 6, * ) 'Unexpected number of band, ', n, '.'
        stop
      end if

      do k2 = 0, kmax
        do k1 = 0, kmax
          do j = js, je
            xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k1,k2) * xyrr_TransCloud(:,j,k1,k2)

            xyrra_TransSave(:,j,k1,k2,n) = xyrr_Trans(:,j,k1,k2)
          end do
        end do
      end do

    end do


  end subroutine RadiationRL78CalcTotTrans
Subroutine :
a_TransFlag(1:nbmax) :logical , intent(in )
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xy_SurfTemp(0:imax-1, 1:jmax) :real(DP), intent(in )
xyr_RadLFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) :real(DP), intent(out)
js :integer , intent(in )
je :integer , intent(in )

[Source]

  subroutine RadiationRL78IntegRTE( a_TransFlag, xyz_Temp, xy_SurfTemp, xyr_RadLFlux, xyra_DelRadLFlux, js, je )


    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    ! 
    ! Physical constants settings
    !
    use constants, only: PI      ! $ \pi $ .
                                 ! Circular constant

    logical , intent(in ):: a_TransFlag     (1:nbmax)
    real(DP), intent(in ):: xyz_Temp        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xy_SurfTemp     (0:imax-1, 1:jmax)
    real(DP), intent(out):: xyr_RadLFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    integer , intent(in ):: js
    integer , intent(in ):: je



    ! 
    ! Work variables
    !
    real(DP):: xyrr_Trans            (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    real(DP):: xyz_PlankFunc         (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xy_SurfPlankFunc      (0:imax-1, 1:jmax)
    real(DP):: xy_DPFDT0             (0:imax-1, 1:jmax)
    real(DP):: xy_DPFDT1             (0:imax-1, 1:jmax)


    integer :: j
    integer :: k, kk, k1, k2
    integer :: l
    integer :: n


    do k = 0, kmax
      do j = js, je
        xyr_RadLFlux(:,j,k)   = 0.0d0
      end do
    end do
    do l = 0, 1
      do k = 0, kmax
        do j = js, je
          xyra_DelRadLFlux(:,j,k,l) = 0.0d0
        end do
      end do
    end do


    LOOP_band: do n = 1, nbmax

      do k2 = 0, kmax
        do k1 = 0, kmax
          do j = js, je
            xyrr_Trans(:,j,k1,k2) = xyrra_TransSave(:,j,k1,k2,n)
          end do
        end do
      end do


      if ( a_TransFlag(n) ) then

        call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_SurfPlankFunc, js, je )

        do k = 0, kmax
          do j = js, je
            xyr_RadLFlux(:,j,k) = xyr_RadLFlux(:,j,k) + PI * xy_SurfPlankFunc(:,j) ! * xyrr_Trans(:,:,0,k)
          end do
        end do

        call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_DPFDT0, js, je, .true. )

        do k = 0, kmax
          do j = js, je
            xyra_DelRadLFlux(:,j,k,0) = xyra_DelRadLFlux(:,j,k,0) + PI * xy_DPFDT0(:,j) ! * xyrr_Trans(:,:,0,k)
          end do
        end do

      else

        call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_SurfPlankFunc, js, je )
        call CalcIntegratedPFWithTable3D( n, kmax, xyz_Temp, xyz_PlankFunc, js, je )

        do k = 0, kmax
          do j = js, je
            xyr_RadLFlux(:,j,k) = xyr_RadLFlux(:,j,k) + PI * xy_SurfPlankFunc(:,j) * xyrr_Trans(:,j,0,k)
          end do
          do kk = 1, kmax
            do j = js, je
              xyr_RadLFlux(:,j,k) = xyr_RadLFlux(:,j,k) - PI * xyz_PlankFunc(:,j,kk) * ( xyrr_Trans(:,j,k,kk-1) - xyrr_Trans(:,j,k,kk) )
            end do
          end do
        end do

        call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_DPFDT0, js, je, .true. )
        call CalcIntegratedPFWithTable2D( n, xyz_Temp(:,:,1), xy_DPFDT1, js, je, .true. )

        do k = 0, kmax
          do j = js, je
            xyra_DelRadLFlux(:,j,k,0) = xyra_DelRadLFlux(:,j,k,0) + PI * xy_DPFDT0(:,j) * xyrr_Trans(:,j,0,k)
            xyra_DelRadLFlux(:,j,k,1) = xyra_DelRadLFlux(:,j,k,1) - PI * xy_DPFDT1(:,j) * ( xyrr_Trans(:,j,k,0) - xyrr_Trans(:,j,k,1) )
          end do
        end do

      end if


    end do LOOP_band


  end subroutine RadiationRL78IntegRTE
Subroutine :
iband :integer , intent(in )
xyz_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QVap(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_Press(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) :real(DP), intent(inout)

[Source]

  subroutine RadiationRLCalcH2OContTrans( iband, xyz_Temp, xyz_QVap, xyz_Press, xyz_H2ODelAbsAmt, xyrr_Trans )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    integer , intent(in   ):: iband
    real(DP), intent(in   ):: xyz_Temp        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ):: xyz_QVap        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ):: xyz_Press       (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ):: xyz_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout):: xyrr_Trans      (0:imax-1, 1:jmax, 0:kmax, 0:kmax)


    ! 
    ! Work variables
    !
    real(DP):: H2OContConstant
    real(DP):: a_H2OContParam(3)
    real(DP):: RefTemp
    real(DP):: gamma
    real(DP):: xyz_PressH2O      (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DelOptDep     (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyrr_H2OContOptDep(0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    integer :: k1
    integer :: k2
    integer :: k3
    integer :: n


    if ( ( iband < 9 ) .and. ( iband > 11 ) ) then
      write( 6, * ) 'Unexpected number of band in RadiationRLCalcH2OContTrans.'
      stop
    end if


    ! These constants are obtained from Roberts et al. (1976) [Applied Optics]
    a_H2OContParam(1) = 1.25d-22 ! mol-1 cm2 atm-1
    a_H2OContParam(2) = 1.67d-19 ! mol-1 cm2 atm-1
    a_H2OContParam(3) = 7.87d-3  ! (cm)

    ! unit conversion
    a_H2OContParam(1) = a_H2OContParam(1) / ( H2OMolWeight / AvogNum ) * 1.0d-4 / 101325.0d0 ! (kg-1 m2 Pa-1)
    a_H2OContParam(2) = a_H2OContParam(2) / ( H2OMolWeight / AvogNum ) * 1.0d-4 / 101325.0d0 ! (kg-1 m2 Pa-1)
    a_H2OContParam(3) = a_H2OContParam(3) * 1.0d-2    ! (m)


    n = iband


    ! This constant is deribed from Roberts et al. (1976) [Applied Optics]
!!$    gamma = 0.005d0
    gamma = 0.0d0

    xyz_PressH2O = xyz_Press * xyz_QVap * MeanMolWeight / H2OMolWeight

    RefTemp = 1800.0d0

    H2OContConstant = a_H2OContParam(1) + a_H2OContParam(2) * exp( - a_H2OContParam(3) * ( a_BandParam(1,n) + a_BandParam(2,n) ) * 0.5d0 )

    xyz_DelOptDep = H2OContConstant * exp( RefTemp * ( 1.0d0 / xyz_Temp - 1.0d0 / 296.0d0 ) ) * ( xyz_PressH2O + gamma * ( xyz_Press - xyz_PressH2O ) ) * xyz_H2ODelAbsAmt

    xyrr_H2OContOptDep(:,:,:,:) = 0.0d0
    do k1 = 0, kmax
      do k2 = k1, kmax

        do k3 = k1+1, k2
          xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) + xyz_DelOptDep(:,:,k3)
        end do

        xyrr_H2OContOptDep(:,:,k1,k2) = xyrr_H2OContOptDep(:,:,k1,k2) * DiffFactor

        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k1,k2) * exp( - xyrr_H2OContOptDep(:,:,k1,k2) )

      end do
    end do

    do k1 = 0, kmax
      do k2 = 0, k1-1
        xyrr_Trans(:,:,k1,k2) = xyrr_Trans(:,:,k2,k1)
      end do
    end do
    do k1 = 0, kmax
      k2 = k1
      xyrr_Trans(:,:,k1,k2) = 1.0d0
    end do


  end subroutine RadiationRLCalcH2OContTrans
Subroutine :
iband :integer , intent(in )
xyz_DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) :real(DP), intent(out)
js :integer , intent(in )
je :integer , intent(in )

[Source]

  subroutine RadiationRoeweLiouCalcBandTrans( iband, xyz_DelScaledAbsAmt, xyrr_Trans, js, je )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    integer , intent(in ):: iband
    real(DP), intent(in ):: xyz_DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyrr_Trans         (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    integer , intent(in ):: js
    integer , intent(in ):: je

    ! 
    ! Work variables
    !
    real(DP):: xyrr_AbsAmt(0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    integer :: j
    integer :: k1
    integer :: k2
    integer :: k3
    integer :: n


    n = iband


    do k2 = 0, kmax
      do k1 = 0, kmax
        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = 0.0d0
        end do
      end do
    end do

    do k1 = 0, kmax
      do k2 = k1, kmax

        do k3 = k1+1, k2
          do j = js, je
            xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) + xyz_DelScaledAbsAmt(:,j,k3)
          end do
        end do

        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) * DiffFactor
        end do

        do j = js, je
          xyrr_Trans(:,j,k1,k2) = - a_BandParam(3,n) * xyrr_AbsAmt(:,j,k1,k2) / sqrt( 1.0d0 + a_BandParam(3,n) / a_BandParam(4,n) * xyrr_AbsAmt(:,j,k1,k2) )
          xyrr_Trans(:,j,k1,k2) = exp( xyrr_Trans(:,j,k1,k2) )
        end do

      end do
    end do


    do k1 = 0, kmax
      do k2 = 0, k1-1
        do j = js, je
          xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k2,k1)
        end do
      end do
    end do
    do k1 = 0, kmax
      k2 = k1
      do j = js, je
        xyrr_Trans(:,j,k1,k2) = 1.0d0
      end do
    end do


  end subroutine RadiationRoeweLiouCalcBandTrans
Subroutine :
xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyrr_Trans(0:imax-1, 1:jmax, 0:kmax, 0:kmax) :real(DP), intent(out)
js :integer , intent(in )
je :integer , intent(in )

[Source]

  subroutine RadiationRoeweLiouCalcCO2Trans( xyz_CO2DelScaledAbsAmt, xyrr_Trans, js, je )

    ! USE statements
    !

    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

!!$    integer , intent(in ):: iband
    real(DP), intent(in ):: xyz_CO2DelScaledAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyrr_Trans            (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
    integer , intent(in ):: js
    integer , intent(in ):: je


    ! 
    ! Work variables
    !
    real(DP):: xyrr_AbsAmt       (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    integer :: j
    integer :: k1
    integer :: k2
    integer :: k3
!!$    integer :: n


!!$    n = iband
!!$    n = 2


    do k2 = 0, kmax
      do k1 = 0, kmax
        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = 0.0d0
        end do
      end do
    end do

    do k1 = 0, kmax
      do k2 = k1, kmax

        do k3 = k1+1, k2
          do j = js, je
            xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) + xyz_CO2DelScaledAbsAmt(:,j,k3)
          end do
        end do

        do j = js, je
          xyrr_AbsAmt(:,j,k1,k2) = xyrr_AbsAmt(:,j,k1,k2) * DiffFactor
        end do

        do j = js, je
          xyrr_Trans(:,j,k1,k2) = - a_CO2BandParam(3) * xyrr_AbsAmt(:,j,k1,k2) / sqrt( 1.0d0 + a_CO2BandParam(3) / a_CO2BandParam(4) * xyrr_AbsAmt(:,j,k1,k2) )
          xyrr_Trans(:,j,k1,k2) = exp( xyrr_Trans(:,j,k1,k2) )
        end do

      end do
    end do


    do k1 = 0, kmax
      do k2 = 0, k1-1
        do j = js, je
          xyrr_Trans(:,j,k1,k2) = xyrr_Trans(:,j,k2,k1)
        end do
      end do
    end do
    do k1 = 0, kmax
      k2 = k1
      do j = js, je
        xyrr_Trans(:,j,k1,k2) = 1.0d0
      end do
    end do


  end subroutine RadiationRoeweLiouCalcCO2Trans
Subroutine :

This procedure input/output NAMELIST#radiation_RL78_nml .

[Source]

  subroutine RadiationRoeweLiouInit


    ! 
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 
                               ! Number of vertical level

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 日付および時刻の取り扱い
    ! Date and time handler
    !
    use dc_date, only: DCDiffTimeCreate

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    real(DP)          :: DelTimeCalcTransValue
    character(STRING) :: DelTimeCalcTransUnit

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    logical            :: FlagCheckLoopExit
    real(DP)           :: xy_TempTMP   (0:imax-1, 1:jmax)
    real(DP)           :: xy_PF        (0:imax-1, 1:jmax)
    real(DP)           :: xy_DPFDT     (0:imax-1, 1:jmax)
    real(DP)           :: xy_PFTable   (0:imax-1, 1:jmax)
    real(DP)           :: xy_DPFDTTable(0:imax-1, 1:jmax)
    real(DP)           :: ErrorPFInteg
    real(DP), parameter:: ThresholdErrorPFInteg = 1.0d-3
                              ! Threshold for checking accuracy of calculation of 
                              ! integrated Planc function by using a pre-calculated 
                              ! table.

    ! Variables for preparation for calculation of Plank function
    !
    real(DP)      , allocatable :: GQP(:)
    real(DP)      , allocatable :: GQW(:)


    integer:: i
    integer:: j
    integer:: l
    integer:: m
    integer:: n


    namelist /radiation_RL78_nml/ VMRCO2, DelTimeCalcTransValue, DelTimeCalcTransUnit, flag_save_time

    VMRCO2                = 382.0d-6

    DelTimeCalcTransValue = 3.0
    DelTimeCalcTransUnit  = 'hrs'
    flag_save_time        = .false.


    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = radiation_RL78_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if

    ! Handle interval time
    !
    call DCDiffTimeCreate( IntTimeSave, DelTimeCalcTransValue, DelTimeCalcTransUnit )   ! (in)


!!$    do n = 1, nbmax
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
!!$      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
!!$    end do



    do n = 1, 12-1
      ! unit conversion from (cm-1) to (m-1)
      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
      ! unit conversion from (g-1 cm2) to (kg-1 m2)
      a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
    end do
    ! O3 unit conversion
    do n = 12, 21
      ! unit conversion from (cm-1) to (m-1)
      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
      ! unit conversion from "wrong (g-1 cm2)" to (kg-1 m2)
      !   The values written by Roewe and Liou (1978) seems to be wrong.
      !   In order to convert cm-1 atm-1 to g-1 cm2, the value has to be MULTIPLIED
      !   by ~400.
      !   The factor ~400 is
      !   1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0 * 1.0d4 * 1.0d-3.
      !   However, it seems that the value was DIVIDED by ~400 to convert the value
      !   in cm-1 atm-1 to g-1 cm2 by Roewe and Liou (1978).
      !   So, first, the value written in Roewe and Liou (1978) is multiplied by ~400
      !   to obtain the correct value in cm-1 atm-1.
      !   Then the value is multiplied by ~40 to convert unit from cm-1 atm-1 to m2 kg-1.
      !   The factor ~40 is
      !   1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0. 
      a_BandParam(3,n) = a_BandParam(3,n) * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0 * 1.0d4 * 1.0d-3 * 1.0d2 / 101325.0d0 * 8.31432d0 / ( 48.0d-3 ) * 233.0d0
    end do
    do n = 21+1, nbmax
      ! unit conversion from (cm-1) to (m-1)
      a_BandParam(1,n) = a_BandParam(1,n) * 1.0d2
      a_BandParam(2,n) = a_BandParam(2,n) * 1.0d2
      ! unit conversion from (g-1 cm2) to (kg-1 m2)
      a_BandParam(3,n) = a_BandParam(3,n) * 1.0d-4 * 1.0d3
    end do


!!$    do n = 1, nbH2Omax
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_H2OBandParam(1,n) = a_H2OBandParam(1,n) * 1.0d2
!!$      a_H2OBandParam(2,n) = a_H2OBandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_H2OBandParam(3,n) = a_H2OBandParam(3,n) * 1.0d-4 * 1.0d3
!!$
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_H2OContParam(1,n) = a_H2OContParam(1,n) * 1.0d2
!!$      a_H2OContParam(2,n) = a_H2OContParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2 atm-1) to (kg-1 m2 Pa-1)
!!$      a_H2OContParam(3,n) = a_H2OContParam(3,n) * 1.0d-4 * 1.0d3 * 1.0d-5
!!$    end do

    ! unit conversion from (cm-1) to (m-1)
    a_CO2BandParam(1) = a_CO2BandParam(1) * 1.0d2
    a_CO2BandParam(2) = a_CO2BandParam(2) * 1.0d2
    ! unit conversion from (g-1 cm2) to (kg-1 m2)
    a_CO2BandParam(3) = a_CO2BandParam(3) * 1.0d-4 * 1.0d3

!!$    do n = 1, nbO3max
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_O3BandParam(1,n) = a_O3BandParam(1,n) * 1.0d2
!!$      a_O3BandParam(2,n) = a_O3BandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_O3BandParam(3,n) = a_O3BandParam(3,n) * 1.0d-4 * 1.0d3
!!$    end do

!!$    do n = 1, 3
!!$      ! unit conversion from (cm-1) to (m-1)
!!$      a_TMPFORO3BandParam(1,n) = a_TMPFORO3BandParam(1,n) * 1.0d2
!!$      a_TMPFORO3BandParam(2,n) = a_TMPFORO3BandParam(2,n) * 1.0d2
!!$      ! unit conversion from (g-1 cm2) to (kg-1 m2)
!!$      a_TMPFORO3BandParam(3,n) = a_TMPFORO3BandParam(3,n) * 1.0d-4 * 1.0d3
!!$    end do


    ! Values by Miyoshi and Morita (19??)
    !
    H2OScaleIndex = 1.0d0
    CO2ScaleIndex = 0.9d0
    O3ScaleIndex  = 0.4d0


    NGaussQuad = 5


    MeanMolWeight = 28.0d-3
    H2OMolWeight  = 18.0d-3
    CO2MolWeight  = 44.0d-3


    allocate( xyrra_TransSave  (0:imax-1, 1:jmax, 0:kmax, 0:kmax, 1:nbmax) )


    ! Preparation of tables for calculation of Plank function
    !
    TableTempMin       =  50.0d0
    TableTempMax       = 400.0d0
    TableTempIncrement =   0.1d0
    ntmax              = ( TableTempMax - TableTempMin ) / TableTempIncrement + 1
    allocate( a_TableTemp   (1:ntmax) )
    allocate( aa_TableIPF   (1:ntmax, 1:nbmax) )
    allocate( aa_TableIDPFDT(1:ntmax, 1:nbmax) )

    do m = 1, ntmax
      a_TableTemp(m) = TableTempMin + TableTempIncrement * ( m - 1 )
    end do


    aa_TableIPF   (:,:) = 0.0d0
    aa_TableIDPFDT(:,:) = 0.0d0

    allocate( GQP(1:NGaussQuad) )
    allocate( GQW(1:NGaussQuad) )
    do n = 1, nbmax
      call gauleg( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, GQP, GQW )
      do m = 1, ntmax
        do l = 1, NGaussQuad
          aa_TableIPF   (m,n) = aa_TableIPF   (m,n) + pf   ( GQP(l), a_TableTemp(m) ) * GQW(l)
          aa_TableIDPFDT(m,n) = aa_TableIDPFDT(m,n) + DPFDT( GQP(l), a_TableTemp(m) ) * GQW(l)
        end do
      end do
    end do
    deallocate( GQP )
    deallocate( GQW )


    !----------------------------------------------------
    ! Check accuracy of integration of Planc function by using a pre-calculated table.
    !

    !      This routine is called once here, to initialize a pre-calculated table. 
    n = 1
    xy_TempTMP = 300.0d0
    call CalcIntegratedPFWithTable2D( n, xy_TempTMP, xy_PFTable, 1, jmax, .false. )

    do n = 1, nbmax

      FlagCheckLoopExit = .false.
      l = 1
      do

        do j = 1, jmax
          do i = 0, imax-1
            xy_TempTMP(i,j) = a_TableTemp(1) + ( a_TableTemp(2) - a_TableTemp(1) ) * 0.5d0 + ( a_TableTemp(2) - a_TableTemp(1) ) * ( imax * jmax * ( l - 1 ) + imax * ( j - 1 ) + i )
          end do
        end do

        do j = 1, jmax
          do i = 0, imax-1
            if ( xy_TempTMP(i,j) > a_TableTemp(ntmax) ) then
              xy_TempTMP(i,j) = a_TableTemp(ntmax)
              FlagCheckLoopExit = .true.
            end if
          end do
        end do


        call CalcIntegratedPF2D( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, imax, jmax, xy_TempTMP, xy_PF )
        call Integ_DPFDT_GQ_Array2D( a_BandParam(1,n), a_BandParam(2,n), NGaussQuad, xy_TempTMP, xy_DPFDT )

        call CalcIntegratedPFWithTable2D( n, xy_TempTMP, xy_PFTable, 1, jmax, .false. )
        call CalcIntegratedPFWithTable2D( n, xy_TempTMP, xy_DPFDTTable, 1, jmax, .true. )

        do j = 1, jmax
          do i = 0, imax-1
            ErrorPFInteg = abs( xy_PF   (i,j) - xy_PFTable   (i,j) ) / xy_PF   (i,j)
            if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
              call MessageNotify( 'E', module_name, 'Error of integrated PF, %f, is greater than threshold, %f, in band %d.', d = (/ ErrorPFInteg, ThresholdErrorPFInteg /), i = (/n/) )
            end if
            ErrorPFInteg = abs( xy_DPFDT(i,j) - xy_DPFDTTable(i,j) ) / xy_DPFDT(i,j)
            if ( ErrorPFInteg > ThresholdErrorPFInteg ) then
              call MessageNotify( 'E', module_name, 'Error of integrated DPFDT, %f, is greater than threshold, %f, in band %d', d = (/ ErrorPFInteg, ThresholdErrorPFInteg /), i = (/n/) )
            end if
          end do
        end do

        if ( FlagCheckLoopExit ) exit
        l = l + 1
      end do

    end do

    !----------------------------------------------------


    do n = 1, nRefAtm
      RefAtmPro(2,n) = RefAtmPro(2,n) * 1.0d2
!!$      RefAtmPro(5,n) = RefAtmPro(5,n) * 1.0d-3 / ( 48.0d0 * 1.66d-27 ) &
!!$        & / ( RefAtmPro(2,n) / ( 1.38d-23 * RefAtmPro(3,n) ) )
      RefAtmPro(5,n) = RefAtmPro(5,n) * 1.0d-3 / ( RefAtmPro(2,n) / ( 8.314d0 / 48.0d-3 * RefAtmPro(3,n) ) )
    end do



    sw_fs = .false.


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  DelTimeCalcTrans  = %f [%c]', d = (/ DelTimeCalcTransValue /), c1 = trim( DelTimeCalcTransUnit ) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


  end subroutine RadiationRoeweLiouInit
RefAtmPro
Variable :
RefAtmPro(1:5, 1:nRefAtm) :real(DP), save
Subroutine :
index :integer , intent(in )
sw_log :logical , intent(in )
xyz_Press(0:imax-1,1:jmax,1:kmax) :real(DP), intent(in )
xyz_Array(0:imax-1,1:jmax,1:kmax) :real(DP), intent(out)

[Source]

  subroutine SetRefAtmPro( index, sw_log, xyz_Press, xyz_Array )

    use dc_types
    use gridset, only: imax, jmax, kmax

    integer , intent(in ) :: index
    logical , intent(in ) :: sw_log
    real(DP), intent(in ) :: xyz_Press(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_Array(0:imax-1,1:jmax,1:kmax)

    !
    ! local variables
    !
    real(DP):: z_RefAtmProArray(1:nRefAtm)
    integer :: k, kk


    if ( sw_log ) then
      do k = 1, nRefAtm
        z_RefAtmProArray(k) = log( RefAtmPro(index,k) + 1.0d-100 )
      end do
    else
      do k = 1, nRefAtm
        z_RefAtmProArray(k) =      RefAtmPro(index,k)
      end do
    end if


    do k = 1, kmax
      if( xyz_Press(0,1,k) <= RefAtmPro(2,nRefAtm) ) then
        xyz_Array(0,1,k) = z_RefAtmProArray(k)
      else
        search_loop : do kk = 2, nRefAtm
          if( RefAtmPro(2,kk) < xyz_Press(0,1,k) ) exit search_loop
        end do search_loop
        if( kk > nRefAtm ) stop 'Unexpected error in setting temperature profile'
        xyz_Array(0,1,k) = ( z_RefAtmProArray( kk ) - z_RefAtmProArray( kk-1 ) ) / ( log( RefAtmPro(2,kk)  / RefAtmPro(2,kk-1) ) ) * ( log( xyz_Press(0,1,k) / RefAtmPro(2,kk-1) ) ) + z_RefAtmProArray( kk-1 )
      end if
    end do


    if ( sw_log ) then
      do k = 1, kmax
        xyz_Array(0,1,k) = exp( xyz_Array(0,1,k) )
      end do
    else
      do k = 1, kmax
        xyz_Array(0,1,k) =      xyz_Array(0,1,k)
      end do
    end if

    do k = 1, kmax
      xyz_Array(:,:,k) = xyz_Array(0,1,k)
    end do


  end subroutine SetRefAtmPro
TableTempIncrement
Variable :
TableTempIncrement :real(DP), save
TableTempMax
Variable :
TableTempMax :real(DP), save
TableTempMin
Variable :
TableTempMin :real(DP), save
VMRCO2
Variable :
VMRCO2 :real(DP), save
: Volume mixing ratio of CO2
a_BandParam
Variable :
a_BandParam(1:4,1:nbmax) :real(DP), save
a_CO2BandParam
Variable :
a_CO2BandParam(1:4) :real(DP), save
a_TableTemp
Variable :
a_TableTemp(:) :real(DP), save, allocatable
aa_TableIDPFDT
Variable :
aa_TableIDPFDT(:,:) :real(DP), save, allocatable
aa_TableIPF
Variable :
aa_TableIPF(:,:) :real(DP), save, allocatable
flag_save_time
Variable :
flag_save_time :logical , save
Subroutine :
x1 :real(DP), intent(in )
x2 :real(DP), intent(in )
n :integer , intent(in )
x(n) :real(DP), intent(out)
w(n) :real(DP), intent(out)

[Source]

  subroutine gauleg( x1, x2, n, x, w )

    ! 
    ! Physical constants settings
    !
    use constants, only: PI   ! $ \pi $ .
                              ! Circular constant

    real(DP), intent(in ) :: x1,x2
    integer , intent(in ) :: n
    real(DP), intent(out) :: x(n),w(n)


    real(DP):: eps
    real(DP):: z1, z, xm, xl, pp, p3, p2, p1
    integer :: i, j, m


    !
    ! changed at 2005/09/14
    !
!!$      eps = 1.0e-11
    eps = 1.0e-15

    m  = ( n+1 ) / 2
    xm = ( x2+x1 ) / 2.0d0
    xl = ( x2-x1 ) / 2.0d0

    do i = 1, m
      z = cos( pi*dble(i-0.25)/(dble(n+0.5)) )

100   continue
      p1 = 1.0d0
      p2 = 0.0d0
!!$         do j = 1, n
      do j = 1, n
        p3 = p2
        p2 = p1
        p1 = ( (2.0d0*dble(j)-1.0d0)*z*p2-(dble(j)-1.0d0)*p3 ) / dble(j)
      enddo

      pp = dble(n) * (z*p1-p2) / (z*z-1.0d0)
      z1 = z
      z  = z1 - p1 / pp

      !
      ! changed at 2005/09/14
      !
!!$         if( abs( z-z1 ) .gt. eps ) go to 100
!!$         if( abs( z-z1 ) / z1 .gt. eps ) go to 100
      !
      ! changed before 2008/08/03
      !
      if( abs( z-z1 ) / abs( z1 + 1.0d-200 ) .gt. eps ) go to 100

      x( i     ) = xm-xl*z
      x( n+1-i ) = xm+xl*z
      w( i     ) = 2.0d0*xl/((1.0d0-z*z)*pp*pp)
      w( n+1-i ) = w(i)
    end do


  end subroutine gauleg
module_name
Constant :
module_name = ‘radiation_RL78 :character(*), parameter
: モジュールの名称. Module name
nRefAtm
Constant :
nRefAtm = 31 :integer , parameter
nbmax
Constant :
nbmax = 32 :integer , parameter
ntmax
Variable :
ntmax :integer , save
Function :
pf :real(DP)
x :real(DP), intent(in)
t :real(DP), intent(in)

[Source]

  real(DP) function pf(x,t)

    real(DP), intent(in) :: x, t


    real(dp), parameter :: c      = 2.99792458d8           , planc  = 6.6260755d-34          , boltz  = 1.380658d-23


    pf=2.0d0*planc*c*c*x*x*x /(exp(planc*c*(x+1.0d-10)/(boltz*t))-1.0d0)

    return
  end function pf
sw_fs
Variable :
sw_fs :logical , save
version
Constant :
version = ’$Name: dcpam5-20100424 $’ // ’$Id: radiation_RL78.f90,v 1.6 2010-04-13 13:14:45 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version
xyrra_TransSave
Variable :
xyrra_TransSave(:,:,:,:,:) :real(DP), allocatable, save