scale_atmos_sub_solarins.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_atmos_solarins
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_solarins_setup
  public :: atmos_solarins_orbit
  public :: atmos_solarins_insolation

  interface atmos_solarins_insolation
     module procedure atmos_solarins_insolation_0d
     module procedure atmos_solarins_insolation_2d
  end interface atmos_solarins_insolation

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public :: atmos_solarins_constant    = 1360.250117_rp ! Solar constant [W/m2]
  logical,  public :: atmos_solarins_fixedlatlon = .false.        ! Latitude/Longitude is fixed?
  logical,  public :: atmos_solarins_fixeddate   = .false.        ! Date is fixed?
  real(RP), public :: atmos_solarins_lon                          ! Longitude for radiation [rad]
  real(RP), public :: atmos_solarins_lat                          ! Latitude  for radiation [rad]
  integer,  public :: atmos_solarins_date(6)                      ! Date      for radiation [Y,M,D,H,M,S]

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private :: obliquity ! obliquity [rad]
  real(RP), private :: e         ! eccentricity
  real(RP), private :: omega     ! longitude of perigee [rad]
  real(RP), private :: lambda_m0 ! longitude at the vernal equinox [rad]

  integer,  private, parameter :: year_ref = 1950              ! reference year [year]
  integer,  private            :: ve_date(6)                   ! reference date of vernal equinox
  data ve_date / 1950, 3, 21, 0, 0, 0 /

  real(RP), private, parameter :: obliquity_ref = 23.320556_rp ! initial condition of obliquity (epsilon_star)
  real(RP), private, parameter :: psi_bar       = 50.439273_rp ! parameter for general precession [second of arc]
  real(RP), private, parameter :: zeta          =  3.392506_rp ! parameter for general precession [degree]

  !-----< Parameter tables from Berger(1978b) >-----
  integer,  private, parameter :: nobliq = 47         ! # of terms of the series expansion of epsilon
  real(RP), private            :: obliq_amp  (nobliq) ! amplitude [second of arc]
  real(RP), private            :: obliq_rate (nobliq) ! mean rate [second of arc/year]
  real(RP), private            :: obliq_phase(nobliq) ! phase     [degree of arc]

  data obliq_amp / &
        -2462.2214466_rp, & !  1
         -857.3232075_rp, & !  2
         -629.3231835_rp, & !  3
         -414.2804924_rp, & !  4
         -311.7632587_rp, & !  5
          308.9408604_rp, & !  6
         -162.5533601_rp, & !  7
         -116.1077911_rp, & !  8
          101.1189923_rp, & !  9
          -67.6856209_rp, & ! 10
           24.9079067_rp, & ! 11
           22.5811241_rp, & ! 12
          -21.1648355_rp, & ! 13
          -15.6549876_rp, & ! 14
           15.3936813_rp, & ! 15
           14.6660938_rp, & ! 16
          -11.7273029_rp, & ! 17
           10.2742696_rp, & ! 18
            6.4914588_rp, & ! 19
            5.8539148_rp, & ! 20
           -5.4872205_rp, & ! 21
           -5.4290191_rp, & ! 22
            5.1609570_rp, & ! 23
            5.0786314_rp, & ! 24
           -4.0735782_rp, & ! 25
            3.7227167_rp, & ! 26
            3.3971932_rp, & ! 27
           -2.8347004_rp, & ! 28
           -2.6550721_rp, & ! 29
           -2.5717867_rp, & ! 30
           -2.4712188_rp, & ! 31
            2.4625410_rp, & ! 32
            2.2464112_rp, & ! 33
           -2.0755511_rp, & ! 34
           -1.9713669_rp, & ! 35
           -1.8813061_rp, & ! 36
           -1.8468785_rp, & ! 37
            1.8186742_rp, & ! 38
            1.7601888_rp, & ! 39
           -1.5428851_rp, & ! 40
            1.4738838_rp, & ! 41
           -1.4593669_rp, & ! 42
            1.4192259_rp, & ! 43
           -1.1818980_rp, & ! 44
            1.1756474_rp, & ! 45
           -1.1316126_rp, & ! 46
            1.0896928_rp  / ! 47

  data obliq_rate / &
        31.609974_rp, & !  1
        32.620504_rp, & !  2
        24.172203_rp, & !  3
        31.983787_rp, & !  4
        44.828336_rp, & !  5
        30.973257_rp, & !  6
        43.668246_rp, & !  7
        32.246691_rp, & !  8
        30.599444_rp, & !  9
        42.681324_rp, & ! 10
        43.836462_rp, & ! 11!
        47.439436_rp, & ! 12
        63.219948_rp, & ! 13
        64.230478_rp, & ! 14
         1.010530_rp, & ! 15
         7.437771_rp, & ! 16
        55.782177_rp, & ! 17
         0.373813_rp, & ! 18
        13.218362_rp, & ! 19
        62.583231_rp, & ! 20
        63.593761_rp, & ! 21
        76.438310_rp, & ! 22
        45.815258_rp, & ! 23
         8.448301_rp, & ! 24
        56.792707_rp, & ! 25
        49.747842_rp, & ! 26
        12.058272_rp, & ! 27
        75.278220_rp, & ! 28
        65.241008_rp, & ! 29
        64.604291_rp, & ! 30
         1.647247_rp, & ! 31
         7.811584_rp, & ! 32
        12.207832_rp, & ! 33
        63.856665_rp, & ! 34
        56.155990_rp, & ! 35
        77.448840_rp, & ! 36
         6.801054_rp, & ! 37
        62.209418_rp, & ! 38
        20.656133_rp, & ! 39
        48.344406_rp, & ! 40
        55.145460_rp, & ! 41
        69.000539_rp, & ! 42
        11.071350_rp, & ! 43
        74.291298_rp, & ! 44
        11.047742_rp, & ! 45
         0.636717_rp, & ! 46
        12.844549_rp  / ! 47

  data obliq_phase / &
        251.9025_rp, & !  1
        280.8325_rp, & !  2
        128.3057_rp, & !  3
        292.7252_rp, & !  4
         15.3747_rp, & !  5
        263.7951_rp, & !  6
        308.4258_rp, & !  7
        240.0099_rp, & !  8
        222.9725_rp, & !  9
        268.7809_rp, & ! 10
        316.7998_rp, & ! 11
        319.6024_rp, & ! 12
        143.8050_rp, & ! 13
        172.7351_rp, & ! 14
         28.9300_rp, & ! 15
        123.5968_rp, & ! 16
         20.2082_rp, & ! 17
         40.8226_rp, & ! 18
        123.4722_rp, & ! 19
        155.6977_rp, & ! 20
        184.6277_rp, & ! 21
        267.2772_rp, & ! 22
         55.0196_rp, & ! 23
        152.5268_rp, & ! 24
         49.1382_rp, & ! 25
        204.6609_rp, & ! 26
         56.5233_rp, & ! 27
        200.3284_rp, & ! 28
        201.6651_rp, & ! 29
        213.5577_rp, & ! 30
         17.0374_rp, & ! 31
        164.4194_rp, & ! 32
         94.5422_rp, & ! 33
        131.9124_rp, & ! 34
         61.0309_rp, & ! 35
        296.2073_rp, & ! 36
        135.4894_rp, & ! 37
        114.8750_rp, & ! 38
        247.0691_rp, & ! 39
        256.6114_rp, & ! 40
         32.1008_rp, & ! 41
        143.6804_rp, & ! 42
         16.8784_rp, & ! 43
        160.6835_rp, & ! 44
         27.5932_rp, & ! 45
        348.1074_rp, & ! 46
         82.6496_rp  / ! 47

  integer,  private, parameter :: neclip = 19         ! # of terms of the series expansion of ecliptic
  real(RP), private            :: eclip_amp  (neclip) ! amplitude
  real(RP), private            :: eclip_rate (neclip) ! mean rate [second of arc/year]
  real(RP), private            :: eclip_phase(neclip) ! phase     [degree of arc]

  data eclip_amp / &
         0.01860798_rp, & !  1
         0.01627522_rp, & !  2
        -0.01300660_rp, & !  3
         0.00988829_rp, & !  4
        -0.00336700_rp, & !  5
         0.00333077_rp, & !  6
        -0.00235400_rp, & !  7
         0.00140015_rp, & !  8
         0.00100700_rp, & !  9
         0.00085700_rp, & ! 10
         0.00064990_rp, & ! 11
         0.00059900_rp, & ! 12
         0.00037800_rp, & ! 13
        -0.00033700_rp, & ! 14
         0.00027600_rp, & ! 15
         0.00018200_rp, & ! 16
        -0.00017400_rp, & ! 17
        -0.00012400_rp, & ! 18
         0.00001250_rp  / ! 19

  data eclip_rate / &
         4.2072050_rp, & !  1
         7.3460910_rp, & !  2
        17.8572630_rp, & !  3
        17.2205460_rp, & !  4
        16.8467330_rp, & !  5
         5.1990790_rp, & !  6
        18.2310760_rp, & !  7
        26.2167580_rp, & !  8
         6.3591690_rp, & !  9
        16.2100160_rp, & ! 10
         3.0651810_rp, & ! 11
        16.5838290_rp, & ! 12
        18.4939800_rp, & ! 13
         6.1909530_rp, & ! 14
        18.8677930_rp, & ! 15
        17.4255670_rp, & ! 16
         6.1860010_rp, & ! 17
        18.4174410_rp, & ! 18
         0.6678630_rp  / ! 19

  data eclip_phase / &
         28.620089_rp, & !  1
        193.788772_rp, & !  2
        308.307024_rp, & !  3
        320.199637_rp, & !  4
        279.376984_rp, & !  5
         87.195000_rp, & !  6
        349.129677_rp, & !  7
        128.443387_rp, & !  8
        154.143880_rp, & !  9
        291.269597_rp, & ! 10
        114.860583_rp, & ! 11
        332.092251_rp, & ! 12
        296.414411_rp, & ! 13
        145.769910_rp, & ! 14
        337.237063_rp, & ! 15
        152.092288_rp, & ! 16
        126.839891_rp, & ! 17
        210.667199_rp, & ! 18
         72.108838_rp  / ! 19

  integer,  private, parameter :: nprece = 78         ! # of terms of the series expansion of general precession
  real(RP), private            :: prece_amp  (nprece) ! amplitude
  real(RP), private            :: prece_rate (nprece) ! mean rate [second of arc/year]
  real(RP), private            :: prece_phase(nprece) ! phase     [degree of arc]

  data prece_amp / &
         7391.0225890_rp, & !  1
         2555.1526947_rp, & !  2
         2022.7629188_rp, & !  3
        -1973.6517951_rp, & !  4
         1240.2321818_rp, & !  5
          953.8679112_rp, & !  6
         -931.7537108_rp, & !  7
          872.3795383_rp, & !  8
          606.3544732_rp, & !  9
         -496.0274038_rp, & ! 10
          456.9608039_rp, & ! 11
          346.9462320_rp, & ! 12
         -305.8412902_rp, & ! 13
          249.6173246_rp, & ! 14
         -199.1027200_rp, & ! 15
          191.0560889_rp, & ! 16
         -175.2936572_rp, & ! 17
          165.9068833_rp, & ! 18
          161.1285917_rp, & ! 19
          139.7878093_rp, & ! 20
         -133.5228399_rp, & ! 21
          117.0673811_rp, & ! 22
          104.6907281_rp, & ! 23
           95.3227476_rp, & ! 24
           86.7824524_rp, & ! 25
           86.0857729_rp, & ! 26
           70.5893698_rp, & ! 27
          -69.9719343_rp, & ! 28
          -62.5817473_rp, & ! 29
           61.5450059_rp, & ! 30
          -57.9364011_rp, & ! 31
           57.1899832_rp, & ! 32
          -57.0236109_rp, & ! 33
          -54.2119253_rp, & ! 34
           53.2834147_rp, & ! 35
           52.1223575_rp, & ! 36
          -49.0059908_rp, & ! 37
          -48.3118757_rp, & ! 38
          -45.4191685_rp, & ! 39
          -42.2357920_rp, & ! 40
          -34.7971099_rp, & ! 41
           34.4623613_rp, & ! 42
          -33.8356643_rp, & ! 43
           33.6689362_rp, & ! 44
          -31.2521586_rp, & ! 45
          -30.8798701_rp, & ! 46
           28.4640769_rp, & ! 47
          -27.1960802_rp, & ! 48
           27.0860736_rp, & ! 49
          -26.3437456_rp, & ! 50
           24.7253740_rp, & ! 51
           24.6732126_rp, & ! 52
           24.4272733_rp, & ! 53
           24.0127327_rp, & ! 54
           21.7150294_rp, & ! 55
          -21.5375347_rp, & ! 56
           18.1148363_rp, & ! 57
          -16.9603104_rp, & ! 58
          -16.1765215_rp, & ! 59
           15.5567653_rp, & ! 60
           15.4846529_rp, & ! 61
           15.2150632_rp, & ! 62
           14.5047426_rp, & ! 63
          -14.3873316_rp, & ! 64
           13.1351419_rp, & ! 65
           12.8776311_rp, & ! 66
           11.9867234_rp, & ! 67
           11.9385578_rp, & ! 68
           11.7030822_rp, & ! 69
           11.6018181_rp, & ! 70
          -11.2617293_rp, & ! 71
          -10.4664199_rp, & ! 72
           10.4333970_rp, & ! 73
          -10.2377466_rp, & ! 74
           10.1934446_rp, & ! 75
          -10.1280191_rp, & ! 76
           10.0289441_rp, & ! 77
          -10.0034259_rp  / ! 78

  data prece_rate / &
        31.609974_rp, & !  1
        32.620504_rp, & !  2
        24.172203_rp, & !  3
         0.636717_rp, & !  4
        31.983787_rp, & !  5
         3.138886_rp, & !  6
        30.973257_rp, & !  7
        44.828336_rp, & !  8
         0.991874_rp, & !  9
         0.373813_rp, & ! 10
        43.668246_rp, & ! 11
        32.246691_rp, & ! 12
        30.599444_rp, & ! 13
         2.147012_rp, & ! 14
        10.511172_rp, & ! 15
        42.681324_rp, & ! 16
        13.650058_rp, & ! 17
         0.986922_rp, & ! 18
         9.874455_rp, & ! 19
        13.013341_rp, & ! 20
         0.262904_rp, & ! 21
         0.004952_rp, & ! 22
         1.142024_rp, & ! 23
        63.219948_rp, & ! 24
         0.205021_rp, & ! 25
         2.151964_rp, & ! 26
        64.230478_rp, & ! 27
        43.836462_rp, & ! 28
        47.439436_rp, & ! 29
         1.384343_rp, & ! 30
         7.437771_rp, & ! 31
        18.829299_rp, & ! 32
         9.500642_rp, & ! 33
         0.431696_rp, & ! 34
         1.160090_rp, & ! 35
        55.782177_rp, & ! 36
        12.639528_rp, & ! 37
         1.155138_rp, & ! 38
         0.168216_rp, & ! 39
         1.647247_rp, & ! 40
        10.884985_rp, & ! 41
         5.610937_rp, & ! 42
        12.658184_rp, & ! 43
         1.010530_rp, & ! 44
         1.983748_rp, & ! 45
        14.023871_rp, & ! 46
         0.560178_rp, & ! 47
         1.273434_rp, & ! 48
        12.021467_rp, & ! 49
        62.583231_rp, & ! 50
        63.593761_rp, & ! 51
        76.438310_rp, & ! 52
         4.280910_rp, & ! 53
        13.218362_rp, & ! 54
        17.818769_rp, & ! 55
         8.359495_rp, & ! 56
        56.792707_rp, & ! 57
         8.448301_rp, & ! 58
         1.978796_rp, & ! 59
         8.863925_rp, & ! 60
         0.186365_rp, & ! 61
         8.996212_rp, & ! 62
         6.771027_rp, & ! 63
        45.815258_rp, & ! 64
        12.002811_rp, & ! 65
        75.278220_rp, & ! 66
        65.241008_rp, & ! 67
        18.870667_rp, & ! 68
        22.009553_rp, & ! 69
        64.604291_rp, & ! 70
        11.498094_rp, & ! 71
         0.578834_rp, & ! 72
         9.237738_rp, & ! 73
        49.747842_rp, & ! 74
         2.147012_rp, & ! 75
         1.196895_rp, & ! 76
         2.133898_rp, & ! 77
         0.173168_rp  / ! 78

  data prece_phase / &
        251.9025_rp, & !  1
        280.8325_rp, & !  2
        128.3057_rp, & !  3
        348.1074_rp, & !  4
        292.7252_rp, & !  5
        165.1686_rp, & !  6
        263.7951_rp, & !  7
         15.3747_rp, & !  8
         58.5749_rp, & !  9
         40.8226_rp, & ! 10
        308.4258_rp, & ! 11
        240.0099_rp, & ! 12
        222.9725_rp, & ! 13
        106.5937_rp, & ! 14
        114.5182_rp, & ! 15
        268.7809_rp, & ! 16
        279.6869_rp, & ! 17
         39.6448_rp, & ! 18
        126.4108_rp, & ! 19
        291.5795_rp, & ! 20
        307.2848_rp, & ! 21
         18.9300_rp, & ! 22
        273.7596_rp, & ! 23
        143.8050_rp, & ! 24
        191.8927_rp, & ! 25
        125.5237_rp, & ! 26
        172.7351_rp, & ! 27
        316.7998_rp, & ! 28
        319.6024_rp, & ! 29
         69.7526_rp, & ! 30
        123.5968_rp, & ! 31
        217.6432_rp, & ! 32
         85.5882_rp, & ! 33
        156.2147_rp, & ! 34
         66.9489_rp, & ! 35
         20.2082_rp, & ! 36
        250.7568_rp, & ! 37
         48.0188_rp, & ! 38
          8.3739_rp, & ! 39
         17.0374_rp, & ! 40
        155.3409_rp, & ! 41
         94.1709_rp, & ! 42
        221.1120_rp, & ! 43
         28.9300_rp, & ! 44
        117.1498_rp, & ! 45
        320.5095_rp, & ! 46
        262.3602_rp, & ! 47
        336.2148_rp, & ! 48
        233.0046_rp, & ! 49
        155.6977_rp, & ! 50
        184.6277_rp, & ! 51
        267.2772_rp, & ! 52
         78.9281_rp, & ! 53
        123.4722_rp, & ! 54
        188.7132_rp, & ! 55
        180.1364_rp, & ! 56
         49.1382_rp, & ! 57
        152.5268_rp, & ! 58
         98.2198_rp, & ! 59
         97.4808_rp, & ! 60
        221.5376_rp, & ! 61
        168.2438_rp, & ! 62
        161.1199_rp, & ! 63
         55.0196_rp, & ! 64
        262.6495_rp, & ! 65
        200.3284_rp, & ! 66
        201.6651_rp, & ! 67
        294.6547_rp, & ! 68
         99.8233_rp, & ! 69
        213.5577_rp, & ! 70
        154.1631_rp, & ! 71
        232.7153_rp, & ! 72
        138.3034_rp, & ! 73
        204.6609_rp, & ! 74
        106.5938_rp, & ! 75
        250.4676_rp, & ! 76
        332.3345_rp, & ! 77
         27.3039_rp  / ! 78

  real(RP), private :: arcsec2d, arcsec2r ! unit converter

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_setup( &
       iyear )
    use scale_process, only: &
       prc_mpistop
    use scale_const, only: &
       const_d2r
    use scale_grid_real, only: &
       real_basepoint_lon, &
       real_basepoint_lat
    implicit none

    integer, intent(in) :: iyear ! year at setup

    namelist / param_atmos_solarins / &
       atmos_solarins_constant,    &
       atmos_solarins_fixedlatlon, &
       atmos_solarins_fixeddate,   &
       atmos_solarins_lon,         &
       atmos_solarins_lat,         &
       atmos_solarins_date

    real(RP) :: dyear ! delta t [year]
    integer  :: year  ! used year

    integer  :: ierr
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[SOLARINS] / Categ[ATMOS SHARE] / Origin[SCALElib]'

    atmos_solarins_lon     = real_basepoint_lon
    atmos_solarins_lat     = real_basepoint_lat
    atmos_solarins_date(:) = -1

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_solarins,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_SOLARINS. Check!'
       call prc_mpistop
    endif
    if( io_nml ) write(io_fid_nml,nml=param_atmos_solarins)

    arcsec2d = 1.0_rp / (60.0_rp*60.0_rp)
    arcsec2r = arcsec2d * const_d2r

    year = iyear

    if ( atmos_solarins_fixeddate ) then
       if( atmos_solarins_date(1) >= 0 ) year = atmos_solarins_date(1)
    endif

    call atmos_solarins_orbit( year )

    dyear = real( year-year_ref, kind=rp )

    !----- report data -----
    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,'(1x,A)')       '*** Solar insolation parameters ***'
    if( io_l ) write(io_fid_log,'(1x,A,I7)')    '*** Reference year                        : ', year_ref
    if( io_l ) write(io_fid_log,'(1x,A,I7)')    '*** Current   year                        : ', iyear
    if( io_l ) write(io_fid_log,'(1x,A,I7)')    '*** Difference from ref.                  : ', int(dyear)
    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** Solar constant                 [W/m2] : ', atmos_solarins_constant
    if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** Obliquity                       [deg] : ', obliquity / const_d2r
    if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** Eccentricity                          : ', e
    if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** Longitude of perihelion         [deg] : ', omega     / const_d2r
    if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** Longitude at the vernal equinox [deg] : ', lambda_m0 / const_d2r
    if( io_l ) write(io_fid_log,*)

    if( io_l ) write(io_fid_log,*)              '*** Latitude/Longitude is fixed?          : ', atmos_solarins_fixedlatlon
    if ( atmos_solarins_fixedlatlon ) then
       if( io_l ) write(io_fid_log,*)           '*** Longitude                       [deg] : ', atmos_solarins_lon
       if( io_l ) write(io_fid_log,*)           '*** Latitude                        [deg] : ', atmos_solarins_lat
       atmos_solarins_lon = atmos_solarins_lon * const_d2r
       atmos_solarins_lat = atmos_solarins_lat * const_d2r
    endif

    if( io_l ) write(io_fid_log,*)              '*** Date is fixed?                        : ', atmos_solarins_fixeddate
    if ( atmos_solarins_fixeddate ) then
       if( io_l ) write(io_fid_log,*)           '*** Date                                  : ', atmos_solarins_date
    endif

    return
  end subroutine atmos_solarins_setup

  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_orbit( &
       iyear )
    use scale_const, only: &
         pi => const_pi, &
         const_d2r
    implicit none

    integer, intent(in) :: iyear ! year at setup

    real(RP) :: dyear ! delta t [year]

    real(RP) :: esinomg     ! e * sin(w) : ecliptic parameter
    real(RP) :: ecosomg     ! e * cos(w) : ecliptic parameter
    real(RP) :: perih_pi    ! pi         : longitude of fixed perihelion [rad]
    real(RP) :: perih_psi   ! psi        : general precession [degree]
    real(RP) :: perih_omega ! omega_bar  : longitude of perihelion [degree]

    real(RP) :: temp, beta, ecosomg_mod
    integer  :: i
    !---------------------------------------------------------------------------

    ! time from reference year(1950.0 AD)
    dyear = real( iyear-year_ref, kind=rp )

    ! obliquity
    temp = 0.0_rp
    do i = 1, nobliq
       temp = temp + obliq_amp(i)*arcsec2d * cos( obliq_rate(i)*arcsec2r*dyear + obliq_phase(i)*const_d2r )
    enddo
    obliquity = ( obliquity_ref + temp ) * const_d2r

    ! eccentricity
    esinomg = 0.0_rp
    ecosomg = 0.0_rp
    do i = 1, neclip
       esinomg = esinomg + eclip_amp(i) * sin( eclip_rate(i)*arcsec2r*dyear + eclip_phase(i)*const_d2r )
       ecosomg = ecosomg + eclip_amp(i) * cos( eclip_rate(i)*arcsec2r*dyear + eclip_phase(i)*const_d2r )
    enddo
    e = sqrt( esinomg*esinomg + ecosomg*ecosomg )

    ! longitude of fixed perihelion
    if ( ecosomg == 0.0_rp ) then
       ecosomg_mod  = 1.e-30_rp
    else
       ecosomg_mod  = ecosomg
    endif
    perih_pi = atan(esinomg/ecosomg_mod) + pi * ( 0.5_rp - sign(0.5_rp, ecosomg) )

