Class rad_short_income
In: radiation/rad_short_income.f90

短波入射 (太陽入射)

Short wave (insolation) incoming

Note that Japanese and English are described in parallel.

短波入射 (太陽入射) を計算します.

Calculate short wave (insolation) incoming radiation.

Procedures List

ShortIncoming :短波入射 (太陽入射) の計算
———— :————
ShortIncoming :Calculate short wave (insolation) incoming radiation.

NAMELIST

NAMELIST#rad_short_income_nml

Methods

Included Modules

dc_types dc_message constants0 gridset timeset axesset dc_calendar gtool_historyauto dc_iounit namelist_util

Public Instance methods

Subroutine :
xy_InAngle(0:imax-1, 1:jmax) :real(DP), intent(out), optional
: sec ( 入射角 ). sec ( Incident angle )
DistFromStarScld :real(DP), intent(out), optional
: 軌道長半径でスケーリングした恒星からの距離 Distance from the star scaled by semimajor axis of the planet‘s orbit
xy_CosZet(0:imax-1, 1:jmax) :real(DP), intent(out), optional
: cos( 入射角 ) cos( Incident angle )
DiurnalMeanFactor :real(DP), intent(out), optional
PlanetLonFromVE :real(DP), intent(out), optional
: 中心星を中心とする座標における春分点から測った惑星の経度 Longitude of the planet measured from the vernal equinox in the coordinate that the central star is located on the origin.
FlagOutput :logical , intent(in ), optional

短波入射 (太陽入射) を計算します.

Calculate short wave (insolation) incoming radiation.

[Source]

  subroutine RadShortIncome( xy_InAngle, DistFromStarScld, xy_CosZet, DiurnalMeanFactor, PlanetLonFromVE, FlagOutput )
    !
    ! 短波入射 (太陽入射) を計算します.
    !
    ! Calculate short wave (insolation) incoming radiation. 
    !

    ! モジュール引用 ; USE statements
    !

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

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


    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(out), optional :: xy_InAngle (0:imax-1, 1:jmax)
                              ! sec ( 入射角 ). 
                              ! sec ( Incident angle )
    real(DP), intent(out), optional :: DistFromStarScld
                              ! 軌道長半径でスケーリングした恒星からの距離
                              ! Distance from the star scaled
                              ! by semimajor axis of the planet's orbit
    real(DP), intent(out), optional :: xy_CosZet(0:imax-1, 1:jmax)
                              ! cos( 入射角 )
                              ! cos( Incident angle )
    real(DP), intent(out), optional :: DiurnalMeanFactor
    real(DP), intent(out), optional :: PlanetLonFromVE
                              ! 中心星を中心とする座標における春分点から測った惑星の経度
                              ! Longitude of the planet measured from the vernal equinox
                              ! in the coordinate that the central star is located on 
                              ! the origin.
    logical , intent(in ), optional :: FlagOutput


    ! 作業変数
    ! Work variables
    !
    integer:: i                ! 経度方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in longitude
    integer:: j                ! 緯度方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in latitude
    real(DP):: SinDel          ! 赤緯
                               ! Declination
    real(DP):: CosZet          ! 入射角
                               ! Incidence angle
    real(DP):: AngleMaxLon     ! 入射が最大となる緯度

    real(DP):: HourAngle       ! 時角
                               ! Hour angle
    real(DP):: ClockByDay      ! 時刻を日で表記したもの (0.0 - 1.0)
                               ! Clock expressed by day (0.0 - 1.0)


    real(DP) :: xy_InAngleLV       (0:imax-1, 1:jmax)
                               ! sec ( 入射角 ). 
                               ! sec ( Incident angle )
                               ! (local variable)
    real(DP) :: DistFromStarScldLV
                               ! Distance between the central star and the planet
                               ! (local variable)
    real(DP) :: xy_CosZetLV        (0:imax-1, 1:jmax)
                               ! cos( 入射角 )
                               ! cos( Incident angle )
                               ! (local variable)
    real(DP) :: PlanetLonFromVELV

    integer         :: hour_in_a_day, min_in_a_hour
    real(DP)        :: sec_in_a_min, sec_in_a_day


    ! 実行文 ; Executable statement
    !

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


    ! Set flag for diurnally averaged insolation : This is temporary one.
    !
    if ( present( DiurnalMeanFactor ) ) then
      if ( ( .not. FlagAnnualMean ) .and. FlagDiurnalMean ) then
