!= 短波入射 (太陽入射)
!
!= Short wave (insolation) incoming
!
! Authors::   Yoshiyuki O. TAKAHASHI, Satoshi NODA, Yasuhiro MORIKAWA, Yukiko YAMADA
! Version::   $Id: radiation_short_income.f90,v 1.13 2010-03-29 07:26:29 noda Exp $
! Tag Name::  $Name: dcpam5-20100424 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module radiation_short_income
  !
  != 短波入射 (太陽入射)
  !
  != Short wave (insolation) incoming
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 短波入射 (太陽入射) を計算します.
  !
  ! Calculate short wave (insolation) incoming radiation. 
  !
  !== Procedures List
  !
  ! ShortIncoming      :: 短波入射 (太陽入射) の計算
  ! ------------       :: ------------
  ! ShortIncoming      :: Calculate short wave (insolation) incoming radiation. 
  !
  !== NAMELIST
  !
  ! NAMELIST#radiation_short_income_nml
  !

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

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数. 
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数. 
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数. 
                             ! Number of vertical level

  ! 日付および時刻の取り扱い
  ! Date and time handler
  !
  use dc_date_types, only: DC_DIFFTIME
                              ! 日時の差を表現するデータ型. 
                              ! Data type for difference about date and time
  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING     ! 文字列.       Strings. 

  ! NAMELIST ファイル入力に関するユーティリティ
  ! Utilities for NAMELIST file input
  !
  use namelist_util, only: MaxNmlArySize
                              ! NAMELIST から読み込む配列の最大サイズ. 
                              ! Maximum size of arrays loaded from NAMELIST

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

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public:: ShortIncoming

  ! 公開変数
  ! Public variables
  !
  logical, save, public:: radiation_short_income_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag.


  ! 非公開変数
  ! Private variables
  !

  real(DP), save:: SolarConst
                              ! 太陽定数. 
                              ! Solar constant.
  integer,  save:: AnnualMean
                              ! 年平均入射フラグ.
                              ! Flag for annual mean incoming radiation. 
  integer,  save:: DiurnalMean
                              ! 日平均入射フラグ.
                              ! Flag for diurnal mean incoming radiation. 
  real(DP), save:: YearLengthScld
                              ! モデルの 1 年でスケーリングした 1 年の長さ.
                              ! Length of a year scaled by a year of model.
  real(DP), save:: DayLengthScld
                              ! モデルの 1 日でスケーリングした 1 日の長さ.
                              ! Length of a day scaled by a day of model.
  real(DP), save:: EpsOrb
                              ! 赤道傾斜角.
                              ! Inclination of equator to orbit.
  real(DP), save:: EqnOrb
                              ! 昇降点黄経. 
                              ! 
  real(DP), save:: Eccentricity
                              ! 離心率.
                              ! Eccentricity.
  real(DP), save:: MeanAnomary_init
                              ! 計算開始時の平均近点角 (0.0 - 2*PI). 
                              ! Mean anomary (0.0 - 2*PI).
  integer,  save:: MaxItrEccAnomary
                              ! 離心近点角を計算する時の最大繰り返し回数.
                              ! Maximum iteration number of times
                              ! to calculate eccentric anomary.
  real(DP), save:: ThreEccAnomaryError
                              ! 離心近点角を計算する時の誤差の許容しきい値.
                              ! Threshold of error
                              ! to calculate eccentric anomary.
  real(DP), save:: IncomAIns
                              ! $ A_{ins} $ . 年平均入射の係数. 
                              ! Coefficient of annual mean incoming radiation. 
  real(DP), save:: IncomBIns
                              ! $ B_{ins} $ . 年平均入射の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 
  real(DP), save:: IncomAZet
                              ! $ A_{\zeta} $ . 年平均入射角の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 
  real(DP), save:: IncomBZet
                              ! $ B_{\zeta} $ . 年平均入射角の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 


  real(DP), allocatable, save:: xy_IncomRadSFlux (:,:)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux
  real(DP), allocatable, save:: xy_InAngle (:,:)
                              ! sec (入射角). 
                              ! sec (angle of incidence)

  logical:: FlagRadiationSynchronous
                              ! 同期回転日射のフラグ
                              ! Flag for synchronous rotation

  character(*), parameter:: module_name = 'radiation_short_income'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20100424 $' // &
    & '$Id: radiation_short_income.f90,v 1.13 2010-03-29 07:26:29 noda Exp $'
                              ! モジュールのバージョン
                              ! Module version

  ! INTERFACE 文 ; INTERFACE statements
  !
  interface ShortIncoming
    module procedure ShortIncoming
  end interface

contains

  subroutine ShortIncoming( &
    & xy_IncomRadSFlux, xy_InAngle, DistFromStarScld, xy_CosZet  & ! (out)
    & )
    !
    ! 短波入射 (太陽入射) を計算します.
    !
    ! 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, & ! $ \varphi $ [rad.] . 緯度. Latitude
                     & 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 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, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, &                ! (in)
        & nml = radiation_short_income_nml, &  ! (out)
        & 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

  !-------------------------------------------------------------------

  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

end module radiation_short_income