!    if( abs(EcosOMG) < 1.E-8_RP ) then
!       if    ( EsinOMG == 0.0_RP ) then
!          Perih_pi = 0.0_RP
!       elseif( EsinOMG <  0.0_RP ) then
!          Perih_pi = 1.5_RP * PI
!       elseif( EsinOMG >  0.0_RP ) then
!          Perih_pi = 0.5_RP * PI
!       endif
!
!    elseif( EcosOMG < 0.0_RP ) then
!
!       Perih_pi = atan(EsinOMG/EcosOMG) + PI
!
!    elseif( EcosOMG > 0.0_RP ) then
!
!       if ( EsinOMG <  0.0_RP ) then
!          Perih_pi = atan(EsinOMG/EcosOMG) + 2.0_RP * PI
!       else
!          Perih_pi = atan(EsinOMG/EcosOMG)
!       endif
!
!    endif
    perih_pi = perih_pi / const_d2r ! [rad]->[degree]

    ! general precession
    temp = 0.0_rp
    do i = 1, nprece
       temp = temp + prece_amp(i)*arcsec2d * sin( prece_rate(i)*arcsec2r*dyear + prece_phase(i)*const_d2r )
    enddo
    perih_psi = psi_bar * arcsec2d * dyear + zeta + temp

    ! longitude of perihelion
    perih_omega = perih_pi + perih_psi

    do i = 1, 1000
       if    ( perih_omega + 180.0_rp <    0.0_rp ) then
          perih_omega = perih_omega + 360.0_rp
       elseif( perih_omega + 180.0_rp >= 360.0_rp ) then
          perih_omega = perih_omega - 360.0_rp
       else
          exit
       endif
    enddo

    ! longitude of perigee (see berger et al.(1993))
    omega = ( perih_omega + 180.0_rp ) * const_d2r

    beta = sqrt( 1.0_rp - e*e )

    ! longitude at the vernal equinox
    lambda_m0 = 2.0_rp * ( ( 1.0_rp/2.0_rp*e + 1.0_rp/8.0_rp*e*e*e ) * ( 1.0_rp        + beta ) * sin(       omega) &
                         - (                   1.0_rp/4.0_rp*e*e   ) * ( 1.0_rp/2.0_rp + beta ) * sin(2.0_rp*omega) &
                         + (                   1.0_rp/8.0_rp*e*e*e ) * ( 1.0_rp/3.0_rp + beta ) * sin(3.0_rp*omega) )

    return
  end subroutine atmos_solarins_orbit

  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_insolation_0d( &
       solins,     &
       cosSZA,     &
       Re_factor,  &
       real_lon,   &
       real_lat,   &
       now_date,   &
       offset_year )
    use scale_const, only: &
       pi  => const_pi,  &
       eps => const_eps
    use scale_calendar, only: &
       calendar_getdayofyear,  &
       calendar_ymd2absday,    &
       calendar_hms2abssec,    &
       i_year, i_month, i_day, &
       i_hour, i_min, i_sec
    implicit none

    real(RP), intent(out) :: solins      ! solar insolation
    real(RP), intent(out) :: cosSZA      ! cos(Solar Zenith Angle)
    real(RP), intent(out) :: Re_factor   ! factor of the distance of Earth from the sun (1/rho2)
    real(RP), intent(in)  :: real_lon    ! longitude
    real(RP), intent(in)  :: real_lat    ! latitude
    integer,  intent(in)  :: now_date(6) ! date(yyyy,mm,dd,hh,mm,ss)
    integer,  intent(in)  :: offset_year ! year offset

    real(RP) :: lon
    real(RP) :: lat
    integer  :: date(6)
    integer  :: oyear

    real(RP) :: lambda_m       ! mean longitude from vernal equinox
    real(RP) :: lambda         !
    real(RP) :: sinDEC, cosDEC ! sin/cos(solar declination)
    real(RP) :: hourangle      ! hour angle: relative longitude of subsolar point

    integer  :: absday, absday_ve
    real(DP) :: DayOfYear, abssec
    real(RP) :: nu
    !---------------------------------------------------------------------------

    lon = real_lon
    lat = real_lat
    if ( atmos_solarins_fixedlatlon ) then
       lon = atmos_solarins_lon
       lat = atmos_solarins_lat
    endif

    date(:) = now_date(:)
    oyear   = offset_year
    if ( atmos_solarins_fixeddate ) then
       if( atmos_solarins_date(1) >= 0 ) oyear   = 0
       if( atmos_solarins_date(1) >= 0 ) date(1) = atmos_solarins_date(1)
       if( atmos_solarins_date(2) >= 1 ) date(2) = atmos_solarins_date(2)
       if( atmos_solarins_date(3) >= 1 ) date(3) = atmos_solarins_date(3)
       if( atmos_solarins_date(4) >= 0 ) date(4) = atmos_solarins_date(4)
       if( atmos_solarins_date(5) >= 0 ) date(5) = atmos_solarins_date(5)
       if( atmos_solarins_date(6) >= 0 ) date(6) = atmos_solarins_date(6)
    endif

    call calendar_getdayofyear( dayofyear, date(i_year) )

    call calendar_ymd2absday( absday,        & ! [OUT]
                              date(i_year),  & ! [IN]
                              date(i_month), & ! [IN]
                              date(i_day),   & ! [IN]
                              oyear          ) ! [IN]

    call calendar_ymd2absday( absday_ve,        & ! [OUT]
                              date(i_year),  & ! [IN]
                              ve_date(i_month), & ! [IN]
                              ve_date(i_day),   & ! [IN]
                              oyear             ) ! [IN]

    call calendar_hms2abssec( abssec,        & ! [OUT]
                              date(i_hour),  & ! [IN]
                              date(i_min),   & ! [IN]
                              date(i_sec),   & ! [IN]
                              0.0_dp         ) ! [IN]

    lambda_m = lambda_m0 + 2.0_rp * pi * real(absday-absday_ve,kind=RP) / DayOfYear

    nu = lambda_m - omega

    ! actual longitude from vernal equinox
    lambda = lambda_m &
           + ( 2.0_rp*e -  1.0_rp/ 4.0_rp*e*e*e ) * sin(        nu ) &
           + (             5.0_rp/ 4.0_rp*e*e   ) * sin( 2.0_rp*nu ) &
           + (            13.0_rp/12.0_rp*e*e*e ) * sin( 3.0_rp*nu )

    ! solar declination
    sindec = sin(lambda) * sin(obliquity)
    cosdec = sqrt( 1.0_rp - sindec*sindec )

    ! hour angle
    hourangle = lon + 2.0_rp * pi * abssec / (24.0_rp*60.0_rp*60.0_rp)

    ! cos(Solar Zenith Angle)
    cossza = sin(lat)*sindec - cos(lat)*cosdec*cos(hourangle)

    ! 1 / (rho*rho)
    re_factor = 1.0_rp                          &
              + 2.0_rp*e     * cos(        nu ) &
              + 0.5_rp*e*e   * cos( 2.0_rp*nu ) &
              + 2.5_rp*e*e                      &
              + 4.0_rp*e*e*e * cos( 3.0_rp*nu )