!!$        DiurnalMeanFactor = 1.0_DP / PI
        DiurnalMeanFactor = -1.0e100_DP
      else
        DiurnalMeanFactor = 1.0_DP
      end if
    end if


    ! 同期回転日射のフラグ
    ! Flag for synchronous rotation
    if ( .not. FlagRadSynchronous ) then

      ! 年, 日平均日射の計算
      ! Calculate annual mean, daily mean insolation
      !

      if ( FlagAnnualMean .and. FlagDiurnalMean ) then

        do i = 0, imax - 1
          do j = 1, jmax
            xy_CosZetLV(i,j) = IncomAIns + IncomBIns * cos( y_Lat(j) )**2

            if ( xy_CosZetLV(i,j) > 0.0_DP ) then
              xy_InAngleLV(i,j) = 1.0_DP / xy_CosZetLV(i,j)
            else
              xy_InAngleLV(i,j) = 0.
            end if

          end do
        end do

        DistFromStarScldLV = 1.0_DP
        SinDel             = 0.0_DP
        PlanetLonFromVELV  = 0.0_DP

      else if ( .not. FlagAnnualMean ) then

        ! Set sine of declination and distance between a planet and a central star 
        ! scaled with semi-major axis
        !
        if ( FlagPerpetual ) then
          SinDel             = PerpSinDel
          DistFromStarScldLV = PerpDistFromStarScld
        else
          call ShortIncomCalcOrbParam( SinDel, DistFromStarScldLV, PlanetLonFromVELV )
          if ( present( PlanetLonFromVE ) ) PlanetLonFromVE = PlanetLonFromVELV
        end if


        call DCCalInquire( hour_in_day      = hour_in_a_day, min_in_hour      = min_in_a_hour, sec_in_min       = sec_in_a_min )
        sec_in_a_day  = hour_in_a_day * min_in_a_hour * sec_in_a_min


        if ( FlagSpecifySolDay ) then
          ! case with solar day which is different from sec_in_a_day (rotation
          ! period)
          ClockByDay = mod( TimeN / SolDay, 1.0_DP )
        else
          ! case with solar day which is the same as sec_in_a_day (rotation
          ! period)
          ClockByDay = DCCalDateEvalSecOfDay( TimeN, date = InitialDate )
          ClockByDay = ClockByDay / sec_in_a_day
        end if


!!$        write( 6, * ) '###', TimeN, ClockByDay, mod( TimeN / sec_in_a_day, 1.0_DP )
!!$        write( 60, * ) TimeN, ClockByDay, mod( TimeN / sec_in_a_day, 1.0_DP )
!!$        call flush( 60 )


        do i = 0, imax - 1
          do j = 1, jmax


            if ( FlagDiurnalMean ) then
              HourAngle = 0.0_DP
            else
              HourAngle = ClockByDay * 2.0_DP * PI - PI + x_Lon(i)
            end if

            CosZet = sin( y_Lat(j) ) * SinDel + cos( y_Lat(j) ) * sqrt( 1.0_DP - SinDel**2 ) * cos( HourAngle )

            xy_CosZetLV(i,j) = CosZet

            if ( CosZet > 0.0_DP ) then
              xy_InAngleLV(i,j)       = 1.0_DP / CosZet
            else
              xy_InAngleLV(i,j)       = 0.
            end if

          end do
        end do

      else

        ! 対応していない入射タイプ
        ! not implemented insolation type
        !
        call MessageNotify( 'E', module_name, 'This type of insolation is not impremented' )

      end if

    else

      ! 短波入射 (太陽入射) を計算します.
      ! 同期回転惑星を想定しており,
      ! 常に経度 90 度で最大入射, 経度 180-360 度で入射ゼロとなっています.
      !
      ! Calculate short wave (insolation) incoming radiation.
      ! A planet with synchronized rotation are assumed.
      ! Incoming is max at latitude 90 deg., and zero between 180 and 360 deg. constantly.
      do i = 0, imax - 1
         AngleMaxLon  = - PI / 2.0_DP + x_Lon( i )
         do j = 1, jmax
            CosZet = cos( y_Lat(j) ) * cos( AngleMaxLon )

            xy_CosZetLV(i,j) = CosZet

            if ( CosZet > 0.0_DP ) then
               xy_InAngleLV(i,j) = 1.0_DP / CosZet
            else
               xy_InAngleLV(i,j) = 0.0_DP
            end if

            DistFromStarScldLV = 1.0_DP
         end do
      end do

    end if


    if ( present( xy_InAngle ) ) then
      xy_InAngle = xy_InAngleLV
    end if
    if ( present( xy_CosZet ) ) then
      xy_CosZet = xy_CosZetLV
    end if
    if ( present( DistFromStarScld ) ) then
      DistFromStarScld = DistFromStarScldLV
    end if

    ! ヒストリデータ出力
    ! History data output
    !
    if ( present( FlagOutput ) ) then
      if ( FlagOutput ) then
        call HistoryAutoPut( TimeN, 'Decl'            , asin(SinDel)*180.0_DP/PI      )
        call HistoryAutoPut( TimeN, 'DistFromStarScld', DistFromStarScldLV            )
        call HistoryAutoPut( TimeN, 'PlanetLonFromVE' , PlanetLonFromVELV*180.0_DP/PI )
      end if
    else
      call HistoryAutoPut( TimeN, 'Decl'            , asin(SinDel)*180.0_DP/PI      )
      call HistoryAutoPut( TimeN, 'DistFromStarScld', DistFromStarScldLV            )
      call HistoryAutoPut( TimeN, 'PlanetLonFromVE' , PlanetLonFromVELV*180.0_DP/PI )
    end if


  end subroutine RadShortIncome
