scale_atmos_solarins.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_solarins
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_solarins_setup
  public :: atmos_solarins_orbit
  public :: atmos_solarins_insolation
  public :: atmos_solarins_ecliptic_longitude
 
  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_set_ve       = .false.        ! Set vernal equinox condition?
 
  logical,  public :: atmos_solarins_set_ideal    = .false.        ! Set obliquity and eccentricity?
  real(rp), public :: atmos_solarins_obliquity    = 23.44_rp       ! Obliquity [deg] cf. https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
  real(rp), public :: atmos_solarins_eccentricity = 0.0167_rp      ! Eccentricity
  real(rp), public :: atmos_solarins_perihelion_lon = 282.94719_rp ! Longitude of perihelion [deg].
  integer,  public :: atmos_solarins_ve_date(6)                    ! Date of first vernal equinox
  data atmos_solarins_ve_date / 2000, 3, 20, 7, 35, 0 /
  real(rp), public :: atmos_solarins_diurnal_sec                   ! Seconds of the diurnal period [sec]
  real(rp), public :: atmos_solarins_annual_sec                    ! Seconds of the annual period [sec]
 
  logical,  public :: atmos_solarins_fixedlatlon  = .false.        ! Latitude/Longitude is fixed?
  real(rp), public :: atmos_solarins_lon                           ! Longitude for radiation [rad]
  real(rp), public :: atmos_solarins_lat                           ! Latitude  for radiation [rad]
 
  logical,  public :: atmos_solarins_fixeddate    = .false.        ! Date is fixed?
  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 perihelion [rad]
  real(rp), private :: lambda_m0 ! longitude at the vernal equinox associated with circle orbit [rad]
 
  integer,  private, parameter :: year_ref = 1950              ! reference year [year]
  integer,  private            :: ve_date(6)                   ! reference date of vernal equinox
  data ve_date / 2000, 3, 20, 7, 35, 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]
 
  real(dp), private            :: caldaysec                     ! seconds of a calendar day
  real(dp), private            :: calyearday                    ! days of a calendar year
  real(dp), private            :: soldaysec                     ! seconds of a solar day
  real(dp), private            :: solyearsec                    ! seconds of a solar year
  logical,  private            :: day_calender_synced  = .true. ! solar diurnal cycle is synced with calendar?
  logical,  private            :: year_calender_synced = .true. ! solar annual cycle is synced with calendar?
 
 
  !-----< 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
 
  logical,  private :: debug = .false.
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_setup( &
       basepoint_lon, basepoint_lat, &
       iyear                         )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       const_d2r,          &
       const_undef
    use scale_calendar, only: &
       calendar_doi,          &
       calendar_hour,         &
       calendar_min,          &
       calendar_sec,          &
       calendar_daysec2date,  &
       calendar_date2daysec,  &
       calendar_date2char
    use scale_time, only: &
       time_offset_year,  &
       time_nowday,       &
       time_nowsec
    implicit none
    integer,  intent(in) :: iyear ! year at setup
    real(rp), intent(in) :: basepoint_lon
    real(rp), intent(in) :: basepoint_lat
 
    namelist / param_atmos_solarins / &
       atmos_solarins_constant,       &
       atmos_solarins_set_ve,         &
       atmos_solarins_set_ideal,      &
       atmos_solarins_obliquity,      &
       atmos_solarins_eccentricity,   &
       atmos_solarins_fixedlatlon,    &
       atmos_solarins_lon,            &
       atmos_solarins_lat,            &
       atmos_solarins_fixeddate,      &
       atmos_solarins_date,           &
       atmos_solarins_perihelion_lon, &
       atmos_solarins_ve_date,        &
       atmos_solarins_annual_sec,     &
       atmos_solarins_diurnal_sec,    &
       debug
 
    real(rp) :: dyear ! delta t [year]
    integer  :: year  ! used year
    integer  :: d
    integer  :: date(6)
    real(rp) :: re_factor, sindec, cosdec, hourangle
    real(dp) :: lambda, subsec
    character(len=27) :: datechar
 
    integer  :: ierr
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_SOLARINS_setup",*) 'Setup'
 
    atmos_solarins_lon         = basepoint_lon
    atmos_solarins_lat         = basepoint_lat
    atmos_solarins_date(:)     = -1
    atmos_solarins_diurnal_sec = const_undef
    atmos_solarins_annual_sec  = const_undef
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_solarins,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_SOLARINS_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_SOLARINS_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_SOLARINS. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_solarins)
 
    if ( atmos_solarins_set_ve ) then
       if ( .NOT. atmos_solarins_set_ideal ) then
          atmos_solarins_set_ideal    = .true.
          atmos_solarins_obliquity    = 0.0_rp
          atmos_solarins_eccentricity = 0.0_rp
          atmos_solarins_perihelion_lon  = 0.0_rp
       endif
       if ( .NOT. atmos_solarins_fixedlatlon ) then
          atmos_solarins_fixedlatlon  = .true.
          atmos_solarins_lon          = 0.0_rp
          atmos_solarins_lat          = 0.0_rp
       endif
       if ( .NOT. atmos_solarins_fixeddate ) then
          atmos_solarins_fixeddate    = .true.
          atmos_solarins_date(1)      = ve_date(1)
          atmos_solarins_date(2)      = ve_date(2)
          atmos_solarins_date(3)      = ve_date(3)
          atmos_solarins_date(4)      = 12
          atmos_solarins_date(5)      = 0
          atmos_solarins_date(6)      = 0
       endif
    endif
 
    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
 
    if ( atmos_solarins_set_ideal ) then
       ve_date = atmos_solarins_ve_date
       calyearday = calendar_doi
       caldaysec  = calendar_hour*calendar_min*calendar_sec
 
       if (atmos_solarins_annual_sec == const_undef ) then
         solyearsec = calyearday*caldaysec
         atmos_solarins_annual_sec = solyearsec
       else
         solyearsec = atmos_solarins_annual_sec
         year_calender_synced = .false.
       endif
 
       if (atmos_solarins_diurnal_sec == const_undef ) then
         soldaysec  = caldaysec
         atmos_solarins_diurnal_sec = soldaysec
       else
         soldaysec  = atmos_solarins_diurnal_sec
         day_calender_synced = .false.
       endif
 
    else
       calyearday = 365.2425_dp
       caldaysec  = calendar_hour*calendar_min*calendar_sec
       solyearsec = calyearday*caldaysec
       soldaysec  = caldaysec
    endif
 
    call atmos_solarins_orbit( year )
 
    dyear = real( year-year_ref, kind=rp )
 
    !----- report data -----
    log_newline
    log_info("ATMOS_SOLARINS_setup",'(1x,A)')       'Solar insolation parameters '
    log_info("ATMOS_SOLARINS_setup",'(1x,A,I7)')    'Reference year                        : ', year_ref
    log_info("ATMOS_SOLARINS_setup",'(1x,A,I7)')    'Current   year                        : ', year
    log_info("ATMOS_SOLARINS_setup",'(1x,A,I7)')    'Difference from ref.                  : ', int(dyear)
    log_newline
    log_info("ATMOS_SOLARINS_setup",'(1x,A,F12.7)') 'Solar constant                 [W/m2] : ', atmos_solarins_constant
    log_info("ATMOS_SOLARINS_setup",'(1x,A,F12.7)') 'Obliquity                       [deg] : ', obliquity / const_d2r
    log_info("ATMOS_SOLARINS_setup",'(1x,A,F12.7)') 'Eccentricity                          : ', e
    log_info("ATMOS_SOLARINS_setup",'(1x,A,F12.7)') 'Longitude of perihelion         [deg] : ', omega     / const_d2r