!    solins = ATMOS_SOLARINS_constant * Re_factor * max(cosSZA,0.0_RP)
    solins = atmos_solarins_constant * re_factor * ( 0.5_rp + sign(0.5_rp,cossza-eps) )

    return
  end subroutine atmos_solarins_insolation_0d

  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_insolation_2d( &
       solins,     &
       cosSZA,     &
       real_lon,   &
       real_lat,   &
       now_date,   &
       offset_year )
    use scale_grid_index
    use scale_const, only: &
       pi  => const_pi,  &
       eps => const_eps
    use scale_calendar, only: &
       calendar_getdayofyear,  &
       calendar_ymd2absday,    &
       calendar_hms2abssec,    &
       i_year, i_month, i_day, &
       i_hour, i_min, i_sec
    implicit none

    real(RP), intent(out) :: solins  (IA,JA) ! solar insolation
    real(RP), intent(out) :: cosSZA  (IA,JA) ! cos(Solar Zenith Angle)
    real(RP), intent(in)  :: real_lon(IA,JA) ! longitude [rad]
    real(RP), intent(in)  :: real_lat(IA,JA) ! latitude  [rad]
    integer,  intent(in)  :: now_date(6)     ! date(yyyy,mm,dd,hh,mm,ss)
    integer,  intent(in)  :: offset_year     ! year offset

    real(RP) :: lon(IA,JA)
    real(RP) :: lat(IA,JA)
    integer  :: date(6)
    integer  :: oyear

    real(RP) :: lambda_m       ! mean longitude from vernal equinox
    real(RP) :: lambda         !
    real(RP) :: sinDEC, cosDEC ! sin/cos(solar declination)
    real(RP) :: hourangle(IA,JA) ! hour angle: relative longitude of subsolar point

    integer  :: absday, absday_ve
    real(DP) :: DayOfYear, abssec
    real(RP) :: nu
    real(RP) :: Re_factor   ! factor of the distance of Earth from the sun (1/rho2)

    integer  :: i, j
    !---------------------------------------------------------------------------

    lon(:,:) = real_lon(:,:)
    lat(:,:) = real_lat(:,:)
    if ( atmos_solarins_fixedlatlon ) then
       lon(:,:) = atmos_solarins_lon
       lat(:,:) = atmos_solarins_lat
    endif

    date(:) = now_date(:)
    oyear   = offset_year
    if ( atmos_solarins_fixeddate ) then
       if( atmos_solarins_date(1) >= 0 ) oyear   = 0
       if( atmos_solarins_date(1) >= 0 ) date(1) = atmos_solarins_date(1)
       if( atmos_solarins_date(2) >= 1 ) date(2) = atmos_solarins_date(2)
       if( atmos_solarins_date(3) >= 1 ) date(3) = atmos_solarins_date(3)
       if( atmos_solarins_date(4) >= 0 ) date(4) = atmos_solarins_date(4)
       if( atmos_solarins_date(5) >= 0 ) date(5) = atmos_solarins_date(5)
       if( atmos_solarins_date(6) >= 0 ) date(6) = atmos_solarins_date(6)
    endif

    call calendar_getdayofyear( dayofyear, date(i_year) )

    call calendar_ymd2absday( absday,        & ! [OUT]
                              date(i_year),  & ! [IN]
                              date(i_month), & ! [IN]
                              date(i_day),   & ! [IN]
                              oyear          ) ! [IN]

    call calendar_ymd2absday( absday_ve,        & ! [OUT]
                              date(i_year),  & ! [IN]
                              ve_date(i_month), & ! [IN]
                              ve_date(i_day),   & ! [IN]
                              oyear             ) ! [IN]

    call calendar_hms2abssec( abssec,       & ! [OUT]
                              date(i_hour), & ! [IN]
                              date(i_min),  & ! [IN]
                              date(i_sec),  & ! [IN]
                              0.0_dp        ) ! [IN]

    lambda_m = lambda_m0 + 2.0_rp * pi * real(absday-absday_ve,kind=RP) / DayOfYear

    nu = lambda_m - omega

    ! 1 / (rho*rho)
    re_factor = 1.0_rp                          &
              + 2.0_rp*e     * cos(        nu ) &
              + 0.5_rp*e*e   * cos( 2.0_rp*nu ) &
              + 2.5_rp*e*e                      &
              + 4.0_rp*e*e*e * cos( 3.0_rp*nu )

    ! actual longitude from vernal equinox
    lambda = lambda_m &
           + ( 2.0_rp*e -  1.0_rp/ 4.0_rp*e*e*e ) * sin(        nu ) &
           + (             5.0_rp/ 4.0_rp*e*e   ) * sin( 2.0_rp*nu ) &
           + (            13.0_rp/12.0_rp*e*e*e ) * sin( 3.0_rp*nu )

    ! solar declination
    sindec = sin(lambda) * sin(obliquity)
    cosdec = sqrt( 1.0_rp - sindec*sindec )

    do j = js, je
    do i = is, ie
       ! hour angle
       hourangle(i,j) = lon(i,j) + 2.0_rp * pi * abssec / (24.0_rp*60.0_rp*60.0_rp)
    enddo
    enddo

    do j = js, je
    do i = is, ie
       ! cos(Solar Zenith Angle)
       cossza(i,j) = sin(lat(i,j))*sindec - cos(lat(i,j))*cosdec*cos(hourangle(i,j))
    enddo
    enddo

    do j = js, je
    do i = is, ie
!       solins(i,j) = ATMOS_SOLARINS_constant * Re_factor * max(cosSZA(i,j),0.0_RP)
       solins(i,j) = atmos_solarins_constant * re_factor * ( 0.5_rp + sign(0.5_rp,cossza(i,j)-eps) )
    enddo
    enddo

    return
  end subroutine atmos_solarins_insolation_2d

end module scale_atmos_solarins