Subroutine :

rad_short_income モジュールの初期化を行います. NAMELIST#rad_short_income_nml の読み込みはこの手続きで行われます.

"rad_short_income" module is initialized. "NAMELIST#rad_short_income_nml" is loaded in this procedure.

This procedure input/output NAMELIST#rad_short_income_nml .

[Source]

  subroutine RadShortIncomeInit
    !
    ! rad_short_income モジュールの初期化を行います. 
    ! NAMELIST#rad_short_income_nml の読み込みはこの手続きで行われます. 
    !
    ! "rad_short_income" module is initialized. 
    ! "NAMELIST#rad_short_income_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

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

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

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    integer:: EpochYear, EpochMonth, EpochDay, EpochHour, EpochMin
                              ! 元期日時 (年月日時分).
                              ! "TimeAtEpoch" が負の場合にこちらが使用される.
                              !
                              ! Date at epoch (year, month, day, hour, minute)
                              ! These are used when "TimeAtEpoch" is negative.
    real(DP):: EpochSec
                              ! 元期日時 (秒).
                              ! "TimeAtEpoch" が負の場合にこちらが使用される.
                              !
                              ! Date at epoch (second)
                              ! These are used when "TimeAtEpoch" is negative.

    type(DC_CAL_DATE):: EpochDate
                              ! 元期の日時
                              ! Date at epoch

    real(DP)         :: PerpDelDeg
                              ! Declination angle in degree used for perpetual experiment

    real(DP)         :: SolDayValue
    character(TOKEN) :: SolDayUnit

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

    logical         :: FlagUseOfEpochDate
    character(TOKEN):: date_print


    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /rad_short_income_nml/ FlagRadSynchronous, FlagAnnualMean, FlagDiurnalMean, FlagPerpetual, PerpDelDeg, PerpDistFromStarScld, EpsOrb, PerLonFromVE, LonFromVEAtEpoch, Eccentricity, TimeAtEpoch, EpochYear, EpochMonth, EpochDay, EpochHour, EpochMin, EpochSec, MaxItrEccAnomaly, ThreEccAnomalyError, IncomAIns, IncomBIns, IncomAZet, IncomBZet, FlagSpecifySolDay, SolDayValue, SolDayUnit
          !
          ! デフォルト値については初期化手続 "rad_short_income#RadShortIncomeInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_short_income#RadShortIncomeInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( rad_short_income_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    FlagRadSynchronous       = .false.
    FlagAnnualMean           = .true.
    FlagDiurnalMean          = .true.
    FlagPerpetual            = .false.

    !---

    PerpDelDeg               = 0.0_DP
    PerpDistFromStarScld     = 1.0_DP

    !---

    EpsOrb               =   23.5_DP    ! Earth-like value
    PerLonFromVE         =    0.0_DP
    LonFromVEAtEpoch     =  280.0_DP    ! This value results in the fact that the planet
                                        ! is located at the position of vernal equinox 
                                        ! at 80 days after calculation with the use of 
                                        ! "360day" calendar.
    Eccentricity         =    0.0_DP
    TimeAtEpoch          =    0.0_DP
    EpochYear            =   -1
    EpochMonth           =   -1
    EpochDay             =   -1
    EpochHour            =   -1
    EpochMin             =   -1
    EpochSec             =   -1.0_DP

    !---

    ! Sample values for the Earth
    !  References: 
    !    Duffett-Smith, P., Practical astronomy with your calculator Third Edition, 
    !    Cambridge University Press, pp.185, 1988.
    !
