LCOV - code coverage report
Current view: top level - columnphysics - icepack_orbital.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 84 175 48.00 %
Date: 2020-07-03 20:11:40 Functions: 4 5 80.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : 
       3             : ! Orbital parameters computed from date
       4             : ! author:  Bruce P. Briegleb, NCAR 
       5             : !
       6             : ! 2006 ECH: Converted to free source form (F90)
       7             : ! 2014 ECH: Moved routines from csm_share/shr_orb_mod.F90
       8             : 
       9             :       module icepack_orbital
      10             : 
      11             :       use icepack_kinds
      12             :       use icepack_parameters, only: c2, p5, pi, secday
      13             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      14             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      15             : 
      16             :       implicit none
      17             :       private
      18             :       public :: compute_coszen, &
      19             :                 icepack_init_orbit, &
      20             :                 icepack_query_orbit
      21             : 
      22             :       ! orbital parameters
      23             :       integer (kind=int_kind)  :: iyear_AD  ! Year to calculate orbit for
      24             :       real(kind=dbl_kind) :: eccen  ! Earth's orbital eccentricity
      25             :       real(kind=dbl_kind) :: obliqr ! Earth's obliquity in radians
      26             :       real(kind=dbl_kind) :: lambm0 ! Mean longitude of perihelion at the
      27             :                                            ! vernal equinox (radians)
      28             :       real(kind=dbl_kind) :: mvelpp ! Earth's moving vernal equinox longitude
      29             :                                            ! of perihelion + pi (radians)
      30             :       real(kind=dbl_kind) :: obliq  ! obliquity in degrees
      31             :       real(kind=dbl_kind) :: mvelp  ! moving vernal equinox long
      32             :       real(kind=dbl_kind) :: decln  ! solar declination angle in radians
      33             :       real(kind=dbl_kind) :: eccf   ! earth orbit eccentricity factor
      34             :       logical(kind=log_kind) :: log_print ! Flags print of status/error
      35             : 
      36             : !=======================================================================
      37             :  
      38             :       contains
      39             : 
      40             : !=======================================================================
      41             : 
      42             : !autodocument_start icepack_init_orbit
      43             : ! Compute orbital parameters for the specified date.
      44             : 
      45         156 :       subroutine icepack_init_orbit(iyear_AD_in, eccen_in, obliqr_in, &
      46             :          lambm0_in, mvelpp_in, obliq_in, mvelp_in, decln_in, eccf_in, &
      47             :          log_print_in)
      48             : 
      49             :       integer(kind=int_kind), optional, intent(in) :: iyear_AD_in  ! Year to calculate orbit for
      50             :       real(kind=dbl_kind), optional, intent(in) :: eccen_in  ! Earth's orbital eccentricity
      51             :       real(kind=dbl_kind), optional, intent(in) :: obliqr_in ! Earth's obliquity in radians
      52             :       real(kind=dbl_kind), optional, intent(in) :: lambm0_in ! Mean longitude of perihelion at the
      53             :                                                              ! vernal equinox (radians)
      54             :       real(kind=dbl_kind), optional, intent(in) :: mvelpp_in ! Earth's moving vernal equinox longitude
      55             :                                                              ! of perihelion + pi (radians)
      56             :       real(kind=dbl_kind), optional, intent(in) :: obliq_in  ! obliquity in degrees
      57             :       real(kind=dbl_kind), optional, intent(in) :: mvelp_in  ! moving vernal equinox long
      58             :       real(kind=dbl_kind), optional, intent(in) :: decln_in  ! solar declination angle in radians
      59             :       real(kind=dbl_kind), optional, intent(in) :: eccf_in   ! earth orbit eccentricity factor
      60             :       logical(kind=log_kind), optional, intent(in) :: log_print_in ! Flags print of status/error
      61             : 
      62             : !autodocument_end
      63             : 
      64             :       character(len=*),parameter :: subname='(icepack_init_orbit)'
      65             : 
      66             :       !call icepack_warnings_add(subname//'  ')
      67         156 :       iyear_AD  = 1950
      68         156 :       log_print = .false.   ! if true, write out orbital parameters
      69             : 
      70             : #ifndef CESMCOUPLED
      71             :       call shr_orb_params( iyear_AD, eccen , obliq , mvelp    , &
      72         156 :                            obliqr  , lambm0, mvelpp, log_print)
      73         156 :       if (icepack_warnings_aborted(subname)) return
      74             : #endif
      75             : 
      76         156 :       if (present(iyear_AD_in)) iyear_AD = iyear_AD_in
      77         156 :       if (present(eccen_in)) eccen = eccen_in
      78         156 :       if (present(obliqr_in)) obliqr = obliqr_in
      79         156 :       if (present(lambm0_in)) lambm0 = lambm0_in
      80         156 :       if (present(mvelpp_in)) mvelpp = mvelpp_in
      81         156 :       if (present(obliq_in)) obliq = obliq_in
      82         156 :       if (present(mvelp_in)) mvelp = mvelp_in
      83         156 :       if (present(decln_in)) decln = decln_in
      84         156 :       if (present(eccf_in)) eccf = eccf_in
      85         156 :       if (present(log_print_in)) log_print = log_print_in
      86             : 
      87             :       end subroutine icepack_init_orbit
      88             :  
      89             : !=======================================================================
      90             : 
      91             : !autodocument_start icepack_query_orbit
      92             : ! Compute orbital parameters for the specified date.
      93             : 
      94           0 :       subroutine icepack_query_orbit(iyear_AD_out, eccen_out, obliqr_out, &
      95             :          lambm0_out, mvelpp_out, obliq_out, mvelp_out, decln_out, eccf_out, &
      96             :          log_print_out)
      97             : 
      98             :       integer(kind=int_kind), optional, intent(out) :: iyear_AD_out  ! Year to calculate orbit for
      99             :       real(kind=dbl_kind), optional, intent(out) :: eccen_out  ! Earth's orbital eccentricity
     100             :       real(kind=dbl_kind), optional, intent(out) :: obliqr_out ! Earth's obliquity in radians
     101             :       real(kind=dbl_kind), optional, intent(out) :: lambm0_out ! Mean longitude of perihelion at the
     102             :                                                              ! vernal equinox (radians)
     103             :       real(kind=dbl_kind), optional, intent(out) :: mvelpp_out ! Earth's moving vernal equinox longitude
     104             :                                                              ! of perihelion + pi (radians)
     105             :       real(kind=dbl_kind), optional, intent(out) :: obliq_out  ! obliquity in degrees
     106             :       real(kind=dbl_kind), optional, intent(out) :: mvelp_out  ! moving vernal equinox long
     107             :       real(kind=dbl_kind), optional, intent(out) :: decln_out  ! solar declination angle in radians
     108             :       real(kind=dbl_kind), optional, intent(out) :: eccf_out   ! earth orbit eccentricity factor
     109             :       logical(kind=log_kind), optional, intent(out) :: log_print_out ! Flags print of status/error
     110             : 
     111             : !autodocument_end
     112             : 
     113             :       character(len=*),parameter :: subname='(icepack_query_orbit)'
     114             : 
     115             :       !call icepack_warnings_add(subname//'  ')
     116           0 :       iyear_AD  = 1950
     117           0 :       log_print = .false.   ! if true, write out orbital parameters
     118             : 
     119             : #ifndef CESMCOUPLED
     120             :       call shr_orb_params( iyear_AD, eccen , obliq , mvelp    , &
     121           0 :                            obliqr  , lambm0, mvelpp, log_print)
     122           0 :       if (icepack_warnings_aborted(subname)) return
     123             : #endif
     124             : 
     125           0 :       if (present(iyear_AD_out)) iyear_AD_out = iyear_AD
     126           0 :       if (present(eccen_out)) eccen_out = eccen
     127           0 :       if (present(obliqr_out)) obliqr_out = obliqr
     128           0 :       if (present(lambm0_out)) lambm0_out = lambm0
     129           0 :       if (present(mvelpp_out)) mvelpp_out = mvelpp
     130           0 :       if (present(obliq_out)) obliq_out = obliq
     131           0 :       if (present(mvelp_out)) mvelp_out = mvelp
     132           0 :       if (present(decln_out)) decln_out = decln
     133           0 :       if (present(eccf_out)) eccf_out = eccf
     134           0 :       if (present(log_print_out)) log_print_out = log_print
     135             : 
     136             :       end subroutine icepack_query_orbit
     137             :  
     138             : !=======================================================================
     139             : 
     140             : ! Uses orbital and lat/lon info to compute cosine solar zenith angle
     141             : ! for the specified date.
     142             : !
     143             : ! author:  Bruce P. Briegleb, NCAR 
     144             : 
     145      882945 :       subroutine compute_coszen (tlat,          tlon,        &
     146             :                                  yday,  sec, coszen,     &
     147             :                                  days_per_year, nextsw_cday, &
     148           0 :                                  calendar_type)
     149             : 
     150             : #ifdef CESMCOUPLED
     151             :       use shr_orb_mod, only: shr_orb_decl
     152             : #endif
     153             :  
     154             :       real (kind=dbl_kind), intent(in) :: &
     155             :          tlat, tlon          ! latitude and longitude (radians)
     156             : 
     157             :       integer (kind=int_kind), intent(in) :: &
     158             :          sec                 ! elapsed seconds into date
     159             : 
     160             :       real (kind=dbl_kind), intent(in) :: &
     161             :          yday                ! day of the year
     162             : 
     163             :       real (kind=dbl_kind), intent(inout) :: &
     164             :          coszen              ! cosine solar zenith angle 
     165             :                              ! negative for sun below horizon
     166             :  
     167             :       integer (kind=int_kind), intent(in), optional :: &
     168             :          days_per_year       ! number of days in one year
     169             : 
     170             :       real (kind=dbl_kind), intent(in), optional :: &
     171             :          nextsw_cday         ! julian day of next shortwave calculation
     172             : 
     173             :       character (len=char_len), intent(in), optional :: &
     174             :          calendar_type       ! differentiates Gregorian from other calendars
     175             : 
     176             :       ! local variables
     177             : 
     178      303198 :       real (kind=dbl_kind) :: ydayp1 ! day of year plus one time step
     179             :  
     180             :       character(len=*),parameter :: subname='(compute_coszen)'
     181             : 
     182             : ! Solar declination for next time step
     183             :  
     184             : #ifdef CESMCOUPLED
     185             :       if (calendar_type == "GREGORIAN") then
     186             :          ydayp1 = min(nextsw_cday, real(days_per_year,kind=dbl_kind))
     187             :       else
     188             :          ydayp1 = nextsw_cday
     189             :       endif
     190             : 
     191             :       !--- update coszen when nextsw_cday valid
     192             :       if (ydayp1 > -0.5_dbl_kind) then
     193             : #else
     194      882945 :       ydayp1 = yday + sec/secday
     195             : #endif
     196             :  
     197             :       call shr_orb_decl(ydayp1, eccen, mvelpp, lambm0, &
     198      882945 :                         obliqr, decln, eccf)
     199      882945 :       if (icepack_warnings_aborted(subname)) return
     200             : 
     201             :       coszen = sin(tlat)*sin(decln) &
     202             :              + cos(tlat)*cos(decln) &
     203      882945 :              *cos((sec/secday-p5)*c2*pi + tlon) !cos(hour angle)
     204             :  
     205             : #ifdef CESMCOUPLED
     206             :       endif
     207             : #endif
     208             : 
     209             :       end subroutine compute_coszen
     210             :  
     211             : !===============================================================================
     212             : 
     213             : #ifndef CESMCOUPLED
     214         156 : SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp    , &
     215             :            &               obliqr   , lambm0, mvelpp, log_print)
     216             : 
     217             : !-------------------------------------------------------------------------------
     218             : !
     219             : ! Calculate earths orbital parameters using Dave Threshers formula which 
     220             : ! came from Berger, Andre.  1978  A Simple Algorithm to Compute Long-Term 
     221             : ! Variations of Daily Insolation.  Contribution 18, Institute of Astronomy 
     222             : ! and Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium
     223             : !
     224             : !------------------------------Code history-------------------------------------
     225             : !
     226             : ! Original Author: Erik Kluzek
     227             : ! Date:            Oct/97
     228             : !
     229             : !-------------------------------------------------------------------------------
     230             : 
     231             :    !----------------------------- Arguments ------------------------------------
     232             :    integer(int_kind),intent(in)    :: iyear_AD  ! Year to calculate orbit for
     233             :    real   (dbl_kind),intent(inout) :: eccen     ! orbital eccentricity
     234             :    real   (dbl_kind),intent(inout) :: obliq     ! obliquity in degrees
     235             :    real   (dbl_kind),intent(inout) :: mvelp     ! moving vernal equinox long
     236             :    real   (dbl_kind),intent(out)   :: obliqr    ! Earths obliquity in rad
     237             :    real   (dbl_kind),intent(out)   :: lambm0    ! Mean long of perihelion at
     238             :                                                    ! vernal equinox (radians)
     239             :    real   (dbl_kind),intent(out)   :: mvelpp    ! moving vernal equinox long
     240             :                                                    ! of perihelion plus pi (rad)
     241             :    logical(log_kind),intent(in)    :: log_print ! Flags print of status/error
     242             : 
     243             :    !------------------------------ Parameters ----------------------------------
     244             :    real   (dbl_kind),parameter :: SHR_ORB_UNDEF_REAL = 1.e36_dbl_kind ! undefined real 
     245             :    integer(int_kind),parameter :: SHR_ORB_UNDEF_INT  = 2000000000        ! undefined int
     246             :    integer(int_kind),parameter :: poblen =47 ! # of elements in series wrt obliquity
     247             :    integer(int_kind),parameter :: pecclen=19 ! # of elements in series wrt eccentricity
     248             :    integer(int_kind),parameter :: pmvelen=78 ! # of elements in series wrt vernal equinox
     249             :    real   (dbl_kind),parameter :: psecdeg = 1.0_dbl_kind/3600.0_dbl_kind ! arc sec to deg conversion
     250             : 
     251          56 :    real   (dbl_kind) :: yb4_1950AD         ! number of years before 1950 AD
     252             : 
     253             :    real   (dbl_kind),parameter :: SHR_ORB_ECCEN_MIN  =   0.0_dbl_kind ! min value for eccen
     254             :    real   (dbl_kind),parameter :: SHR_ORB_ECCEN_MAX  =   0.1_dbl_kind ! max value for eccen
     255             :    real   (dbl_kind),parameter :: SHR_ORB_OBLIQ_MIN  = -90.0_dbl_kind ! min value for obliq
     256             :    real   (dbl_kind),parameter :: SHR_ORB_OBLIQ_MAX  = +90.0_dbl_kind ! max value for obliq
     257             :    real   (dbl_kind),parameter :: SHR_ORB_MVELP_MIN  =   0.0_dbl_kind ! min value for mvelp
     258             :    real   (dbl_kind),parameter :: SHR_ORB_MVELP_MAX  = 360.0_dbl_kind ! max value for mvelp
     259             : 
     260             :    ! Cosine series data for computation of obliquity: amplitude (arc seconds),
     261             :    ! rate (arc seconds/year), phase (degrees).
     262             :  
     263             :    real   (dbl_kind), parameter :: obamp(poblen) =  & ! amplitudes for obliquity cos series
     264             :    &      (/   -2462.2214466_dbl_kind, -857.3232075_dbl_kind, -629.3231835_dbl_kind,   &
     265             :    &            -414.2804924_dbl_kind, -311.7632587_dbl_kind,  308.9408604_dbl_kind,   &
     266             :    &            -162.5533601_dbl_kind, -116.1077911_dbl_kind,  101.1189923_dbl_kind,   &
     267             :    &             -67.6856209_dbl_kind,   24.9079067_dbl_kind,   22.5811241_dbl_kind,   &
     268             :    &             -21.1648355_dbl_kind,  -15.6549876_dbl_kind,   15.3936813_dbl_kind,   &
     269             :    &              14.6660938_dbl_kind,  -11.7273029_dbl_kind,   10.2742696_dbl_kind,   &
     270             :    &               6.4914588_dbl_kind,    5.8539148_dbl_kind,   -5.4872205_dbl_kind,   &
     271             :    &              -5.4290191_dbl_kind,    5.1609570_dbl_kind,    5.0786314_dbl_kind,   &
     272             :    &              -4.0735782_dbl_kind,    3.7227167_dbl_kind,    3.3971932_dbl_kind,   &
     273             :    &              -2.8347004_dbl_kind,   -2.6550721_dbl_kind,   -2.5717867_dbl_kind,   &
     274             :    &              -2.4712188_dbl_kind,    2.4625410_dbl_kind,    2.2464112_dbl_kind,   &
     275             :    &              -2.0755511_dbl_kind,   -1.9713669_dbl_kind,   -1.8813061_dbl_kind,   &
     276             :    &              -1.8468785_dbl_kind,    1.8186742_dbl_kind,    1.7601888_dbl_kind,   &
     277             :    &              -1.5428851_dbl_kind,    1.4738838_dbl_kind,   -1.4593669_dbl_kind,   &
     278             :    &               1.4192259_dbl_kind,   -1.1818980_dbl_kind,    1.1756474_dbl_kind,   &
     279             :    &              -1.1316126_dbl_kind,    1.0896928_dbl_kind/)
     280             :  
     281             :    real   (dbl_kind), parameter :: obrate(poblen) = & ! rates for obliquity cosine series
     282             :    &        (/  31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind,   &
     283             :    &            31.983787_dbl_kind, 44.828336_dbl_kind, 30.973257_dbl_kind,   &
     284             :    &            43.668246_dbl_kind, 32.246691_dbl_kind, 30.599444_dbl_kind,   &
     285             :    &            42.681324_dbl_kind, 43.836462_dbl_kind, 47.439436_dbl_kind,   &
     286             :    &            63.219948_dbl_kind, 64.230478_dbl_kind,  1.010530_dbl_kind,   &
     287             :    &             7.437771_dbl_kind, 55.782177_dbl_kind,  0.373813_dbl_kind,   &
     288             :    &            13.218362_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind,   &
     289             :    &            76.438310_dbl_kind, 45.815258_dbl_kind,  8.448301_dbl_kind,   &
     290             :    &            56.792707_dbl_kind, 49.747842_dbl_kind, 12.058272_dbl_kind,   &
     291             :    &            75.278220_dbl_kind, 65.241008_dbl_kind, 64.604291_dbl_kind,   &
     292             :    &             1.647247_dbl_kind,  7.811584_dbl_kind, 12.207832_dbl_kind,   &
     293             :    &            63.856665_dbl_kind, 56.155990_dbl_kind, 77.448840_dbl_kind,   &
     294             :    &             6.801054_dbl_kind, 62.209418_dbl_kind, 20.656133_dbl_kind,   &
     295             :    &            48.344406_dbl_kind, 55.145460_dbl_kind, 69.000539_dbl_kind,   &
     296             :    &            11.071350_dbl_kind, 74.291298_dbl_kind, 11.047742_dbl_kind,   &
     297             :    &             0.636717_dbl_kind, 12.844549_dbl_kind/)
     298             :  
     299             :    real   (dbl_kind), parameter :: obphas(poblen) = & ! phases for obliquity cosine series
     300             :    &      (/    251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind,   &
     301             :    &            292.7252_dbl_kind,  15.3747_dbl_kind, 263.7951_dbl_kind,   &
     302             :    &            308.4258_dbl_kind, 240.0099_dbl_kind, 222.9725_dbl_kind,   &
     303             :    &            268.7809_dbl_kind, 316.7998_dbl_kind, 319.6024_dbl_kind,   &
     304             :    &            143.8050_dbl_kind, 172.7351_dbl_kind,  28.9300_dbl_kind,   &
     305             :    &            123.5968_dbl_kind,  20.2082_dbl_kind,  40.8226_dbl_kind,   &
     306             :    &            123.4722_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind,   &
     307             :    &            267.2772_dbl_kind,  55.0196_dbl_kind, 152.5268_dbl_kind,   &
     308             :    &             49.1382_dbl_kind, 204.6609_dbl_kind,  56.5233_dbl_kind,   &
     309             :    &            200.3284_dbl_kind, 201.6651_dbl_kind, 213.5577_dbl_kind,   &
     310             :    &             17.0374_dbl_kind, 164.4194_dbl_kind,  94.5422_dbl_kind,   &
     311             :    &            131.9124_dbl_kind,  61.0309_dbl_kind, 296.2073_dbl_kind,   &
     312             :    &            135.4894_dbl_kind, 114.8750_dbl_kind, 247.0691_dbl_kind,   &
     313             :    &            256.6114_dbl_kind,  32.1008_dbl_kind, 143.6804_dbl_kind,   &
     314             :    &             16.8784_dbl_kind, 160.6835_dbl_kind,  27.5932_dbl_kind,   &
     315             :    &            348.1074_dbl_kind,  82.6496_dbl_kind/)
     316             :  
     317             :    ! Cosine/sine series data for computation of eccentricity and fixed vernal 
     318             :    ! equinox longitude of perihelion (fvelp): amplitude, 
     319             :    ! rate (arc seconds/year), phase (degrees).
     320             :  
     321             :    real   (dbl_kind), parameter :: ecamp (pecclen) = & ! ampl for eccen/fvelp cos/sin series
     322             :    &      (/   0.01860798_dbl_kind,  0.01627522_dbl_kind, -0.01300660_dbl_kind,   &
     323             :    &           0.00988829_dbl_kind, -0.00336700_dbl_kind,  0.00333077_dbl_kind,   &
     324             :    &          -0.00235400_dbl_kind,  0.00140015_dbl_kind,  0.00100700_dbl_kind,   &
     325             :    &           0.00085700_dbl_kind,  0.00064990_dbl_kind,  0.00059900_dbl_kind,   &
     326             :    &           0.00037800_dbl_kind, -0.00033700_dbl_kind,  0.00027600_dbl_kind,   &
     327             :    &           0.00018200_dbl_kind, -0.00017400_dbl_kind, -0.00012400_dbl_kind,   &
     328             :    &           0.00001250_dbl_kind/)
     329             :  
     330             :    real   (dbl_kind), parameter :: ecrate(pecclen) = & ! rates for eccen/fvelp cos/sin series
     331             :    &      (/    4.2072050_dbl_kind,  7.3460910_dbl_kind, 17.8572630_dbl_kind,  &
     332             :    &           17.2205460_dbl_kind, 16.8467330_dbl_kind,  5.1990790_dbl_kind,  &
     333             :    &           18.2310760_dbl_kind, 26.2167580_dbl_kind,  6.3591690_dbl_kind,  &
     334             :    &           16.2100160_dbl_kind,  3.0651810_dbl_kind, 16.5838290_dbl_kind,  &
     335             :    &           18.4939800_dbl_kind,  6.1909530_dbl_kind, 18.8677930_dbl_kind,  &
     336             :    &           17.4255670_dbl_kind,  6.1860010_dbl_kind, 18.4174410_dbl_kind,  &
     337             :    &            0.6678630_dbl_kind/)
     338             :  
     339             :    real   (dbl_kind), parameter :: ecphas(pecclen) = & ! phases for eccen/fvelp cos/sin series
     340             :    &      (/    28.620089_dbl_kind, 193.788772_dbl_kind, 308.307024_dbl_kind,  &
     341             :    &           320.199637_dbl_kind, 279.376984_dbl_kind,  87.195000_dbl_kind,  &
     342             :    &           349.129677_dbl_kind, 128.443387_dbl_kind, 154.143880_dbl_kind,  &
     343             :    &           291.269597_dbl_kind, 114.860583_dbl_kind, 332.092251_dbl_kind,  &
     344             :    &           296.414411_dbl_kind, 145.769910_dbl_kind, 337.237063_dbl_kind,  &
     345             :    &           152.092288_dbl_kind, 126.839891_dbl_kind, 210.667199_dbl_kind,  &
     346             :    &            72.108838_dbl_kind/)
     347             :  
     348             :    ! Sine series data for computation of moving vernal equinox longitude of 
     349             :    ! perihelion: amplitude (arc seconds), rate (arc sec/year), phase (degrees).      
     350             :  
     351             :    real   (dbl_kind), parameter :: mvamp (pmvelen) = & ! amplitudes for mvelp sine series 
     352             :    &      (/   7391.0225890_dbl_kind, 2555.1526947_dbl_kind, 2022.7629188_dbl_kind,  &
     353             :    &          -1973.6517951_dbl_kind, 1240.2321818_dbl_kind,  953.8679112_dbl_kind,  &
     354             :    &           -931.7537108_dbl_kind,  872.3795383_dbl_kind,  606.3544732_dbl_kind,  &
     355             :    &           -496.0274038_dbl_kind,  456.9608039_dbl_kind,  346.9462320_dbl_kind,  &
     356             :    &           -305.8412902_dbl_kind,  249.6173246_dbl_kind, -199.1027200_dbl_kind,  &
     357             :    &            191.0560889_dbl_kind, -175.2936572_dbl_kind,  165.9068833_dbl_kind,  &
     358             :    &            161.1285917_dbl_kind,  139.7878093_dbl_kind, -133.5228399_dbl_kind,  &
     359             :    &            117.0673811_dbl_kind,  104.6907281_dbl_kind,   95.3227476_dbl_kind,  &
     360             :    &             86.7824524_dbl_kind,   86.0857729_dbl_kind,   70.5893698_dbl_kind,  &
     361             :    &            -69.9719343_dbl_kind,  -62.5817473_dbl_kind,   61.5450059_dbl_kind,  &
     362             :    &            -57.9364011_dbl_kind,   57.1899832_dbl_kind,  -57.0236109_dbl_kind,  &
     363             :    &            -54.2119253_dbl_kind,   53.2834147_dbl_kind,   52.1223575_dbl_kind,  &
     364             :    &            -49.0059908_dbl_kind,  -48.3118757_dbl_kind,  -45.4191685_dbl_kind,  &
     365             :    &            -42.2357920_dbl_kind,  -34.7971099_dbl_kind,   34.4623613_dbl_kind,  &
     366             :    &            -33.8356643_dbl_kind,   33.6689362_dbl_kind,  -31.2521586_dbl_kind,  &
     367             :    &            -30.8798701_dbl_kind,   28.4640769_dbl_kind,  -27.1960802_dbl_kind,  &
     368             :    &             27.0860736_dbl_kind,  -26.3437456_dbl_kind,   24.7253740_dbl_kind,  &
     369             :    &             24.6732126_dbl_kind,   24.4272733_dbl_kind,   24.0127327_dbl_kind,  &
     370             :    &             21.7150294_dbl_kind,  -21.5375347_dbl_kind,   18.1148363_dbl_kind,  &
     371             :    &            -16.9603104_dbl_kind,  -16.1765215_dbl_kind,   15.5567653_dbl_kind,  &
     372             :    &             15.4846529_dbl_kind,   15.2150632_dbl_kind,   14.5047426_dbl_kind,  &
     373             :    &            -14.3873316_dbl_kind,   13.1351419_dbl_kind,   12.8776311_dbl_kind,  &
     374             :    &             11.9867234_dbl_kind,   11.9385578_dbl_kind,   11.7030822_dbl_kind,  &
     375             :    &             11.6018181_dbl_kind,  -11.2617293_dbl_kind,  -10.4664199_dbl_kind,  &
     376             :    &             10.4333970_dbl_kind,  -10.2377466_dbl_kind,   10.1934446_dbl_kind,  &
     377             :    &            -10.1280191_dbl_kind,   10.0289441_dbl_kind,  -10.0034259_dbl_kind/)
     378             :  
     379             :    real   (dbl_kind), parameter :: mvrate(pmvelen) = & ! rates for mvelp sine series 
     380             :    &      (/    31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind,   &
     381             :    &             0.636717_dbl_kind, 31.983787_dbl_kind,  3.138886_dbl_kind,   &
     382             :    &            30.973257_dbl_kind, 44.828336_dbl_kind,  0.991874_dbl_kind,   &
     383             :    &             0.373813_dbl_kind, 43.668246_dbl_kind, 32.246691_dbl_kind,   &
     384             :    &            30.599444_dbl_kind,  2.147012_dbl_kind, 10.511172_dbl_kind,   &
     385             :    &            42.681324_dbl_kind, 13.650058_dbl_kind,  0.986922_dbl_kind,   &
     386             :    &             9.874455_dbl_kind, 13.013341_dbl_kind,  0.262904_dbl_kind,   &
     387             :    &             0.004952_dbl_kind,  1.142024_dbl_kind, 63.219948_dbl_kind,   &
     388             :    &             0.205021_dbl_kind,  2.151964_dbl_kind, 64.230478_dbl_kind,   &
     389             :    &            43.836462_dbl_kind, 47.439436_dbl_kind,  1.384343_dbl_kind,   &
     390             :    &             7.437771_dbl_kind, 18.829299_dbl_kind,  9.500642_dbl_kind,   &
     391             :    &             0.431696_dbl_kind,  1.160090_dbl_kind, 55.782177_dbl_kind,   &
     392             :    &            12.639528_dbl_kind,  1.155138_dbl_kind,  0.168216_dbl_kind,   &
     393             :    &             1.647247_dbl_kind, 10.884985_dbl_kind,  5.610937_dbl_kind,   &
     394             :    &            12.658184_dbl_kind,  1.010530_dbl_kind,  1.983748_dbl_kind,   &
     395             :    &            14.023871_dbl_kind,  0.560178_dbl_kind,  1.273434_dbl_kind,   &
     396             :    &            12.021467_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind,   &
     397             :    &            76.438310_dbl_kind,  4.280910_dbl_kind, 13.218362_dbl_kind,   &
     398             :    &            17.818769_dbl_kind,  8.359495_dbl_kind, 56.792707_dbl_kind,   &
     399             :    &            8.448301_dbl_kind,  1.978796_dbl_kind,  8.863925_dbl_kind,   &
     400             :    &             0.186365_dbl_kind,  8.996212_dbl_kind,  6.771027_dbl_kind,   &
     401             :    &            45.815258_dbl_kind, 12.002811_dbl_kind, 75.278220_dbl_kind,   &
     402             :    &            65.241008_dbl_kind, 18.870667_dbl_kind, 22.009553_dbl_kind,   &
     403             :    &            64.604291_dbl_kind, 11.498094_dbl_kind,  0.578834_dbl_kind,   &
     404             :    &             9.237738_dbl_kind, 49.747842_dbl_kind,  2.147012_dbl_kind,   &
     405             :    &             1.196895_dbl_kind,  2.133898_dbl_kind,  0.173168_dbl_kind/)
     406             : 
     407             :    real   (dbl_kind), parameter :: mvphas(pmvelen) = & ! phases for mvelp sine series
     408             :    &      (/    251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind,   &
     409             :    &            348.1074_dbl_kind, 292.7252_dbl_kind, 165.1686_dbl_kind,   &
     410             :    &            263.7951_dbl_kind,  15.3747_dbl_kind,  58.5749_dbl_kind,   &
     411             :    &             40.8226_dbl_kind, 308.4258_dbl_kind, 240.0099_dbl_kind,   &
     412             :    &            222.9725_dbl_kind, 106.5937_dbl_kind, 114.5182_dbl_kind,   &
     413             :    &            268.7809_dbl_kind, 279.6869_dbl_kind,  39.6448_dbl_kind,   &
     414             :    &            126.4108_dbl_kind, 291.5795_dbl_kind, 307.2848_dbl_kind,   &
     415             :    &             18.9300_dbl_kind, 273.7596_dbl_kind, 143.8050_dbl_kind,   &
     416             :    &            191.8927_dbl_kind, 125.5237_dbl_kind, 172.7351_dbl_kind,   &
     417             :    &            316.7998_dbl_kind, 319.6024_dbl_kind,  69.7526_dbl_kind,   &
     418             :    &            123.5968_dbl_kind, 217.6432_dbl_kind,  85.5882_dbl_kind,   &
     419             :    &            156.2147_dbl_kind,  66.9489_dbl_kind,  20.2082_dbl_kind,   &
     420             :    &            250.7568_dbl_kind,  48.0188_dbl_kind,   8.3739_dbl_kind,   &
     421             :    &             17.0374_dbl_kind, 155.3409_dbl_kind,  94.1709_dbl_kind,   &
     422             :    &            221.1120_dbl_kind,  28.9300_dbl_kind, 117.1498_dbl_kind,   &
     423             :    &            320.5095_dbl_kind, 262.3602_dbl_kind, 336.2148_dbl_kind,   &
     424             :    &            233.0046_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind,   &
     425             :    &            267.2772_dbl_kind,  78.9281_dbl_kind, 123.4722_dbl_kind,   &
     426             :    &            188.7132_dbl_kind, 180.1364_dbl_kind,  49.1382_dbl_kind,   &
     427             :    &            152.5268_dbl_kind,  98.2198_dbl_kind,  97.4808_dbl_kind,   &
     428             :    &            221.5376_dbl_kind, 168.2438_dbl_kind, 161.1199_dbl_kind,   &
     429             :    &             55.0196_dbl_kind, 262.6495_dbl_kind, 200.3284_dbl_kind,   &
     430             :    &            201.6651_dbl_kind, 294.6547_dbl_kind,  99.8233_dbl_kind,   &
     431             :    &            213.5577_dbl_kind, 154.1631_dbl_kind, 232.7153_dbl_kind,   &
     432             :    &            138.3034_dbl_kind, 204.6609_dbl_kind, 106.5938_dbl_kind,   &
     433             :    &            250.4676_dbl_kind, 332.3345_dbl_kind,  27.3039_dbl_kind/)
     434             :  
     435             :    !---------------------------Local variables----------------------------------
     436             :    integer(int_kind) :: i       ! Index for series summations
     437          56 :    real   (dbl_kind) :: obsum   ! Obliquity series summation
     438          56 :    real   (dbl_kind) :: cossum  ! Cos series summation for eccentricity/fvelp
     439          56 :    real   (dbl_kind) :: sinsum  ! Sin series summation for eccentricity/fvelp
     440          56 :    real   (dbl_kind) :: fvelp   ! Fixed vernal equinox long of perihelion
     441          56 :    real   (dbl_kind) :: mvsum   ! mvelp series summation
     442          56 :    real   (dbl_kind) :: beta    ! Intermediate argument for lambm0
     443          56 :    real   (dbl_kind) :: years   ! Years to time of interest ( pos <=> future)
     444          56 :    real   (dbl_kind) :: eccen2  ! eccentricity squared
     445          56 :    real   (dbl_kind) :: eccen3  ! eccentricity cubed
     446          56 :    real   (dbl_kind) :: degrad  ! degrees to rad conversion
     447             :    integer (int_kind), parameter :: s_loglev    = 0         
     448             :    character(len=*),parameter :: subname='(shr_orb_params)'
     449             :  
     450             :    !-------------------------- Formats -----------------------------------------
     451             :    character(len=*),parameter :: F00 = "('(shr_orb_params) ',4a)"
     452             :    character(len=*),parameter :: F01 = "('(shr_orb_params) ',a,i9)"
     453             :    character(len=*),parameter :: F02 = "('(shr_orb_params) ',a,f6.3)"
     454             :    character(len=*),parameter :: F03 = "('(shr_orb_params) ',a,es14.6)"
     455             : 
     456             :    !----------------------------------------------------------------------------
     457             :    ! radinp and algorithms below will need a degree to radian conversion factor
     458             : 
     459             :    !call icepack_warnings_add(subname//'  ')
     460         156 :    degrad = pi/180._dbl_kind   ! degree to radian conversion factor
     461             :  
     462             :    if ( log_print .and. s_loglev > 0 ) then
     463             :      write(warnstr,F00) subname//'Calculate characteristics of the orbit:'
     464             :      call icepack_warnings_add(warnstr)
     465             :    end if
     466             :  
     467             :    ! Check for flag to use input orbit parameters
     468             :  
     469         156 :    IF ( iyear_AD == SHR_ORB_UNDEF_INT ) THEN
     470             : 
     471             :       ! Check input obliq, eccen, and mvelp to ensure reasonable
     472             :  
     473           0 :       if( obliq == SHR_ORB_UNDEF_REAL )then
     474           0 :          write(warnstr,F00) subname//' Have to specify orbital parameters:'
     475           0 :          call icepack_warnings_add(warnstr)
     476           0 :          write(warnstr,F00) subname//'Either set: iyear_AD, OR [obliq, eccen, and mvelp]:'
     477           0 :          call icepack_warnings_add(warnstr)
     478           0 :          write(warnstr,F00) subname//'iyear_AD is the year to simulate orbit for (ie. 1950): '
     479           0 :          call icepack_warnings_add(warnstr)
     480           0 :          write(warnstr,F00) subname//'obliq, eccen, mvelp specify the orbit directly:'
     481           0 :          call icepack_warnings_add(warnstr)
     482           0 :          write(warnstr,F00) subname//'The AMIP II settings (for a 1995 orbit) are: '
     483           0 :          call icepack_warnings_add(warnstr)
     484           0 :          write(warnstr,F00) subname//' obliq =  23.4441'
     485           0 :          call icepack_warnings_add(warnstr)
     486           0 :          write(warnstr,F00) subname//' eccen =   0.016715'
     487           0 :          call icepack_warnings_add(warnstr)
     488           0 :          write(warnstr,F00) subname//' mvelp = 102.7'
     489           0 :          call icepack_warnings_add(warnstr)
     490           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     491           0 :          call icepack_warnings_add(subname//' unreasonable oblip')
     492           0 :       else if ( log_print ) then
     493           0 :          write(warnstr,F00) subname//'Use input orbital parameters: '
     494           0 :          call icepack_warnings_add(warnstr)
     495             :       end if
     496           0 :       if( (obliq < SHR_ORB_OBLIQ_MIN).or.(obliq > SHR_ORB_OBLIQ_MAX) ) then
     497           0 :          write(warnstr,F03) subname//'Input obliquity unreasonable: ', obliq
     498           0 :          call icepack_warnings_add(warnstr)
     499           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     500           0 :          call icepack_warnings_add(subname//' unreasonable obliq')
     501             :       end if
     502           0 :       if( (eccen < SHR_ORB_ECCEN_MIN).or.(eccen > SHR_ORB_ECCEN_MAX) ) then
     503           0 :          write(warnstr,F03) subname//'Input eccentricity unreasonable: ', eccen
     504           0 :          call icepack_warnings_add(warnstr)
     505           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     506           0 :          call icepack_warnings_add(subname//' unreasonable eccen')
     507             :       end if
     508           0 :       if( (mvelp < SHR_ORB_MVELP_MIN).or.(mvelp > SHR_ORB_MVELP_MAX) ) then
     509           0 :          write(warnstr,F03) subname//'Input mvelp unreasonable: ' , mvelp
     510           0 :          call icepack_warnings_add(warnstr)
     511           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     512           0 :          call icepack_warnings_add(subname//' unreasonable mvelp')
     513             :       end if
     514           0 :       eccen2 = eccen*eccen
     515           0 :       eccen3 = eccen2*eccen
     516             : 
     517             :    ELSE  ! Otherwise calculate based on years before present
     518             :  
     519             :       if ( log_print .and. s_loglev > 0) then
     520             :          write(warnstr,F01) subname//'Calculate orbit for year: ' , iyear_AD
     521             :          call icepack_warnings_add(warnstr)
     522             :       end if
     523         156 :       yb4_1950AD = 1950.0_dbl_kind - real(iyear_AD,dbl_kind)
     524         156 :       if ( abs(yb4_1950AD) .gt. 1000000.0_dbl_kind )then
     525           0 :          write(warnstr,F00) subname//'orbit only valid for years+-1000000'
     526           0 :          call icepack_warnings_add(warnstr)
     527           0 :          write(warnstr,F00) subname//'Relative to 1950 AD'
     528           0 :          call icepack_warnings_add(warnstr)
     529           0 :          write(warnstr,F03) subname//'# of years before 1950: ',yb4_1950AD
     530           0 :          call icepack_warnings_add(warnstr)
     531           0 :          write(warnstr,F01) subname//'Year to simulate was  : ',iyear_AD
     532           0 :          call icepack_warnings_add(warnstr)
     533           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     534           0 :          call icepack_warnings_add(subname//' unreasonable year')
     535             :       end if
     536             :  
     537             :       ! The following calculates the earths obliquity, orbital eccentricity
     538             :       ! (and various powers of it) and vernal equinox mean longitude of
     539             :       ! perihelion for years in the past (future = negative of years past),
     540             :       ! using constants (see parameter section) given in the program of:
     541             :       !
     542             :       ! Berger, Andre.  1978  A Simple Algorithm to Compute Long-Term Variations
     543             :       ! of Daily Insolation.  Contribution 18, Institute of Astronomy and
     544             :       ! Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium.
     545             :       !
     546             :       ! and formulas given in the paper (where less precise constants are also
     547             :       ! given):
     548             :       !
     549             :       ! Berger, Andre.  1978.  Long-Term Variations of Daily Insolation and
     550             :       ! Quaternary Climatic Changes.  J. of the Atmo. Sci. 35:2362-2367
     551             :       !
     552             :       ! The algorithm is valid only to 1,000,000 years past or hence.
     553             :       ! For a solution valid to 5-10 million years past see the above author.
     554             :       ! Algorithm below is better for years closer to present than is the
     555             :       ! 5-10 million year solution.
     556             :       !
     557             :       ! Years to time of interest must be negative of years before present
     558             :       ! (1950) in formulas that follow. 
     559             :  
     560         156 :       years = - yb4_1950AD
     561             :  
     562             :       ! In the summations below, cosine or sine arguments, which end up in
     563             :       ! degrees, must be converted to radians via multiplication by degrad.
     564             :       !
     565             :       ! Summation of cosine series for obliquity (epsilon in Berger 1978) in
     566             :       ! degrees. Convert the amplitudes and rates, which are in arc secs, into
     567             :       ! degrees via multiplication by psecdeg (arc seconds to degrees conversion
     568             :       ! factor).  For obliq, first term is Berger 1978 epsilon star; second
     569             :       ! term is series summation in degrees.
     570             :   
     571         156 :       obsum = 0.0_dbl_kind
     572        7488 :       do i = 1, poblen
     573        2632 :          obsum = obsum + obamp(i)*psecdeg*cos((obrate(i)*psecdeg*years + &
     574        7488 :          &       obphas(i))*degrad)
     575             :       end do
     576         156 :       obliq = 23.320556_dbl_kind + obsum
     577             :  
     578             :       ! Summation of cosine and sine series for computation of eccentricity 
     579             :       ! (eccen; e in Berger 1978) and fixed vernal equinox longitude of 
     580             :       ! perihelion (fvelp; pi in Berger 1978), which is used for computation 
     581             :       ! of moving vernal equinox longitude of perihelion.  Convert the rates, 
     582             :       ! which are in arc seconds, into degrees via multiplication by psecdeg.
     583             :  
     584         156 :       cossum = 0.0_dbl_kind
     585        3120 :       do i = 1, pecclen
     586        3120 :         cossum = cossum+ecamp(i)*cos((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
     587             :       end do
     588             :  
     589         156 :       sinsum = 0.0_dbl_kind
     590        3120 :       do i = 1, pecclen
     591        3120 :         sinsum = sinsum+ecamp(i)*sin((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
     592             :       end do
     593             :  
     594             :       ! Use summations to calculate eccentricity
     595             :  
     596         156 :       eccen2 = cossum*cossum + sinsum*sinsum
     597         156 :       eccen  = sqrt(eccen2)
     598         156 :       eccen3 = eccen2*eccen
     599             :  
     600             :       ! A series of cases for fvelp, which is in radians.
     601             :          
     602         156 :       if (abs(cossum) .le. 1.0E-8_dbl_kind) then
     603           0 :         if (sinsum .eq. 0.0_dbl_kind) then
     604           0 :           fvelp = 0.0_dbl_kind
     605           0 :         else if (sinsum .lt. 0.0_dbl_kind) then
     606           0 :           fvelp = 1.5_dbl_kind*pi
     607           0 :         else if (sinsum .gt. 0.0_dbl_kind) then
     608           0 :           fvelp = .5_dbl_kind*pi
     609             :         endif
     610         156 :       else if (cossum .lt. 0.0_dbl_kind) then
     611         156 :         fvelp = atan(sinsum/cossum) + pi
     612           0 :       else if (cossum .gt. 0.0_dbl_kind) then
     613           0 :         if (sinsum .lt. 0.0_dbl_kind) then
     614           0 :           fvelp = atan(sinsum/cossum) + 2.0_dbl_kind*pi
     615             :         else
     616           0 :           fvelp = atan(sinsum/cossum)
     617             :         endif
     618             :       endif
     619             :  
     620             :       ! Summation of sin series for computation of moving vernal equinox long
     621             :       ! of perihelion (mvelp; omega bar in Berger 1978) in degrees.  For mvelp,
     622             :       ! first term is fvelp in degrees; second term is Berger 1978 psi bar 
     623             :       ! times years and in degrees; third term is Berger 1978 zeta; fourth 
     624             :       ! term is series summation in degrees.  Convert the amplitudes and rates,
     625             :       ! which are in arc seconds, into degrees via multiplication by psecdeg.  
     626             :       ! Series summation plus second and third terms constitute Berger 1978
     627             :       ! psi, which is the general precession.
     628             :  
     629         156 :       mvsum = 0.0_dbl_kind
     630       12324 :       do i = 1, pmvelen
     631        4368 :         mvsum = mvsum + mvamp(i)*psecdeg*sin((mvrate(i)*psecdeg*years + &
     632       12324 :         &       mvphas(i))*degrad)
     633             :       end do
     634         156 :       mvelp = fvelp/degrad + 50.439273_dbl_kind*psecdeg*years + 3.392506_dbl_kind + mvsum
     635             :  
     636             :       ! Cases to make sure mvelp is between 0 and 360.
     637             :  
     638         156 :       do while (mvelp .lt. 0.0_dbl_kind)
     639         156 :         mvelp = mvelp + 360.0_dbl_kind
     640             :       end do
     641         156 :       do while (mvelp .ge. 360.0_dbl_kind)
     642         156 :         mvelp = mvelp - 360.0_dbl_kind
     643             :       end do
     644             : 
     645             :    END IF  ! end of test on whether to calculate or use input orbital params
     646             :  
     647             :    ! Orbit needs the obliquity in radians
     648             :  
     649         156 :    obliqr = obliq*degrad
     650             :  
     651             :    ! 180 degrees must be added to mvelp since observations are made from the
     652             :    ! earth and the sun is considered (wrongly for the algorithm) to go around
     653             :    ! the earth. For a more graphic explanation see Appendix B in:
     654             :    !
     655             :    ! A. Berger, M. Loutre and C. Tricot. 1993.  Insolation and Earth Orbital
     656             :    ! Periods.  J. of Geophysical Research 98:10,341-10,362.
     657             :    !
     658             :    ! Additionally, orbit will need this value in radians. So mvelp becomes
     659             :    ! mvelpp (mvelp plus pi)
     660             :  
     661         156 :    mvelpp = (mvelp + 180._dbl_kind)*degrad
     662             :  
     663             :    ! Set up an argument used several times in lambm0 calculation ahead.
     664             :  
     665         156 :    beta = sqrt(1._dbl_kind - eccen2)
     666             :  
     667             :    ! The mean longitude at the vernal equinox (lambda m nought in Berger
     668             :    ! 1978; in radians) is calculated from the following formula given in 
     669             :    ! Berger 1978.  At the vernal equinox the true longitude (lambda in Berger
     670             :    ! 1978) is 0.
     671             : 
     672             :    lambm0 = 2._dbl_kind*((.5_dbl_kind*eccen + .125_dbl_kind*eccen3)*(1._dbl_kind + beta)*sin(mvelpp)  &
     673             :    &      - .250_dbl_kind*eccen2*(.5_dbl_kind    + beta)*sin(2._dbl_kind*mvelpp)            &
     674         156 :    &      + .125_dbl_kind*eccen3*(1._dbl_kind/3._dbl_kind + beta)*sin(3._dbl_kind*mvelpp))
     675             :  
     676         156 :    if ( log_print ) then
     677           0 :      write(warnstr,F03) subname//'------ Computed Orbital Parameters ------'
     678           0 :      call icepack_warnings_add(warnstr)
     679           0 :      write(warnstr,F03) subname//'Eccentricity      = ',eccen
     680           0 :      call icepack_warnings_add(warnstr)
     681           0 :      write(warnstr,F03) subname//'Obliquity (deg)   = ',obliq
     682           0 :      call icepack_warnings_add(warnstr)
     683           0 :      write(warnstr,F03) subname//'Obliquity (rad)   = ',obliqr
     684           0 :      call icepack_warnings_add(warnstr)
     685           0 :      write(warnstr,F03) subname//'Long of perh(deg) = ',mvelp
     686           0 :      call icepack_warnings_add(warnstr)
     687           0 :      write(warnstr,F03) subname//'Long of perh(rad) = ',mvelpp
     688           0 :      call icepack_warnings_add(warnstr)
     689           0 :      write(warnstr,F03) subname//'Long at v.e.(rad) = ',lambm0
     690           0 :      call icepack_warnings_add(warnstr)
     691           0 :      write(warnstr,F03) subname//'-----------------------------------------'
     692           0 :      call icepack_warnings_add(warnstr)
     693             :    end if
     694             :  
     695      303298 : END SUBROUTINE shr_orb_params
     696             : 
     697             : !===============================================================================
     698             : 
     699      882945 : SUBROUTINE shr_orb_decl(calday ,eccen ,mvelpp ,lambm0 ,obliqr ,delta ,eccf)
     700             : 
     701             : !-------------------------------------------------------------------------------
     702             : !
     703             : ! Compute earth/orbit parameters using formula suggested by
     704             : ! Duane Thresher.
     705             : !
     706             : !---------------------------Code history----------------------------------------
     707             : !
     708             : ! Original version:  Erik Kluzek
     709             : ! Date:              Oct/1997
     710             : !
     711             : !-------------------------------------------------------------------------------
     712             : 
     713             :    !------------------------------Arguments--------------------------------
     714             :    real   (dbl_kind),intent(in)  :: calday ! Calendar day, including fraction
     715             :    real   (dbl_kind),intent(in)  :: eccen  ! Eccentricity
     716             :    real   (dbl_kind),intent(in)  :: obliqr ! Earths obliquity in radians
     717             :    real   (dbl_kind),intent(in)  :: lambm0 ! Mean long of perihelion at the 
     718             :                                               ! vernal equinox (radians)
     719             :    real   (dbl_kind),intent(in)  :: mvelpp ! moving vernal equinox longitude
     720             :                                               ! of perihelion plus pi (radians)
     721             :    real   (dbl_kind),intent(out) :: delta  ! Solar declination angle in rad
     722             :    real   (dbl_kind),intent(out) :: eccf   ! Earth-sun distance factor (ie. (1/r)**2)
     723             :  
     724             :    !---------------------------Local variables-----------------------------
     725             :    real   (dbl_kind),parameter :: dayspy = 365.0_dbl_kind  ! days per year
     726             :    real   (dbl_kind),parameter :: ve     = 80.5_dbl_kind   ! Calday of vernal equinox
     727             :                                                      ! assumes Jan 1 = calday 1
     728             :  
     729      303198 :    real   (dbl_kind) ::   lambm  ! Lambda m, mean long of perihelion (rad)
     730      303198 :    real   (dbl_kind) ::   lmm    ! Intermediate argument involving lambm
     731      303198 :    real   (dbl_kind) ::   lamb   ! Lambda, the earths long of perihelion
     732      303198 :    real   (dbl_kind) ::   invrho ! Inverse normalized sun/earth distance
     733      303198 :    real   (dbl_kind) ::   sinl   ! Sine of lmm
     734             :  
     735             :    character(len=*),parameter :: subname='(shr_orb_decl)'
     736             : 
     737             :    ! Compute eccentricity factor and solar declination using
     738             :    ! day value where a round day (such as 213.0) refers to 0z at
     739             :    ! Greenwich longitude.
     740             :    !
     741             :    ! Use formulas from Berger, Andre 1978: Long-Term Variations of Daily
     742             :    ! Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci.
     743             :    ! 35:2362-2367.
     744             :    !
     745             :    ! To get the earths true longitude (position in orbit; lambda in Berger 
     746             :    ! 1978) which is necessary to find the eccentricity factor and declination,
     747             :    ! must first calculate the mean longitude (lambda m in Berger 1978) at
     748             :    ! the present day.  This is done by adding to lambm0 (the mean longitude
     749             :    ! at the vernal equinox, set as March 21 at noon, when lambda=0; in radians)
     750             :    ! an increment (delta lambda m in Berger 1978) that is the number of
     751             :    ! days past or before (a negative increment) the vernal equinox divided by
     752             :    ! the days in a model year times the 2*pi radians in a complete orbit.
     753             :  
     754      882945 :    lambm = lambm0 + (calday - ve)*2._dbl_kind*pi/dayspy
     755      882945 :    lmm   = lambm  - mvelpp
     756             :  
     757             :    ! The earths true longitude, in radians, is then found from
     758             :    ! the formula in Berger 1978:
     759             :  
     760      882945 :    sinl  = sin(lmm)
     761             :    lamb  = lambm  + eccen*(2._dbl_kind*sinl + eccen*(1.25_dbl_kind*sin(2._dbl_kind*lmm)  &
     762      882945 :    &     + eccen*((13.0_dbl_kind/12.0_dbl_kind)*sin(3._dbl_kind*lmm) - 0.25_dbl_kind*sinl)))
     763             :  
     764             :    ! Using the obliquity, eccentricity, moving vernal equinox longitude of
     765             :    ! perihelion (plus), and earths true longitude, the declination (delta)
     766             :    ! and the normalized earth/sun distance (rho in Berger 1978; actually inverse
     767             :    ! rho will be used), and thus the eccentricity factor (eccf), can be 
     768             :    ! calculated from formulas given in Berger 1978.
     769             :  
     770      882945 :    invrho = (1._dbl_kind + eccen*cos(lamb - mvelpp)) / (1._dbl_kind - eccen*eccen)
     771             :  
     772             :    ! Set solar declination and eccentricity factor
     773             :  
     774      882945 :    delta  = asin(sin(obliqr)*sin(lamb))
     775      882945 :    eccf   = invrho*invrho
     776             :  
     777      882945 :    return
     778             :  
     779             : END SUBROUTINE shr_orb_decl
     780             : #endif
     781             : 
     782             : !=======================================================================
     783             :  
     784             :       end module icepack_orbital
     785             :  
     786             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd