Class rad_Earth_LW_V2_4
In: radiation/rad_Earth_LW_V2_4.f90

地球大気向け長波放射モデル Ver. 2.4

long wave radiation model for the Earth‘s atmosphere Ver. 2.4

Note that Japanese and English are described in parallel.

長波放射モデル.

This is a model of long wave radiation for the Earth‘s atmospehre. Radiation in the wavenumber range from 0 to 3000 cm-1 is calculated following the scheme by Chou et al. (2001). But absorptions by CFC and weak bands desinated as band 10 in Chou et al. (2001) are neglected.

References

 Chou, M.-D., M. J. Suarez, X.-Z. Liang, and M. M.-H. Yan,
   A thermal infrared radiation parameterization for atmospheric studies,
   NASA Technical Report Series on Global Modeling and Data Assimilation,
   19, NASA/TM-2001-104606, 2001.

Procedures List

RadEarthLWV24Flux :放射フラックスの計算
———— :————
RadEarthLWV24Flux :Calculate radiation flux

NAMELIST

NAMELIST#rad_Earth_LW_V2_4_nml

Methods

Included Modules

dc_types constants0 gridset dc_message timeset rad_CK1991 rad_C2001 rad_rte_nonscat cloud_utils dc_iounit dc_calendar namelist_util gauss_quad planck_func

Public Instance methods