!!$    EpsOrb               =   23.44_DP                ! Rika nenpyo (Chronological 
!!$                                                     ! Scientific Tables 2010)
!!$    PerLonFromVE         =  102.768413_DP + 180.0_DP ! Duffett-Smith (1988), p.105
!!$                                                     ! modified (plus 180 degrees)
!!$    LonFromVEAtEpoch     =   99.403308_DP + 180.0_DP ! Duffett-Smith (1988), p.105
!!$                                                     ! modified (plus 180 degrees)
!!$    Eccentricity         =    0.016713_DP            ! Duffett-Smith (1988), p.105
!!$    TimeAtEpoch          =   -1.0_DP                 ! EpochXXX written below are used 
!!$                                                     ! because this is negative.
!!$    EpochYear            = 1990                      ! Duffett-Smith (1988), p.105
!!$    EpochMonth           =    1
!!$    EpochDay             =    1
!!$    EpochHour            =    0
!!$    EpochMin             =    0
!!$    EpochSec             =    0.0_DP
    !---

    ! Sample values for Mars
    !  References: 
    !    Allison, M., Geophys. Res. Lett., 24, 1967-1970, 1997.
    !
!!$    EpsOrb               =   25.19_DP              ! Allison (1997)
!!$    PerLonFromVE         =  250.98_DP              ! Allison (1997) (modified)
!!$    LonFromVEAtEpoch     =  -10.342_DP             ! Arbitrarily set for clarity
!!$                                                   ! This results in Ls ~ 0 at Time = 0
!!$    Eccentricity         =    0.0934_DP            ! Allison (1997), value at epoch J2000
!!$    TimeAtEpoch          =    0.0_DP               ! Arbitrarily set for clarity
!!$    EpochYear            =   -1                    ! not used
!!$    EpochMonth           =   -1
!!$    EpochDay             =   -1
!!$    EpochHour            =   -1
!!$    EpochMin             =   -1
!!$    EpochSec             =   -1.0_DP

    !---

    MaxItrEccAnomaly     = 20
    ThreEccAnomalyError  = 1e-6_DP

    IncomAIns            = 0.127_DP   ! see dcpam document for reference
    IncomBIns            = 0.183_DP   ! see dcpam document for reference
    IncomAZet            = 0.410_DP   ! see dcpam document for reference
    IncomBZet            = 0.590_DP   ! see dcpam document for reference

    FlagSpecifySolDay = .false.
    SolDayValue       = 0.0_DP
    SolDayUnit        = 'sec'


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

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


    if ( ( .not. FlagAnnualMean ) .and. FlagDiurnalMean ) then
      call MessageNotify( 'E', module_name, 'FlagDiurnalMean cannot be true if FlagAnnualMean is false.' )
    end if


    PerpSinDel = sin( PerpDelDeg * PI / 180.0_DP )


    FlagUseOfEpochDate = .false.

    if ( TimeAtEpoch < 0.0_DP ) then
      call DCCalDateCreate( year  = EpochYear, month = EpochMonth, day   = EpochDay, hour  = EpochHour, min   = EpochMin, sec   = EpochSec, date  = EpochDate )    ! (out) optional

      TimeAtEpoch = DCCalDateDifference( start_date = InitialDate, end_date   = EpochDate )

      FlagUseOfEpochDate = .true.
    end if

    SolDay = DCCalConvertByUnit( SolDayValue, SolDayUnit, 'sec' )

    ! 保存用の変数の割り付け
    ! Allocate variables for saving
    !

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
!!$    call HistoryAutoAddVariable( 'xxxxx' , &
!!$      & (/ 'lon ', 'lat ', 'sig', 'time'/), &
!!$      & 'xxxx', 'W m-2' )

    call HistoryAutoAddVariable( 'ISR' , (/ 'lon ', 'lat ', 'time'/), 'incoming shortwave', 'W m-2' )

    call HistoryAutoAddVariable( 'Decl' , (/ 'time'/), 'declination angle', 'degree' )

    call HistoryAutoAddVariable( 'DistFromStarScld' , (/ 'time'/), 'distance between the central star and the planet', '1' )

    call HistoryAutoAddVariable( 'PlanetLonFromVE' , (/ 'time'/), 'planetary longitude from the vernal equinox', 'degree' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'ShortIncomming:' )
    call MessageNotify( 'M', module_name, '  FlagRadSynchronous       = %b', l = (/ FlagRadSynchronous /) )
    call MessageNotify( 'M', module_name, '  FlagAnnualMean           = %b', l = (/ FlagAnnualMean            /) )
    call MessageNotify( 'M', module_name, '  FlagDiurnalMean          = %b', l = (/ FlagDiurnalMean           /) )
    call MessageNotify( 'M', module_name, '  FlagPerpetual            = %b', l = (/ FlagPerpetual             /) )
    call MessageNotify( 'M', module_name, '  PerpDelDeg               = %f', d = (/ PerpDelDeg                /) )
    call MessageNotify( 'M', module_name, '  PerpDistFromStarScld     = %f', d = (/ PerpDistFromStarScld      /) )
    call MessageNotify( 'M', module_name, '  EpsOrb                   = %f', d = (/ EpsOrb                   /) )
    call MessageNotify( 'M', module_name, '  PerLonFromVE             = %f', d = (/ PerLonFromVE             /) )
    call MessageNotify( 'M', module_name, '  Eccentricity             = %f', d = (/ Eccentricity             /) )

    if ( FlagUseOfEpochDate ) then
      call DCCalDateInquire( date_print, date = EpochDate )
      call MessageNotify( 'M', module_name, '  EpochDate  = %c', c1 = trim(date_print) )
    end if
    call MessageNotify( 'M', module_name, '  TimeAtEpoch              = %f', d = (/ TimeAtEpoch              /) )
    call MessageNotify( 'M', module_name, '  LonFromVEAtEpoch         = %f', d = (/ LonFromVEAtEpoch         /) )

    call MessageNotify( 'M', module_name, '  MaxItrEccAnomaly         = %d', i = (/ MaxItrEccAnomaly         /) )
    call MessageNotify( 'M', module_name, '  ThreEccAnomalyError      = %f', d = (/ ThreEccAnomalyError      /) )
    call MessageNotify( 'M', module_name, '  IncomAIns                = %f', d = (/ IncomAIns                /) )
    call MessageNotify( 'M', module_name, '  IncomBIns                = %f', d = (/ IncomBIns                /) )
    call MessageNotify( 'M', module_name, '  IncomAZet                = %f', d = (/ IncomAZet                /) )
    call MessageNotify( 'M', module_name, '  IncomBZet                = %f', d = (/ IncomBZet                /) )

    call MessageNotify( 'M', module_name, '  FlagSpecifySolDay        = %b', l = (/ FlagSpecifySolDay        /) )
    call MessageNotify( 'M', module_name, '  SolDay                   = %f', d = (/ SolDay                   /) )


    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    rad_short_income_inited = .true.

  end subroutine RadShortIncomeInit

