| Class | rad_short_income |
| In: |
radiation/rad_short_income.f90
|
Note that Japanese and English are described in parallel.
短波入射 (太陽入射) を計算します.
Calculate short wave (insolation) incoming radiation.
| ShortIncoming : | 短波入射 (太陽入射) の計算 |
| ———— : | ———— |
| ShortIncoming : | Calculate short wave (insolation) incoming radiation. |
| Subroutine : | |||
| xy_InAngle(0:imax-1, 1:jmax) : | real(DP), intent(out), optional
| ||
| DistFromStarScld : | real(DP), intent(out), optional
| ||
| xy_CosZet(0:imax-1, 1:jmax) : | real(DP), intent(out), optional
| ||
| DiurnalMeanFactor : | real(DP), intent(out), optional | ||
| PlanetLonFromVE : | real(DP), intent(out), optional
| ||
| FlagOutput : | logical , intent(in ), optional |
短波入射 (太陽入射) を計算します.
Calculate short wave (insolation) incoming radiation.
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 .
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
| Variable : | |||
| FlagAnnualMean : | logical, save
|
| Variable : | |||
| FlagDiurnalMean : | logical, save
|
| Variable : | |||
| FlagPerpetual : | logical, save
|
| Variable : | |||
| FlagRadSynchronous : | logical
|
| Variable : | |||
| IncomAIns : | real(DP), save
|
| Variable : | |||
| IncomAZet : | real(DP), save
|
| Variable : | |||
| IncomBIns : | real(DP), save
|
| Variable : | |||
| IncomBZet : | real(DP), save
|
| Variable : | |||
| LonFromVEAtEpoch : | real(DP), save
|
| Variable : | |||
| MaxItrEccAnomaly : | integer, save
|
| Variable : | |||
| PerLonFromVE : | real(DP), save
|
| Variable : | |||
| PerpDistFromStarScld : | real(DP), save
|
| Variable : | |||
| PerpSinDel : | real(DP), save
|
| Subroutine : | |||
| SinDel : | real(DP), intent(out)
| ||
| DistFromStarScld : | real(DP), intent(out)
| ||
| PlanetLonFromVE : | real(DP), intent(out)
|
短波入射 (太陽入射) を計算します.
Calculate short wave (insolation) incoming radiation.
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
| Variable : | |||
| ThreEccAnomalyError : | real(DP), save
|
| Constant : | |||
| module_name = ‘rad_short_income‘ : | character(*), parameter
|
| Variable : | |||
| rad_short_income_inited = .false. : | logical, save
|
| Constant : | |||
| version = ’$Name: $’ // ’$Id: rad_short_income.f90,v 1.6 2013/05/25 06:33:57 yot Exp $’ : | character(*), parameter
|
| Variable : | |||
| xy_InAngle(:,:) : | real(DP), allocatable, save
|