Subroutine :
xyz_DelCO2Mass(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelH2OVapMass(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelH2OLiqMass(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelH2OSolMass(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelO3Mass(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelN2OMass(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_DelCH4Mass(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_Temp(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xy_SurfTemp(0:imax-1, 1:jmax) :real(DP), intent(in )
xyz_QCO2(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QH2OVap(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QN2O(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_QCH4(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_CloudCover(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_CloudWatREff(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyz_CloudIceREff(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyr_RadLUwFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
xyr_RadLDwFlux(0:imax-1, 1:jmax, 0:kmax) :real(DP), intent(out)
xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) :real(DP), intent(out)
xyra_DelRadLDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) :real(DP), intent(out)

[Source]

  subroutine RadEarthLWV24Flux( xyz_DelCO2Mass, xyz_DelH2OVapMass, xyz_DelH2OLiqMass, xyz_DelH2OSolMass, xyz_DelO3Mass, xyz_DelN2OMass, xyz_DelCH4Mass, xyz_Press, xyz_Temp, xy_SurfTemp, xyz_QCO2, xyz_QH2OVap, xyz_QN2O, xyz_QCH4, xyz_CloudCover, xyz_CloudWatREff, xyz_CloudIceREff, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )


    ! USE statements
    !

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

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


!!$    ! Chou et al (1991) による長波放射モデル
!!$    ! Long radiation model described by Chou et al (1991)
!!$    !
!!$    use rad_C1991, only :               &
!!$      & RadC1991CalcTransMAH2O

    ! Chou and Kouvaris (1991) による長波放射モデル
    ! Long radiation model described by Chou and Kouvaris (1991)
    !
    use rad_CK1991, only : RadCK1991CalcTrans

    ! Chou et al. (2001) による長波放射モデル
    ! Long radiation model described by Chou et al. (2001)
    !
    use rad_C2001, only : RadC2001CalcTransBand3CO2, RadC2001CalcTransBand3H2O, RadC2001CalcTrans, RadC2001ReduceCloudOptDep, RadC2001CalcCloudOptProp , RadC2001CalcIntegratedPF2D, RadC2001CalcIntegratedPF3D

    ! 散乱を無視した放射伝達方程式
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScat !, RadRTENonScatAnotherForm

    ! 雲関系ルーチン
    ! Cloud-related routines
    !
    use cloud_utils, only : CloudUtilsLocalizeCloud, CloudUtilsCalcOverlapCloudTrans


    real(DP), intent(in ):: xyz_DelCO2Mass   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelH2OVapMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelH2OLiqMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelH2OSolMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelO3Mass    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelN2OMass   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_DelCH4Mass   (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_Temp         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xy_SurfTemp      (0:imax-1, 1:jmax)
    real(DP), intent(in ):: xyz_QCO2         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QH2OVap      (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QN2O         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_QCH4         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudCover   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudWatREff (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ):: xyz_CloudIceREff (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyr_RadLUwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyr_RadLDwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out):: xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP), intent(out):: xyra_DelRadLDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


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

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

    real(DP) :: xyz_CloudWatSSA      (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudIceSSA      (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudWatAF       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_CloudIceAF       (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_DelCloudWatOptDep(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelCloudIceOptDep(0:imax-1, 1:jmax, 1: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)
    real(DP) :: xyrr_TransEach        (0:imax-1, 1:jmax, 0:kmax, 0:kmax)

    real(DP) :: xyz_IntPF   (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_SurfIntPF(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT0(0:imax-1, 1:jmax)
    real(DP) :: xy_IntDPFDT1(0:imax-1, 1:jmax)


!!$    real(DP) :: xyr_RadFluxMA    (0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP) :: xyra_DelRadFluxMA(0:imax-1, 1:jmax, 0:kmax, 0:1)


    real(DP) :: xyr_RadUwFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_RadDwFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyra_DelRadUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP) :: xyra_DelRadDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)


    integer  :: n


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. rad_Earth_LW_V2_4_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


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


    ! Check whether the transmittance is saved or not
    !
    if (  .not. FlagTransSaved ) then
      if ( TimeN - PrevTimeSave >= IntTimeSave ) then
        call MessageNotify( 'M', module_name, 'Transmittance is not saved, but criterion for transmittance calculation is met.' )
      else
        call MessageNotify( 'M', module_name, 'Transmittance is not saved, and criterion for transmittance calculation ' // 'is not met. However, transmittance will be calculated.' )
      end if
    end if


    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.


      LOOP_band_trans: do n = 1, nbmax

        xyrr_Trans = 1.0_DP

        if ( n == nbmax ) then

          ! Now, nothing is done when n = nbmax.

        else if ( n == 3 ) then
          ! 540-800 cm-1

          !   Calculation of H2O line and continuum transmittance
          call RadC2001CalcTransBand3H2O( xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransEach )
          xyrr_Trans = xyrr_Trans * xyrr_TransEach
          !   Calculation of CO2 transmittance
          if ( FlagHighAlt ) then
            ! Transmittance calculation for middle atmospehre as well as lower atmosphere
            call RadCK1991CalcTrans( xyz_DelCO2Mass, xyz_Press, xyz_Temp, 'CO2', xyrr_TransEach )
          else
            ! Transmittance calculation for lower atmoshere
            call RadC2001CalcTransBand3CO2( xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2, xyrr_TransEach )
          end if
          xyrr_Trans = xyrr_Trans * xyrr_TransEach

        else

          !   Calculation of H2O continuum transmittance
          call RadC2001CalcTrans( 'H2OCont', n, xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransEach )
          xyrr_Trans = xyrr_Trans * xyrr_TransEach
          !   Calculation of H2O line transmittance
          call RadC2001CalcTrans( 'H2OLine', n, xyz_DelH2OVapMass, xyz_Press, xyz_Temp, xyz_QH2OVap, xyrr_TransEach )
          xyrr_Trans = xyrr_Trans * xyrr_TransEach

          ! Calculation of O3 transmittance
          if ( n == 5 ) then
            ! 980-1100 cm-1
            call RadCK1991CalcTrans( xyz_DelO3Mass, xyz_Press, xyz_Temp, 'O3', xyrr_TransEach )
            xyrr_Trans = xyrr_Trans * xyrr_TransEach
          end if

          ! Calculation of N2O transmittance
          if ( ( n == 6 ) .or. ( n == 7 ) .or. ( n == 10 ) ) then
            if ( maxval( xyz_DelN2OMass ) > 0.0_DP ) then
              call RadC2001CalcTrans( 'N2OLine', n, xyz_DelN2OMass, xyz_Press, xyz_Temp, xyz_QN2O, xyrr_TransEach )
              xyrr_Trans = xyrr_Trans * xyrr_TransEach
            end if
          end if

          ! Calculation of CH4 transmittance
          if ( ( n == 6 ) .or. ( n == 7 ) ) then
            if ( maxval( xyz_DelCH4Mass ) > 0.0_DP ) then
              call RadC2001CalcTrans( 'CH4Line', n, xyz_DelCH4Mass, xyz_Press, xyz_Temp, xyz_QCH4, xyrr_TransEach )
              xyrr_Trans = xyrr_Trans * xyrr_TransEach
            end if
          end if

          ! Calculation of CO2 transmittance in weak bands
          if ( ( n == 4 ) .or. ( n == 5 ) .or. ( n == 10 ) ) then
            call RadC2001CalcTrans( 'WeakBandCO2Line', n, xyz_DelCO2Mass, xyz_Press, xyz_Temp, xyz_QCO2, xyrr_TransEach )
            xyrr_Trans = xyrr_Trans * xyrr_TransEach
          end if

        end if

        xyrra_TransSaved(:,:,:,:,n) = xyrr_Trans

      end do LOOP_band_trans

      !
      !   Calculation of transmittance of water vapor by using a method for middle 
      !   atmosphere
      !
!!$      call RadC1991CalcTransMAH2O(      &
!!$        & xyr_Press, xyz_Press, xyz_Temp, xyz_QH2OVap, & ! (in)
!!$        & xyrr_Trans                                & ! (out)
!!$        & )
      xyrr_Trans = -1.0d100

      xyrr_TransMASaved = xyrr_Trans

    end if



    !
    ! Calculate radiative flux
    !

    xy_SurfEmis = 1.0_DP

    xyr_RadLUwFlux     = 0.0_DP
    xyr_RadLDwFlux     = 0.0_DP
    xyra_DelRadLUwFlux = 0.0_DP
    xyra_DelRadLDwFlux = 0.0_DP


    LOOP_band_RTE : do n = 1, nbmax


      xyz_CloudREff = xyz_CloudWatREff
      call RadC2001CalcCloudOptProp( 'Liquid', n, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudWatSSA, xyz_CloudWatAF )
      xyz_DelCloudWatOptDep = xyz_CloudExtCoef * xyz_DelH2OLiqMass
      !
      xyz_CloudREff = xyz_CloudIceREff
      call RadC2001CalcCloudOptProp( 'Ice', n, xyz_CloudREff, xyz_CloudExtCoef, xyz_CloudIceSSA, xyz_CloudIceAF )
      xyz_DelCloudIceOptDep = xyz_CloudExtCoef * xyz_DelH2OSolMass

      call RadC2001ReduceCloudOptDep( xyz_CloudWatSSA, xyz_CloudWatAF, xyz_DelCloudWatOptDep )
      call RadC2001ReduceCloudOptDep( xyz_CloudIceSSA, xyz_CloudIceAF, xyz_DelCloudIceOptDep )
      !
      call CloudUtilsLocalizeCloud( xyz_CloudCover, xyz_DelCloudWatOptDep )
      call CloudUtilsLocalizeCloud( xyz_CloudCover, xyz_DelCloudIceOptDep )
      !
      xyz_TransCloudOneLayer = exp( - ( xyz_DelCloudWatOptDep + xyz_DelCloudIceOptDep ) * DiffFactor )
      !
      call CloudUtilsCalcOverlapCloudTrans( xyz_TransCloudOneLayer, xyz_CloudCover, xyrr_TransCloud )


      ! Now, nothing is done when n = nbmax.
      if ( n == nbmax ) cycle


      xyrr_Trans = xyrra_TransSaved(:,:,:,:,n)

      xyrr_Trans = xyrr_Trans * xyrr_TransCloud



      call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_SurfIntPF, 1, jmax )
      call CalcIntegratedPFWithTable3D( n, kmax, xyz_Temp, xyz_IntPF, 1, jmax )

      call CalcIntegratedPFWithTable2D( n, xy_SurfTemp, xy_IntDPFDT0, 1, jmax, .true. )
      call CalcIntegratedPFWithTable2D( n, xyz_Temp(:,:,1), xy_IntDPFDT1, 1, jmax, .true. )

      xy_SurfIntPF = xy_SurfEmis * PI * xy_SurfIntPF
      xyz_IntPF    =               PI * xyz_IntPF
      xy_IntDPFDT0 = xy_SurfEmis * PI * xy_IntDPFDT0
      xy_IntDPFDT1 =               PI * xy_IntDPFDT1



      ! Lines below are under testing.
      !
!!$      xy_SurfIntPF2 = xy_SurfIntPF
!!$      xyz_IntPF2    = xyz_IntPF
!!$      xy_IntDPFDT02 = xy_IntDPFDT0
!!$      xy_IntDPFDT12 = xy_IntDPFDT1
!!$
!!$
!!$      call RadC2001CalcIntegratedPF2D(       &
!!$        & n, xy_SurfTemp,                    &
!!$        & xy_SurfIntPF                       &
!!$        & )
!!$      call RadC2001CalcIntegratedPF3D(       &
!!$        & n, kmax, xyz_Temp,                 &
!!$        & xyz_IntPF                          &
!!$        & )
!!$
!!$      call RadC2001CalcIntegratedPF2D(       &
!!$        & n, xy_SurfTemp,                    &
!!$        & xy_IntDPFDT0,                      &
!!$        & .true.                             &
!!$        & )
!!$      call RadC2001CalcIntegratedPF2D(       &
!!$        & n, xyz_Temp(:,:,1),                &
!!$        & xy_IntDPFDT1,                      &
!!$        & .true.                             &
!!$        & )
!!$
!!$      xy_SurfIntPF = xy_SurfEmis * xy_SurfIntPF
!!$      xyz_IntPF    =               xyz_IntPF
!!$      xy_IntDPFDT0 = xy_SurfEmis * xy_IntDPFDT0
!!$      xy_IntDPFDT1 =               xy_IntDPFDT1
!!$
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          write( 20+n, * ) xy_SurfTemp(i,j), xy_SurfIntPF2(i,j), xy_SurfIntPF(i,j), &
!!$            & xy_IntDPFDT02(i,j), xy_IntDPFDT02(i,j)
!!$        end do
!!$      end do
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          write( 40+n, * ) xyz_Temp(i,j,1), xy_IntDPFDT12(i,j), xy_IntDPFDT12(i,j)
!!$        end do
!!$      end do
!!$      do k = 1, kmax
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            write( 60+n, * ) xyz_Temp(i,j,k), xyz_IntPF2(i,j,k), xyz_IntPF(i,j,k)
!!$          end do
!!$        end do
!!$      end do
!!$      call flush( 20+n )
!!$      call flush( 40+n )
!!$      call flush( 60+n )



!!$      call OLD_RadRTENonScat(                                      &
!!$        & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0, & ! (in)
!!$        & xyrr_Trans,                                          & ! (in)
!!$        & xyr_RadFlux, xyra_DelRadFlux                         & ! (out)
!!$        & )

      call RadRTENonScat( xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_IntDPFDT0, xyrr_Trans, xyr_RadUwFlux, xyr_RadDwFlux, xyra_DelRadUwFlux, xyra_DelRadDwFlux )



!!$      i = 0
!!$      j = jmax/2+1
!!$      do k = 0, kmax
!!$        write( 6, * ) k, xyr_RadFlux(i,j,k), &
!!$          & xyr_RadFlux(i,j,k) - ( xyr_RadUwFlux(i,j,k) - xyr_RadDwFlux(i,j,k) )
!!$      end do
!!$      do k = 0, kmax
!!$        write( 6, * ) k, &
!!$          & xyra_DelRadFlux(i,j,k,0), &
!!$          & xyra_DelRadFlux(i,j,k,0) - ( xyra_DelRadUwFlux(i,j,k,0) - xyra_DelRadDwFlux(i,j,k,0) ), &
!!$          & xyra_DelRadFlux(i,j,k,1), &
!!$          & xyra_DelRadFlux(i,j,k,1) - ( xyra_DelRadUwFlux(i,j,k,1) - xyra_DelRadDwFlux(i,j,k,1) )
!!$      end do



!!$      if ( ( n == 1 ) .or. ( n == 2 ) .or. ( n == 9 ) ) then
!!$        !
!!$        ! For bands 0-340, 340-540, 1380-1900
!!$        ! merge with flux calculated with a method for middle atmosphere
!!$        !
!!$
!!$        xyrr_Trans = xyrr_TransMASaved
!!$
!!$        call RadELWV22IntegRTE(                               &
!!$          & n,                                                & ! (in )
!!$          & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans,   & ! (in )
!!$          & xyr_RadFluxMA, xyra_DelRadFluxMA                  & ! (out)
!!$          & )
!!$
!!$        call RadDcpamELWV22CutMergeFlux(          &
!!$          & xyr_Press,                            & ! (in)
!!$          & xyr_RadFlux, xyra_DelRadFlux,         & ! (inout)
!!$          & xyr_RadFluxMA, xyra_DelRadFluxMA      & ! (in) optional
!!$          & )
!!$
!!$      else if ( ( n == 4 ) .or. ( n == 6 ) .or. ( n == 8 ) ) then
!!$        !
!!$        ! For bands 800-980, 1100-1380, 1900-3000
!!$        ! flux above a pressure level is modified to be constant
!!$        !
!!$
!!$        call RadDcpamELWV22CutMergeFlux(          &
!!$          & xyr_Press,                            & ! (in)
!!$          & xyr_RadFlux, xyra_DelRadFlux          & ! (inout)
!!$          & )
!!$
!!$      end if


      xyr_RadLUwFlux     = xyr_RadLUwFlux     + xyr_RadUwFlux
      xyr_RadLDwFlux     = xyr_RadLDwFlux     + xyr_RadDwFlux
      xyra_DelRadLUwFlux = xyra_DelRadLUwFlux + xyra_DelRadUwFlux
      xyra_DelRadLDwFlux = xyra_DelRadLDwFlux + xyra_DelRadDwFlux

    end do LOOP_band_RTE





!!$    i = 0
!!$    j = jmax / 2 + 1
!!$    write( 73, * ) xy_SurfTemp(i,j), 0.0d0, 0.0d0, xyr_Press(i,j,0)
!!$    do k = 1, kmax
!!$      write( 73, * ) xyz_Temp(i,j,k), xyz_QH2OVap(i,j,k), xyz_QO3(i,j,k), &
!!$        & xyz_Press(i,j,k)
!!$    end do
!!$    call flush( 73 )
!!$
!!$    i = 0
!!$    j = jmax / 2 + 1
!!$    do k = 1, kmax
!!$      write( 83, * ) &
!!$        & + (     xyr_RadLFlux(i,j,k-1) - xyr_RadLFlux(i,j,k) )  &
!!$        &     / ( xyr_Press(i,j,k-1)    - xyr_Press(i,j,k) )     &
!!$        &     / 1004.6 * 9.8, &
!!$        & xyz_Press(i,j,k)
!!$    end do
!!$    call flush( 83 )
!!$
!!$    i = 0
!!$    j = jmax / 2 + 1
!!$    do k = 0, kmax
!!$      write( 93, * ) xyr_RadLFlux(i,j,k), xyr_Press(i,j,k)
!!$    end do
!!$    call flush( 93 )
!!$    stop


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

  end subroutine RadEarthLWV24Flux
Subroutine :
FlagSnow :logical, intent(in)

This procedure input/output NAMELIST#rad_Earth_LW_V2_4_nml .

[Source]

  subroutine RadEarthLWV24Init( FlagSnow )

    ! USE statements
    !

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

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

    ! 暦と日時の取り扱い
    ! Calendar and Date handler
    !
    use dc_calendar, only: DCCalConvertByUnit

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

    ! Chou and Kouvaris (1991) による長波放射モデル
    ! Long radiation model described by Chou and Kouvaris (1991)
    !
    use rad_CK1991, only : RadCK1991Init

    ! Chou et al. (2001) による長波放射モデル
    ! Long radiation model described by Chou et al. (2001)
    !
    use rad_C2001, only : RadC2001Init

    ! 散乱を無視した放射伝達方程式
    ! Radiative transfer equation without considering scattering
    !
    use rad_rte_nonscat, only : RadRTENonScatInit

    ! 雲関系ルーチン
    ! Cloud-related routines
    !
    use cloud_utils, only : CloudUtilsInit


    logical, intent(in) :: FlagSnow


    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

    integer :: n


    namelist /rad_Earth_LW_V2_4_nml/ FlagHighAlt, DelTimeCalcTransValue, DelTimeCalcTransUnit, flag_save_time


    if ( rad_Earth_LW_V2_4_inited ) return


    FlagHighAlt           = .false.

    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 = rad_Earth_LW_V2_4_nml, iostat = iostat_nml )                  ! (out)
      close( unit_nml )

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

    ! Handle interval time
    !
    IntTimeSave = DCCalConvertByUnit( DelTimeCalcTransValue, DelTimeCalcTransUnit, 'sec' ) ! (in)




    do n = 1, nbmax
      ! unit conversion from (cm-1) to (m-1)
      aa_BandParam(1,n) = aa_BandParam(1,n) * 1.0d2
      aa_BandParam(2,n) = aa_BandParam(2,n) * 1.0d2
    end do


    ! allocate a variable for saving transmittance
    !
    allocate( xyrra_TransSaved (0:imax-1,1:jmax,0:kmax,0:kmax,1:nbmax) )
    allocate( xyrr_TransMASaved(0:imax-1,1:jmax,0:kmax,0:kmax)         )


    call RadEarthLWV24PrepPFTable


    ! Initialization of modules used in this module
    !

    ! Chou and Kouvaris (1991) による長波放射モデル
    ! Long radiation model described by Chou and Kouvaris (1991)
    !
    call RadCK1991Init

    ! Chou et al. (2001) による長波放射モデル
    ! Long radiation model described by Chou et al. (2001)
    !
    call RadC2001Init

    ! 散乱を無視した放射伝達方程式
    ! Radiative transfer equation without considering scattering
    !
    call RadRTENonScatInit

    ! 雲関系ルーチン
    ! Cloud-related routines
    !
    call CloudUtilsInit( FlagSnow )


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


    rad_Earth_LW_V2_4_inited = .true.

  end subroutine RadEarthLWV24Init

Private Instance methods

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
    !

    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
    !

    ! メッセージ出力
    ! 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


    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
DiffFactor
Constant :
DiffFactor = 1.66d0 :real(DP), parameter
FlagHighAlt
Variable :
FlagHighAlt :logical , save
FlagTransSaved
Variable :
FlagTransSaved :logical , save
IntTimeSave
Variable :
IntTimeSave :real(DP), save
: 長波フラックスを計算する時間間隔. Interval time of long wave flux calculation
PrevTimeSave
Variable :
PrevTimeSave :real(DP), save
: 前回長波フラックスを計算した時刻. Time when long wave flux is calculated
Subroutine :

[Source]

  subroutine RadEarthLWV24PrepPFTable

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

    ! ガウス重み, 分点の計算
    ! Calculate Gauss node and Gaussian weight
    !
    use gauss_quad, only : GauLeg

    ! プランク関数の計算
    ! Calculate Planck function
    !
    use planck_func, only : PF, DPFDT, Integ_PF_GQ_Array2D, Integ_DPFDT_GQ_Array2D

    integer , parameter :: NGaussQuad = 5
    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


    ! Preparation of tables for calculation of Plank function
    !
    TableTempMin       =  50.0d0
    TableTempMax       = 600.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( aa_BandParam(1,n), aa_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 Integ_PF_GQ_Array2D( aa_BandParam(1,n), aa_BandParam(2,n), NGaussQuad, 0, imax-1, 1, jmax, xy_TempTMP, xy_PF )
        call Integ_DPFDT_GQ_Array2D( aa_BandParam(1,n), aa_BandParam(2,n), NGaussQuad, 0, imax-1, 1, jmax, 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

  end subroutine RadEarthLWV24PrepPFTable
TableTempIncrement
Variable :
TableTempIncrement :real(DP), save
TableTempMax
Variable :
TableTempMax :real(DP), save
TableTempMin
Variable :
TableTempMin :real(DP), save
a_TableTemp
Variable :
a_TableTemp(:) :real(DP), save, allocatable
aa_BandParam
Variable :
aa_BandParam(1:2, 1:nbmax) :real(DP), save
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
module_name
Constant :
module_name = ‘rad_Earth_LW_V2_4 :character(*), parameter
: モジュールの名称. Module name
nbmax
Constant :
nbmax = 10 :integer , parameter
ntmax
Variable :
ntmax :integer , save
rad_Earth_LW_V2_4_inited
Variable :
rad_Earth_LW_V2_4_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag
version
Constant :
version = ’$Name: $’ // ’$Id: rad_Earth_LW_V2_4.f90,v 1.3 2015/01/29 12:06:43 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version
xyrr_TransMASaved
Variable :
xyrr_TransMASaved(:,:,:,:) :real(DP), allocatable, save
xyrra_TransSaved
Variable :
xyrra_TransSaved(:,:,:,:,:) :real(DP), allocatable, save