Private Instance methods

Eccentricity
Variable :
Eccentricity :real(DP), save
: 離心率. Eccentricity.
EpsOrb
Variable :
EpsOrb :real(DP), save
: 赤道傾斜角. Inclination of equator to orbit.
FlagAnnualMean
Variable :
FlagAnnualMean :logical, save
: 年平均入射フラグ. Flag for annual mean incoming radiation.
FlagDiurnalMean
Variable :
FlagDiurnalMean :logical, save
: 日平均入射フラグ. Flag for diurnal mean incoming radiation.
FlagPerpetual
Variable :
FlagPerpetual :logical, save
: 季節変化なし入射フラグ. Flag for perpetual incoming radiation.
FlagRadSynchronous
Variable :
FlagRadSynchronous :logical
: 同期回転日射のフラグ Flag for synchronous rotation
FlagSpecifySolDay
Variable :
FlagSpecifySolDay :logical, save
: Flag for specifying solar day
IncomAIns
Variable :
IncomAIns :real(DP), save
: $ A_{ins} $ . 年平均入射の係数. Coefficient of annual mean incoming radiation.
IncomAZet
Variable :
IncomAZet :real(DP), save
: $ A_{zeta} $ . 年平均入射角の係数. AIns に同じ. Coefficient of annual mean incoming radiation. Same as "AIns".
IncomBIns
Variable :
IncomBIns :real(DP), save
: $ B_{ins} $ . 年平均入射の係数. AIns に同じ. Coefficient of annual mean incoming radiation. Same as "AIns".
IncomBZet
Variable :
IncomBZet :real(DP), save
: $ B_{zeta} $ . 年平均入射角の係数. AIns に同じ. Coefficient of annual mean incoming radiation. Same as "AIns".
LonFromVEAtEpoch
Variable :
LonFromVEAtEpoch :real(DP), save
: 元期における惑星の経度 (黄経) (degree) Longitude of the planet at epoch (degree)
MaxItrEccAnomaly
Variable :
MaxItrEccAnomaly :integer, save
: 離心近点角を計算する時の最大繰り返し回数. Maximum iteration number of times to calculate eccentric anomaly.
PerLonFromVE
Variable :
PerLonFromVE :real(DP), save
: 春分から測った近日点の経度 (近日点黄経) (degree) Longitude of the perihelion from vernal equinox (degree)
PerpDistFromStarScld
Variable :
PerpDistFromStarScld :real(DP), save
: distance between a planet and a central star for perpetual experiment
PerpSinDel
Variable :
PerpSinDel :real(DP), save
: sine of declination angle for perpetual experiments
Subroutine :
SinDel :real(DP), intent(out)
: 赤緯 Declination
DistFromStarScld :real(DP), intent(out)
: 軌道長半径でスケーリングした恒星からの距離 Distance from the star scaled by semimajor axis of the planet‘s orbit
PlanetLonFromVE :real(DP), intent(out)
: 中心星を中心とする座標における春分点から測った惑星の経度 Longitude of the planet measured from the vernal equinox in the coordinate that the central star is located on the origin.