!    LOG_INFO("ATMOS_SOLARINS_setup",'(1x,A,F12.7)') 'Longitude at the vernal equinox associated with circle orbit [deg] : ', lambda_m0 / CONST_D2R
    log_newline
    if ( year_calender_synced ) then
       log_info("ATMOS_SOLARINS_setup",'(1x,A      )') 'Solar annual cycle is synced with calendar'
    else
       log_info("ATMOS_SOLARINS_setup",'(1x,A      )') 'Solar annual cycle is not synced with calendar'
    endif
    log_info("ATMOS_SOLARINS_setup",'(1x,A,F14.4)') 'Solar annual cycle              [sec] : ', solyearsec
    log_info("ATMOS_SOLARINS_setup",'(1x,A,F14.4)') 'Solar annual cycle              [day] : ', solyearsec/caldaysec
    if ( day_calender_synced ) then
       log_info("ATMOS_SOLARINS_setup",'(1x,A      )') 'Solar diurnal cycle is synced with calendar'
    else
       log_info("ATMOS_SOLARINS_setup",'(1x,A      )') 'Solar diurnal cycle is not synced with calendar'
    endif
    log_info("ATMOS_SOLARINS_setup",'(1x,A,F14.4)') 'Solar diurnal cycle             [sec] : ', soldaysec
    log_newline
    log_info("ATMOS_SOLARINS_setup",*)              'Latitude/Longitude is fixed?          : ', atmos_solarins_fixedlatlon
    if ( atmos_solarins_fixedlatlon ) then
       log_info("ATMOS_SOLARINS_setup",*)           'Longitude                       [deg] : ', atmos_solarins_lon
       log_info("ATMOS_SOLARINS_setup",*)           'Latitude                        [deg] : ', atmos_solarins_lat
       atmos_solarins_lon = atmos_solarins_lon * const_d2r
       atmos_solarins_lat = atmos_solarins_lat * const_d2r
    endif
 
    log_info("ATMOS_SOLARINS_setup",*)              'Date is fixed?                        : ', atmos_solarins_fixeddate
    if ( atmos_solarins_fixeddate ) then
       log_info("ATMOS_SOLARINS_setup",*)           'Date                                  : ', atmos_solarins_date
    endif
 
    log_newline
    log_info("ATMOS_SOLARINS_setup",*)  'Orbit info for next 1 solar year'
    if (io_l) write(io_fid_log, '(A)') 'DATE                         SOLAR_CONST, Scaled dist from star, Longitude'
    do d = time_nowday, time_nowday + int(solyearsec/caldaysec)
       call calendar_daysec2date(date,subsec,d,time_nowsec,time_offset_year)
       call atmos_solarins_ecliptic_longitude(  &
          re_factor,                            & ![OUT]
          sindec,                               & ![OUT]
          cosdec,                               & ![OUT]
          hourangle,                            & ![OUT]
          date,                                 & ![IN]
          time_offset_year,                     & ![IN]
          lambda_out = lambda                   & ![optional OUT]
       )
       call calendar_date2char(datechar, date, subsec)
       if (lambda < 0.0_rp) then
          lambda = lambda/const_d2r + 360.0_rp
       else
          lambda = lambda/const_d2r
       endif
       if (io_l) then
          write(io_fid_log,'(2A,F11.4,F23.6,F11.4)') datechar, "  ", re_factor*atmos_solarins_constant, 1.0_dp/sqrt(re_factor), lambda
       endif
    enddo
 
    return
  end subroutine atmos_solarins_setup
 
  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_orbit( &
       iyear )
    use scale_const, only: &
       pi  => const_pi, &
       d2r => const_d2r,&
       const_undef
    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
    !---------------------------------------------------------------------------
 
    if ( .not. atmos_solarins_set_ideal ) then ! real Earth orbit
       ! 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)*d2r )
       enddo
       obliquity = ( obliquity_ref + temp ) * 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)*d2r )
          ecosomg = ecosomg + eclip_amp(i) * cos( eclip_rate(i)*arcsec2r*dyear + eclip_phase(i)*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 / 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)*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 perihelion (see berger et al.(1993))
       omega = ( perih_omega + 180.0_rp ) * d2r
 
       atmos_solarins_eccentricity    = e
       atmos_solarins_obliquity       = obliquity / d2r
       atmos_solarins_perihelion_lon  = omega / d2r
 
    else ! ideal orbit
 
      obliquity = atmos_solarins_obliquity * d2r
      e         = atmos_solarins_eccentricity
      omega     = atmos_solarins_perihelion_lon * d2r
 
    endif
 
    beta = sqrt( 1.0_rp - e*e )
 
    ! longitude at the vernal equinox associated with circle orbit
    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_ecliptic_longitude( &
       Re_factor,  &
       sinDEC,     &
       cosDEC,     &
       hourangle,  &
       now_date,   &
       offset_year,&
       lambda_out  &
       )
    use scale_const, only: &
       pi => const_pi
    use scale_calendar, only: &
       calendar_ymd2absday,    &
       calendar_hms2abssec,    &
       i_year, i_month, i_day, &
       i_hour, i_min, i_sec
    implicit none
 
    real(rp), intent(out) :: re_factor   ! factor of the distance of Earth from the sun (1/rho2)
    real(rp), intent(out) :: sindec      ! sin/cos(solar declination)
    real(rp), intent(out) :: cosdec      ! sin/cos(solar declination)
    real(rp), intent(out) :: hourangle   ! hour angle: relative longitude of subsolar point (noon=0deg)
    integer,  intent(in)  :: now_date(6) ! date(yyyy,mm,dd,hh,mm,ss)
    integer,  intent(in)  :: offset_year ! year offset
    real(dp), intent(out), optional :: lambda_out   ! actual longitude from vernal equinox (for output)
 
 
    integer  :: date(6)
    integer  :: oyear
    integer  :: absday, absday_ve
    real(dp) :: abssec, abssec_ve
    real(dp) :: diffday, diffyr, diffsec
 
    real(dp) :: lambda_m ! mean   longitude from vernal equinox
    real(dp) :: lambda   ! actual longitude from vernal equinox
    real(dp) :: nu
 
    !---------------------------------------------------------------------------
 
    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_ymd2absday( absday,           & ! [OUT]
                              date(i_year),     & ! [IN]
                              date(i_month),    & ! [IN]
                              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]
 
    call calendar_ymd2absday( absday_ve,        & ! [OUT]
                              ve_date(i_year),  & ! [IN]
                              ve_date(i_month), & ! [IN]
                              ve_date(i_day),   & ! [IN]
                              oyear             ) ! [IN]
 
    call calendar_hms2abssec( abssec_ve,        & ! [OUT]
                              ve_date(i_hour),  & ! [IN]
                              ve_date(i_min),   & ! [IN]
                              ve_date(i_sec),   & ! [IN]
                              0.0_dp            ) ! [IN]
 
    if ( year_calender_synced ) then
       diffday = real(absday-absday_ve,kind=dp) + (abssec-abssec_ve) / caldaysec
       diffyr  = diffday / calyearday - real( nint( diffday / calyearday ),kind=dp)
    else
       diffsec = real(absday-absday_ve,kind=dp)*caldaysec + (abssec-abssec_ve)
       diffyr  = diffsec / solyearsec - real( nint( diffsec / solyearsec ),kind=dp)
    endif
 
    lambda_m = lambda_m0 + 2.0_rp * pi * diffyr
 
    nu = lambda_m - omega
 
    ! 1 / (rho*rho)
    re_factor = 1.0_dp                          &
              + 2.0_dp*e     * cos(        nu ) &
              + 0.5_dp*e*e   * cos( 2.0_dp*nu ) &
              + 2.5_dp*e*e                      &
              + 4.0_dp*e*e*e * cos( 3.0_dp*nu )
 
    ! actual longitude from vernal equinox
    lambda = lambda_m &
           + ( 2.0_dp*e -  1.0_dp/ 4.0_dp*e*e*e ) * sin(        nu ) &
           + (             5.0_dp/ 4.0_dp*e*e   ) * sin( 2.0_dp*nu ) &
           + (            13.0_dp/12.0_dp*e*e*e ) * sin( 3.0_dp*nu )
 
    ! solar declination
    sindec = sin(lambda) * sin(obliquity)
    cosdec = sqrt( 1.0_rp - sindec*sindec )
 
    ! hour angle
    if ( day_calender_synced ) then
       hourangle = pi + 2.0_rp * pi * abssec / caldaysec
    else
       hourangle = pi + 2.0_rp * pi                                                    &
                      * ( (absday*caldaysec + abssec) / atmos_solarins_diurnal_sec     &
                        - real( nint((absday*caldaysec + abssec)/ soldaysec ),kind=dp) )
    endif
 
    if ( debug ) then
       log_info("ATMOS_SOLARINS_ecliptic_longitude",*) lambda_m, nu, lambda, re_factor
    endif
 
    if ( present(lambda_out) ) lambda_out = lambda
 
    return
  end subroutine atmos_solarins_ecliptic_longitude
 
  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_insolation_0d( &
       real_lon, real_lat,    &
       now_date, offset_year, &
       solins, cosSZA,        &
       Re_factor_out )
    use scale_const, only: &
       eps => const_eps
    implicit none
    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), intent(out) :: solins      ! solar insolation
    real(RP), intent(out) :: cosSZA      ! cos(Solar Zenith Angle)
    real(RP), intent(out), optional :: Re_factor_out ! factor of the distance of Earth from the sun (1/rho2) [for output]
 
    real(RP) :: Re_factor      ! factor of the distance of Earth from the sun (1/rho2)
    real(RP) :: sinDEC, cosDEC ! sin/cos(solar declination)
    real(RP) :: hourangle      ! hour angle: relative longitude of subsolar point (noon=0deg)
 
    real(RP) :: lon
    real(RP) :: lat
    !---------------------------------------------------------------------------
 
    call atmos_solarins_ecliptic_longitude( re_factor,   & ! [OUT]
                                            sindec,      & ! [OUT]
                                            cosdec,      & ! [OUT]
                                            hourangle,   & ! [OUT]
                                            now_date(:), & ! [IN]
                                            offset_year  ) ! [IN]
 
    if ( atmos_solarins_fixedlatlon ) then
       lon = atmos_solarins_lon
       lat = atmos_solarins_lat
    else
       lon = real_lon
       lat = real_lat
    endif
 
    cossza = sin(lat)*sindec + cos(lat)*cosdec*cos(lon+hourangle)
    solins = atmos_solarins_constant * re_factor * ( 0.5_rp + sign(0.5_rp,cossza-eps) )
 
    if ( present( re_factor_out ) ) then
       re_factor_out = re_factor
    endif
 
    return
  end subroutine atmos_solarins_insolation_0d
 
  !-----------------------------------------------------------------------------
  subroutine atmos_solarins_insolation_2d( &
       IA, IS, IE, JA, JS, JE, &
       real_lon, real_lat,    &
       now_date, offset_year, &
       solins, cosSZA )
    use scale_const, only: &
       eps => const_eps
    implicit none
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    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), intent(out) :: solins  (IA,JA) ! solar insolation
    real(RP), intent(out) :: cosSZA  (IA,JA) ! cos(Solar Zenith Angle)
 
    real(RP) :: Re_factor      ! factor of the distance of Earth from the sun (1/rho2)
    real(RP) :: sinDEC, cosDEC ! sin/cos(solar declination)
    real(RP) :: hourangle      ! hour angle: relative longitude of subsolar point (noon=0deg)
 
    real(RP) :: lon(IA,JA)
    real(RP) :: lat(IA,JA)
 
    integer  :: i, j
    !---------------------------------------------------------------------------
 
    call atmos_solarins_ecliptic_longitude( re_factor,   & ! [OUT]
                                            sindec,      & ! [OUT]
                                            cosdec,      & ! [OUT]
                                            hourangle,   & ! [OUT]
                                            now_date(:), & ! [IN]
                                            offset_year  ) ! [IN]
 
    !$acc data create(lon, lat)
 
    if ( atmos_solarins_fixedlatlon ) then
       !$omp parallel do OMP_SCHEDULE_ collapse(2)
       !$acc kernels
       do j = js, je
       do i = is, ie
          lon(i,j) = atmos_solarins_lon
          lat(i,j) = atmos_solarins_lat
       enddo
       enddo
       !$acc end kernels
    else
       !$omp parallel do OMP_SCHEDULE_ collapse(2)
       !$acc kernels
       do j = js, je
       do i = is, ie
          lon(i,j) = real_lon(i,j)
          lat(i,j) = real_lat(i,j)
       enddo
       enddo
       !$acc end kernels
    endif
 
    !$omp parallel do OMP_SCHEDULE_ collapse(2)
    !$acc kernels
    do j = js, je
    do i = is, ie
       cossza(i,j) = sin(lat(i,j))*sindec + cos(lat(i,j))*cosdec*cos(lon(i,j)+hourangle)
       solins(i,j) = atmos_solarins_constant * re_factor * ( 0.5_rp + sign(0.5_rp,cossza(i,j)-eps) )
    enddo
    enddo
    !$acc end kernels
 
    !$acc end data
 
    return
  end subroutine atmos_solarins_insolation_2d
 
end module scale_atmos_solarins