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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: PI       ! $ \pi $ .
                                  ! 円周率.  Circular constant
    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimeN

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

    use dc_date_types, only: DAY_SECONDS, YEAR_DAYS

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: y_Lat, x_Lon    ! $ \lambda $ [rad.] . 経度. Longitude

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

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(out):: xy_IncomRadSFlux (0:imax-1, 1:jmax)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux
    real(DP), intent(out):: xy_InAngle (0:imax-1, 1:jmax)
                              ! sec (入射角). 
                              ! sec (angle of incidence)
    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)
                              ! 入射角
                              ! Incidence angle

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

    real(DP):: HourAngle       ! 時角
                               ! Hour angle
    real(DP):: SeasonByYear    ! 季節を年で表記したもの (0.0 - 1.0)
                               ! Season expressed by year (0.0 - 1.0)
    real(DP):: ClockByDay      ! 時刻を日で表記したもの (0.0 - 1.0)
                               ! Clock expressed by day (0.0 - 1.0)
    real(DP):: MeanAnomary     ! 平均近点角
                               ! Mean anomary
    real(DP):: EccAnomary      ! 離心近点角
                               ! eccentric anomary
    real(DP):: EccAnomaryError ! ニュートン法における離心近点角の誤差
                               ! error of eccentric anomary in Newton method
    real(DP):: DistFromStarScldLR


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. radiation_short_income_inited ) call ShtIncomeInit

    ! 同期回転日射のフラグ
    ! Flag for synchronous rotation
    if ( .not. FlagRadiationSynchronous ) then
  
      ! 年, 日平均日射の計算
      ! Calculate annual mean, daily mean insolation
      !

      if ( (AnnualMean == 1) .and. (DiurnalMean == 1) ) then


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

            if ( xy_IncomRadSFlux(i,j) < 0.0_DP ) then
              xy_InAngle(i,j) = 1.0_DP / ( IncomAZet + IncomBZet * cos( y_Lat(j) )**2 )
            else
              xy_IncomRadSFlux(i,j) = 0.
              xy_InAngle(i,j) = 0.
            end if
          end do
        end do 
       
      else if ( (AnnualMean == 0) .and. (DiurnalMean == 0) ) then

        ! 季節変化, 日変化がある場合の計算
        ! Calculate with seasonal change and diurnal change
        !

        SeasonByYear = mod( EvalSec( TimeN ) / (YearLengthScld * YEAR_DAYS * DAY_SECONDS), 1.0_DP )
        ClockByDay   = mod( EvalSec( TimeN ) / (DayLengthScld * DAY_SECONDS), 1.0_DP )
        MeanAnomary  = mod( 2.0_DP * PI * SeasonByYear + MeanAnomary_init, 2.0_DP * PI )

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

        EccAnomary = MeanAnomary
        do itr = 1, MaxItrEccAnomary
           EccAnomaryError = EccAnomary - Eccentricity * sin(EccAnomary) - MeanAnomary
           if ( abs(EccAnomaryError) <= ThreEccAnomaryError ) exit
           EccAnomary = mod( EccAnomary - EccAnomaryError / ( 1.0_DP - Eccentricity * cos(EccAnomary) ), 2.0 * PI )
        end do
        if ( itr > MaxItrEccAnomary ) then
           call MessageNotify( 'E', module_name, 'Calculation for eccentric anomary does not converge.' )
        end if

        DistFromStarScldLR = 1.0_DP - Eccentricity * cos( EccAnomary )

        SinDel = sin( EpsOrb/180.0_DP * PI ) * sin( ( SeasonByYear - EqnOrb/360.0_DP ) * 2.0_DP * PI )

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

              HourAngle = ClockByDay * 2.0_DP * PI - PI + x_Lon(i)
              CosZet = sin( y_Lat(j) ) * SinDel + cos( y_Lat(j) ) * sqrt( 1.0_DP - SinDel**2 ) * cos( HourAngle )

              if ( CosZet > 0.0_DP ) then
                 xy_IncomRadSFlux(i,j) = - ( SolarConst / DistFromStarScldLR ** 2 ) * CosZet
                 xy_InAngle(i,j) = 1.0_DP / CosZet
              else
                 xy_IncomRadSFlux(i,j) = 0.
                 xy_InAngle(i,j) = 0.
              end if

              if ( present( xy_CosZet ) ) then
                 xy_CosZet(i,j) = CosZet
              end if

           end do
        end do


        if ( present( DistFromStarScld ) ) then
           DistFromStarScld = DistFromStarScldLR
        end if

      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 )
            if ( CosZet > 0.0_DP ) then
               xy_IncomRadSFlux( i,j ) = - SolarConst * CosZet
               xy_InAngle( i,j )       = 1.0_DP / CosZet
            else
               xy_IncomRadSFlux( i,j ) = 0.0_DP
               xy_InAngle( i,j )       = 0.0_DP
            end if
         end do
      end do

    end if

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'ISR', xy_IncomRadSFlux )


  end subroutine ShortIncoming

  subroutine InitCheck
    !
    ! 依存モジュールの初期化チェック
    !
    ! Check initialization of dependency modules

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

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: gridset_inited

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: constants_inited

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: axesset_inited

    ! 時刻管理
    ! Time control
    !
    use timeset, only: timeset_inited

    ! 実行文 ; Executable statement
    !

    if ( .not. namelist_util_inited ) call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )

    if ( .not. gridset_inited ) call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )

    if ( .not. constants_inited ) call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )

    if ( .not. axesset_inited ) call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )

    if ( .not. timeset_inited ) call MessageNotify( 'E', module_name, '"timeset" module is not initialized.' )

  end subroutine InitCheck

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

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

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

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

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    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_short_income_nml/ FlagRadiationSynchronous, SolarConst, AnnualMean, DiurnalMean, YearLengthScld, DayLengthScld, EpsOrb, EqnOrb, Eccentricity, MeanAnomary_init, MaxItrEccAnomary, ThreEccAnomaryError, IncomAIns, IncomBIns, IncomAZet, IncomBZet
          !
          ! デフォルト値については初期化手続 "radiation_short_income#ShtIncomeInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_short_income#ShtIncomeInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( radiation_short_income_inited ) return
    call InitCheck

    ! デフォルト値の設定
    ! Default values settings
    !
    FlagRadiationSynchronous = .false.
    SolarConst           = 1380.0_DP
    AnnualMean           = 1
    DiurnalMean          = 1
    YearLengthScld       = 1.0_DP
    DayLengthScld        = 1.0_DP
    EpsOrb               = 23.5_DP
    EqnOrb               = 110.0_DP
    Eccentricity         = 0.0_DP
    MeanAnomary_init     = 0.0_DP
    MaxItrEccAnomary     = 20
    ThreEccAnomaryError  = 1e-6_DP
    IncomAIns            = 0.127_DP
    IncomBIns            = 0.183_DP
    IncomAZet            = 0.410_DP
    IncomBZet            = 0.590_DP

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

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

    ! 保存用の変数の割り付け
    ! 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' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'ShortIncomming:' )
    call MessageNotify( 'M', module_name, '  FlagRadiationSynchronous = <%b>', l = (/ FlagRadiationSynchronous /) )
    call MessageNotify( 'M', module_name, '  SolarConst               = %f', d = (/ SolarConst               /) )
    call MessageNotify( 'M', module_name, '  AnnualMean               = %d', i = (/ AnnualMean               /) )
    call MessageNotify( 'M', module_name, '  DiurnalMean              = %d', i = (/ DiurnalMean              /) )
    call MessageNotify( 'M', module_name, '  YearLengthScld           = %f', d = (/ YearLengthScld           /) )
    call MessageNotify( 'M', module_name, '  DayLengthScld            = %f', d = (/ DayLengthScld            /) )
    call MessageNotify( 'M', module_name, '  EpsOrb                   = %f', d = (/ EpsOrb                   /) )
    call MessageNotify( 'M', module_name, '  EqnOrb                   = %f', d = (/ EqnOrb                   /) )
    call MessageNotify( 'M', module_name, '  Eccentricity             = %f', d = (/ Eccentricity             /) )
    call MessageNotify( 'M', module_name, '  MeanAnomary_init         = %f', d = (/ MeanAnomary_init         /) )
    call MessageNotify( 'M', module_name, '  MaxItrEccAnomary         = %d', i = (/ MaxItrEccAnomary         /) )
    call MessageNotify( 'M', module_name, '  ThreEccAnomaryError      = %f', d = (/ ThreEccAnomaryError      /) )
    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, '-- version = %c', c1 = trim(version) )

    radiation_short_income_inited = .true.
  end subroutine ShtIncomeInit