短波入射 (太陽入射) を計算します.

Calculate short wave (insolation) incoming radiation.

[Source]

  subroutine ShortIncomCalcOrbParam( SinDel, DistFromStarScld, PlanetLonFromVE )
    !
    ! 短波入射 (太陽入射) を計算します.
    !
    ! Calculate short wave (insolation) incoming radiation. 
    !

    ! モジュール引用 ; USE statements
    !

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

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(out) :: SinDel
                              ! 赤緯
                              ! Declination
    real(DP), intent(out) :: DistFromStarScld
                              ! 軌道長半径でスケーリングした恒星からの距離
                              ! Distance from the star scaled
                              ! by semimajor axis of the planet's orbit
    real(DP), intent(out) :: PlanetLonFromVE
                              ! 中心星を中心とする座標における春分点から測った惑星の経度
                              ! Longitude of the planet measured from the vernal equinox
                              ! in the coordinate that the central star is located on 
                              ! the origin.

    ! 作業変数
    ! Work variables
    !
    integer:: itr              ! イテレーション方向に回る DO ループ用作業変数
                               ! Work variables for DO loop in iteration direction

    real(DP):: MeanAnomaly     ! 平均近点角
                               ! Mean anomaly
    real(DP):: EccAnomaly      ! 離心近点角
                               ! eccentric anomaly
    real(DP):: EccAnomalyError ! ニュートン法における離心近点角の誤差
                               ! error of eccentric anomaly in Newton method
    real(DP):: TrueAnomaly     ! 真点離角
                               ! true anomaly

    integer         :: hour_in_a_day, min_in_a_hour, day_in_a_year
    integer, pointer:: day_in_month_ptr(:) => null()
    real(DP)        :: sec_in_a_min, sec_in_a_day, sec_in_a_year



    ! 実行文 ; Executable statement
    !

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


    ! 季節変化, 日変化がある場合の計算
    ! Calculate with seasonal change and diurnal change
    !
    call DCCalInquire( day_in_month_ptr = day_in_month_ptr , hour_in_day      = hour_in_a_day, min_in_hour      = min_in_a_hour, sec_in_min       = sec_in_a_min )
    ! Add 1 to day_in_month_ptr(2) if it is leap year.
    !
    if ( DCCalDateChkLeapYear( TimeN, InitialDate ) ) then
      day_in_month_ptr(2) = day_in_month_ptr(2) + 1
    end if

    day_in_a_year = sum( day_in_month_ptr )
    deallocate( day_in_month_ptr )
    sec_in_a_day  = hour_in_a_day * min_in_a_hour * sec_in_a_min
    sec_in_a_year = day_in_a_year * sec_in_a_day

    MeanAnomaly  = 2.0_DP * PI * ( TimeN - TimeAtEpoch ) / sec_in_a_year + ( LonFromVEAtEpoch - PerLonFromVE ) * PI / 180.0_DP
    MeanAnomaly  = mod( MeanAnomaly, 2.0_DP * PI )


    ! ニュートン法を用いて平均近点角から離心近点角を求める.
    ! Calculate eccentric anomaly from mean anomaly by using Newton method.

    EccAnomaly = MeanAnomaly
    do itr = 1, MaxItrEccAnomaly
      EccAnomalyError = EccAnomaly - Eccentricity * sin(EccAnomaly) - MeanAnomaly
      if ( abs(EccAnomalyError) <= ThreEccAnomalyError ) exit
      EccAnomaly      = EccAnomaly - EccAnomalyError / ( 1.0_DP - Eccentricity * cos(EccAnomaly) )
      EccAnomaly      = mod( EccAnomaly, 2.0 * PI )
    end do
    if ( itr > MaxItrEccAnomaly ) then
      call MessageNotify( 'E', module_name, 'Calculation for eccentric anomaly does not converge.' )
    end if

    DistFromStarScld = 1.0_DP - Eccentricity * cos( EccAnomaly )

    TrueAnomaly = 2.0_DP * atan( sqrt( ( 1.0d0 + Eccentricity ) / ( 1.0d0 - Eccentricity ) ) * tan( EccAnomaly / 2.0_DP ) )

    PlanetLonFromVE = TrueAnomaly + PerLonFromVE * PI / 180.0_DP
    PlanetLonFromVE = mod( PlanetLonFromVE, 2.0_DP * PI )

    SinDel = sin( EpsOrb * PI / 180.0_DP ) * sin( PlanetLonFromVE )


        ! code for debug
