Class radiation_LH74
In: radiation/radiation_LH74.f90

Lacis and Hansen (1974) による短波放射モデル

Short wave radiation model described by Lacis and Hansen (1974)

Note that Japanese and English are described in parallel.

短波放射モデル.

This is a model of short wave radiation.

References

 Lacis, A. A., and J. E. Hansen, A parameterization for the absorption of solar
   radiation in the Earth's atmosphere, J. Atmos. Sci., 31, 118-133, 1974.

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 gridset constants radiation_short_income read_time_series namelist_util dc_iounit dc_message

Public Instance methods

Subroutine :
xyz_QO3(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xy_SurfAlbedo(0:imax-1, 1:jmax) :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_DTempDtRadS(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(inout)

Lacis and Hansen (1974) の計算方を用いて O3 の短波吸収加熱率を計算し, 短波放射加熱率に加える.

Calculate radiative heating rate due to absorption of short wave radiation of O3 and add to the short wave radiative heating rate

[Source]

  subroutine RadiationLH74AddO3Heating( xyz_QO3, xy_SurfAlbedo, xyr_Press, xyz_Press, xyz_DTempDtRadS )
    !
    ! Lacis and Hansen (1974) の計算方を用いて O3 の短波吸収加熱率を計算し, 
    ! 短波放射加熱率に加える. 
    !
    ! Calculate radiative heating rate due to absorption of short wave radiation of O3
    ! and add to the short wave radiative heating rate
    !


    ! USE statements
    !

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

    ! 
    ! Physical constants settings
    !
    use constants, only: Grav, PI, CpDry   ! $ C_p $ [J kg-1 K-1].
                                 ! 乾燥大気の定圧比熱.
                                 ! Specific heat of air at constant pressure

    ! 短波入射 (太陽入射)
    ! Short wave (insolation) incoming
    !
    use radiation_short_income, only: ShortIncoming
!    use radiation_short_income_sr, only: ShortIncoming

!!$  ! 時刻管理
!!$  ! Time control
!!$  !
!!$  use timeset, only: &
!!$    & TimeN                   ! ステップ $ t $ の時刻.
!!$                              ! Time of step $ t $.

    use read_time_series, only: SetValuesFromTimeSeriesWrapper

!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

    real(DP), intent(in   ):: xyz_QO3        (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ):: xy_SurfAlbedo  (0:imax-1, 1:jmax)
    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(inout):: xyz_DTempDtRadS(0:imax-1, 1:jmax, 1:kmax)


    !
    ! Work variables
    !
    real(DP):: xyz_O3DelAbsAmt (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyr_O3ColDen    (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_O3AbsAmt    (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_RadSFlux    (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xy_IncomRadSFlux(0:imax-1, 1:jmax)
                              ! 短波 (日射) フラックス.
                              ! Short wave (insolation) flux
    real(DP):: xy_InAngle      (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP):: xy_MagFac       (0:imax-1, 1:jmax)

!!$    logical :: flag_dry_atmosphere

    character(STRING):: O3FileName

    integer :: i
    integer :: j
    integer :: k


    if ( .not. radiation_lh74_inited ) then
      call RadiationLH74Init
    end if


!!$    write( O3FileName, '(a,i3.3,a)' ) &
!!$      & "../../../data_Earth/O3_CMIP5_climatology_zonalmean_T", (imax-1)/3, ".nc"
!!$    call SetValuesFromTimeSeriesWrapper( &
!!$      & O3FileName, "O3", &
!!$      & xyz_Press,         &               ! (in)
!!$      & xyz_QO3,           &               ! (inout)
!!$      & "O3"            &
!!$      & )
!!$
!!$    call HistoryAutoPut( TimeN, "O3", xyz_QO3 )


    ! 短波入射の計算
    ! Calculate short wave (insolation) incoming radiation
    !
    call ShortIncoming( xy_IncomRadSFlux, xy_InAngle )


    do k = 1, kmax
      xyz_O3DelAbsAmt(:,:,k) = xyz_QO3(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do

    xyr_O3ColDen(:,:,:) = 0.0d0
    do k = kmax-1, 0, -1
      xyr_O3ColDen(:,:,k) = xyr_O3ColDen(:,:,k+1) + xyz_O3DelAbsAmt(:,:,k+1)
    end do

    if ( FlagSimpleMagFac ) then

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = xy_InAngle(i,j)
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    else

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = 35.0d0 / sqrt( 1224.0d0 * ( 1.0d0 / xy_InAngle(i,j) )**2 + 1.0d0 )
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    end if


    xyr_RadSFlux(:,:,:) = 0.0d0

    ! Downward flux
    do k = 0, kmax
      xyr_O3AbsAmt(:,:,k) = xyr_O3ColDen(:,:,k) * xy_MagFac(:,:)
    end do
    call RadiationLH74CalcO3Absorptivity( xyr_O3AbsAmt, xyr_Absorptivity )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) + xy_IncomRadSFlux(:,:) * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do

    ! Upward flux
    do k = 0, kmax
      xyr_O3AbsAmt(:,:,k) = xyr_O3ColDen(:,:,0) * xy_MagFac(:,:) + ( xyr_O3ColDen(:,:,0) - xyr_O3ColDen(:,:,k) ) * DiffFactorO3
    end do
    call RadiationLH74CalcO3Absorptivity( xyr_O3AbsAmt, xyr_Absorptivity )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) - xy_IncomRadSFlux(:,:) * xy_SurfAlbedo(:,:) * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do

    ! 放射加熱率の計算
    ! Calculate radiation heating rate
    !
    do k = 1, kmax
      xyz_DTempDtRadS(:,:,k) = xyz_DTempDtRadS(:,:,k) + (     xyr_RadSFlux(:,:,k-1) - xyr_RadSFlux(:,:,k) ) / ( xyr_Press(:,:,k-1)    - xyr_Press(:,:,k) ) / CpDry * Grav
    end do

!!$    i = imax / 2
!!$    j = jmax / 2
!!$    j = jmax * 3 / 4
!!$    do k = 0, kmax
!!$      write( 93, * ) xyr_RadSFlux(i,j,k), xyr_Press(i,j,k)
!!$    end do
!!$    call flush( 93 )
!!$
!!$    do k = 1, kmax
!!$      write( 83, * ) &
!!$        & + (     xyr_RadSFlux(i,j,k-1) - xyr_RadSFlux(i,j,k) )  &
!!$        &     / ( xyr_Press(i,j,k-1)    - xyr_Press(i,j,k) )     &
!!$        &     / CpDry * Grav, &
!!$        & xyz_Press(i,j,k)
!!$    end do
!!$    call flush( 83 )
!!$
!!$
!!$    stop

  end subroutine RadiationLH74AddO3Heating
Subroutine :
xy_SurfAlbedo(0:imax-1, 1:jmax) :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 )
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 )
xyr_RadSFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)

[Source]

  subroutine RadiationLH74Flux( xy_SurfAlbedo, xyz_Temp, xyz_QVap, xyr_Press, xyz_Press, xyr_RadSFlux )


    ! 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

    ! 短波入射 (太陽入射)
    ! Short wave (insolation) incoming
    !
    use radiation_short_income, only: ShortIncoming

    real(DP), intent(in ):: xy_SurfAlbedo   (0:imax-1, 1:jmax)
    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 ):: 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(out):: xyr_RadSFlux    (0:imax-1, 1:jmax, 0:kmax)


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

    real(DP):: xyz_H2ODelAbsAmt(0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyr_H2OColDen   (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_H2OAbsAmt   (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xy_IncomRadSFlux(0:imax-1, 1:jmax)
                              ! 短波 (日射) フラックス.
                              ! Short wave (insolation) flux
    real(DP):: xy_InAngle      (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    real(DP):: xy_MagFac       (0:imax-1, 1:jmax)

!!$    logical :: flag_dry_atmosphere

    integer :: i
    integer :: j
    integer :: k


    if ( .not. radiation_lh74_inited ) then
      call RadiationLH74Init
    end if


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


    ! 短波入射の計算
    ! Calculate short wave (insolation) incoming radiation
    !
    call ShortIncoming( xy_IncomRadSFlux, xy_InAngle )

    ! 大気アルベドの考慮
    ! Taking atmospheric albedo into consideration
    !
    xy_IncomRadSFlux = xy_IncomRadSFlux * ( 1.0d0 - ShortAtmosAlbedo )

!!$    if ( flag_dry_atmosphere ) then
!!$      do k = 0, kmax
!!$        xyr_RadSFlux(:,:,k) = - xy_IncomRadSFlux(:,:) + ...
!!$      end do
!!$      return
!!$    end if


    RefPress = 1.013d5
    RefTemp  = 273.0d0

    do k = 1, kmax
      xyz_H2ODelAbsAmt(:,:,k) = ( xyz_Press(:,:,k) / RefPress        )**H2OScaleIndex * ( RefTemp          / xyz_Temp(:,:,k) )**0.5 * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
!!$      xyz_H2ODelAbsAmt(:,:,k) = &
!!$        &   ( xyz_Press(:,:,k) / RefPress ) &
!!$        & * xyz_QVap(:,:,k) * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do


    xyr_H2OColDen(:,:,:) = 0.0d0

    do k = kmax-1, 0, -1
      xyr_H2OColDen(:,:,k) = xyr_H2OColDen(:,:,k+1) + xyz_H2ODelAbsAmt(:,:,k+1)
    end do


    if ( FlagSimpleMagFac ) then

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = xy_InAngle(i,j)
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    else

      do j = 1, jmax
        do i = 0, imax-1
          if ( xy_InAngle(i,j) > 0.0d0 ) then
            xy_MagFac(i,j) = 35.0d0 / sqrt( 1224.0d0 * ( 1.0d0 / xy_InAngle(i,j) )**2 + 1.0d0 )
          else
            xy_MagFac(i,j) = 0.0d0
          end if
        end do
      end do

    end if


    xyr_RadSFlux(:,:,:) = 0.0d0

    ! Downward flux
    do k = 0, kmax
      xyr_H2OAbsAmt(:,:,k) = xyr_H2OColDen(:,:,k) * xy_MagFac(:,:)
    end do
    call RadiationLH74CalcAbsorptivity( xyr_H2OAbsAmt, xyr_Absorptivity )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) + xy_IncomRadSFlux(:,:) * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do

    ! Upward flux
    do k = 0, kmax
      xyr_H2OAbsAmt(:,:,k) = xyr_H2OColDen(:,:,0) * xy_MagFac(:,:) + ( xyr_H2OColDen(:,:,0) - xyr_H2OColDen(:,:,k) ) * DiffFactorH2O
    end do
    call RadiationLH74CalcAbsorptivity( xyr_H2OAbsAmt, xyr_Absorptivity )
    do k = 0, kmax
      xyr_RadSFlux(:,:,k) = xyr_RadSFlux(:,:,k) - xy_IncomRadSFlux(:,:) * xy_SurfAlbedo(:,:) * ( 1.0d0 - xyr_Absorptivity(:,:,k) )
    end do


  end subroutine RadiationLH74Flux

Private Instance methods

DiffFactorH2O
Variable :
DiffFactorH2O :real(DP), save
DiffFactorO3
Variable :
DiffFactorO3 :real(DP), save
FlagSimpleMagFac
Variable :
FlagSimpleMagFac :logical , save
: エアマス計算のためのフラグ Flag for air-mass calculation
H2OScaleIndex
Variable :
H2OScaleIndex :real(DP), save
Subroutine :
xyr_H2OAbsAmt(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)

[Source]

  subroutine RadiationLH74CalcAbsorptivity( xyr_H2OAbsAmt, xyr_Absorptivity )

    ! USE statements
    !

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

    real(DP), intent(in ):: xyr_H2OAbsAmt   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)

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


    xyr_H2Oprcm(:,:,:) = xyr_H2OAbsAmt(:,:,:) / ( 1.0d0 * 1.0d-3 * 1.0d6 ) * 1.0d2

    xyr_Absorptivity(:,:,:) = 2.9d0 * xyr_H2Oprcm(:,:,:) / ( ( 1.0d0 + 141.5d0 * xyr_H2Oprcm(:,:,:) )**0.635 + 5.925d0 * xyr_H2Oprcm(:,:,:) )


  end subroutine RadiationLH74CalcAbsorptivity
Subroutine :
xyr_O3AbsAmt(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(in )
xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)

[Source]

  subroutine RadiationLH74CalcO3Absorptivity( xyr_O3AbsAmt, xyr_Absorptivity )

    ! USE statements
    !

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: GasRUniv
                              ! $ R^{*} $ [J K-1 mol-1].
                              ! 普遍気体定数.  Universal gas constant

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

    real(DP), intent(in ):: xyr_O3AbsAmt    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyr_Absorptivity(0:imax-1, 1:jmax, 0:kmax)

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


    O3DensNTP = 1013.25d2 / ( GasRUniv / 48.0d-3 * 273.15d0 )

    xyr_O3cm(:,:,:) = xyr_O3AbsAmt(:,:,:) / O3DensNTP * 1.0d2

    xyr_Absorptivity(:,:,:) = 0.02118d0 * xyr_O3cm(:,:,:) / (   1.0d0 + 0.042d0 * xyr_O3cm(:,:,:) + 0.000323d0 * xyr_O3cm(:,:,:) * xyr_O3cm(:,:,:) ) + 1.082d0 * xyr_O3cm(:,:,:) / (   1.0d0 + 138.6d0 * xyr_O3cm(:,:,:) )**0.805 + 0.0658d0 * xyr_O3cm(:,:,:) / (   1.0d0 + ( 103.6d0 * xyr_O3cm(:,:,:) )**3 )


  end subroutine RadiationLH74CalcO3Absorptivity
Subroutine :

This procedure input/output NAMELIST#radiation_LH74_nml .

[Source]

  subroutine RadiationLH74Init


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

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

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

!!$    ! ヒストリデータ出力
!!$    ! History data output
!!$    !
!!$    use gtool_historyauto, only: HistoryAutoAddVariable


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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /radiation_LH74_nml/ DiffFactorH2O, DiffFactorO3, ShortAtmosAlbedo, FlagSimpleMagFac
          !
          ! デフォルト値については初期化手続 "radiation_LH74#RadiationLH74Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_LH74#RadiationLH74Init" for the default values.
          !

    ! デフォルト値の設定
    ! Default values settings
    !
    DiffFactorH2O    = 5.0d0 / 3.0d0
    DiffFactorO3     = 1.9_DP

    ShortAtmosAlbedo = 0.2_DP

    FlagSimpleMagFac = .false.

    ! NAMELIST の読み込み
    ! 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_LH74_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

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


    H2OScaleIndex = 1.0d0


!!$    call HistoryAutoAddVariable( "O3",         & ! (in)
!!$      & (/ 'lon ', 'lat ', 'sig ', 'time' /),  & ! (in)
!!$      & "ozone", 'kg kg-1' )                     ! (in)


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'DiffFactorH2O    = %f', d = (/ DiffFactorH2O /) )
    call MessageNotify( 'M', module_name, 'DiffFactorO3     = %f', d = (/ DiffFactorO3 /) )
    call MessageNotify( 'M', module_name, 'ShortAtmosAlbedo = %f', d = (/ ShortAtmosAlbedo /) )
    call MessageNotify( 'M', module_name, 'FlagSimpleMagFac = %b', l = (/ FlagSimpleMagFac /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    radiation_lh74_inited = .true.
  end subroutine RadiationLH74Init
ShortAtmosAlbedo
Variable :
ShortAtmosAlbedo :real(DP), save
: 大気アルベド. Albedo of air.
module_name
Constant :
module_name = ‘radiation_LH74 :character(*), parameter
: モジュールの名称. Module name
radiation_lh74_inited
Variable :
radiation_lh74_inited :logical , save
version
Constant :
version = ’$Name: dcpam5-20101008-1 $’ // ’$Id: radiation_LH74.f90,v 1.5 2010-04-12 02:29:06 noda Exp $’ :character(*), parameter
: モジュールのバージョン Module version