!!$        write( 60, * ) TimeN/sec_in_a_day, DCCalDateChkLeapYear(TimeN,date=InitialDate), day_in_a_year
!!$        write(  6, * ) TimeN/sec_in_a_day, DCCalDateChkLeapYear(TimeN,date=InitialDate), day_in_a_year
!!$        call flush( 60 )


!!$        write( 60, * ) TimeN/sec_in_a_day, asin(SinDel)*180.0/PI, DistFromStarScld, PlanetLonFromVE*180.0_DP/PI
!!$        write(  6, * ) TimeN/sec_in_a_day, asin(SinDel)*180.0/PI, DistFromStarScld, PlanetLonFromVE*180.0_DP/PI
!!$        call flush( 60 )


  end subroutine ShortIncomCalcOrbParam
SolDay
Variable :
SolDay :real(DP), save
: length of solar day
ThreEccAnomalyError
Variable :
ThreEccAnomalyError :real(DP), save
: 離心近点角を計算する時の誤差の許容しきい値. Threshold of error to calculate eccentric anomaly.
TimeAtEpoch
Variable :
TimeAtEpoch :real(DP), save
: 元期における時刻 (sec) Time at epoch (sec)
module_name
Constant :
module_name = ‘rad_short_income :character(*), parameter
: モジュールの名称. Module name
rad_short_income_inited
Variable :
rad_short_income_inited = .false. :logical, save
: 初期設定フラグ. Initialization flag.
version
Constant :
version = ’$Name: $’ // ’$Id: rad_short_income.f90,v 1.6 2013/05/25 06:33:57 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version
xy_InAngle
Variable :
xy_InAngle(:,:) :real(DP), allocatable, save
: sec (入射角). sec (angle of incidence)