MAPL_sun_uc.F90 Source File


This file depends on

sourcefile~~mapl_sun_uc.f90~~EfferentGraph sourcefile~mapl_sun_uc.f90 MAPL_sun_uc.F90 sourcefile~base_base.f90 Base_Base.F90 sourcefile~mapl_sun_uc.f90->sourcefile~base_base.f90 sourcefile~constants.f90 Constants.F90 sourcefile~mapl_sun_uc.f90->sourcefile~constants.f90 sourcefile~mapl_comms.f90 MAPL_Comms.F90 sourcefile~mapl_sun_uc.f90->sourcefile~mapl_comms.f90 sourcefile~mapl_exceptionhandling.f90 MAPL_ExceptionHandling.F90 sourcefile~mapl_sun_uc.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~mapl_io.f90 MAPL_IO.F90 sourcefile~mapl_sun_uc.f90->sourcefile~mapl_io.f90 sourcefile~pflogger_stub.f90 pflogger_stub.F90 sourcefile~mapl_sun_uc.f90->sourcefile~pflogger_stub.f90 sourcefile~base_base.f90->sourcefile~constants.f90 sourcefile~mapl_keywordenforcer.f90 MAPL_KeywordEnforcer.F90 sourcefile~base_base.f90->sourcefile~mapl_keywordenforcer.f90 sourcefile~mapl_range.f90 MAPL_Range.F90 sourcefile~base_base.f90->sourcefile~mapl_range.f90 sourcefile~maplgrid.f90 MaplGrid.F90 sourcefile~base_base.f90->sourcefile~maplgrid.f90 sourcefile~internalconstants.f90 InternalConstants.F90 sourcefile~constants.f90->sourcefile~internalconstants.f90 sourcefile~mathconstants.f90 MathConstants.F90 sourcefile~constants.f90->sourcefile~mathconstants.f90 sourcefile~physicalconstants.f90 PhysicalConstants.F90 sourcefile~constants.f90->sourcefile~physicalconstants.f90 sourcefile~mapl_comms.f90->sourcefile~base_base.f90 sourcefile~mapl_comms.f90->sourcefile~constants.f90 sourcefile~mapl_comms.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~shmem.f90 Shmem.F90 sourcefile~mapl_comms.f90->sourcefile~shmem.f90 sourcefile~mapl_errorhandling.f90 MAPL_ErrorHandling.F90 sourcefile~mapl_exceptionhandling.f90->sourcefile~mapl_errorhandling.f90 sourcefile~mapl_throw.f90 MAPL_Throw.F90 sourcefile~mapl_exceptionhandling.f90->sourcefile~mapl_throw.f90 sourcefile~binio.f90 BinIO.F90 sourcefile~mapl_io.f90->sourcefile~binio.f90 sourcefile~fileioshared.f90 FileIOShared.F90 sourcefile~mapl_io.f90->sourcefile~fileioshared.f90 sourcefile~ncio.f90 NCIO.F90 sourcefile~mapl_io.f90->sourcefile~ncio.f90 sourcefile~pfl_keywordenforcer.f90 PFL_KeywordEnforcer.F90 sourcefile~pflogger_stub.f90->sourcefile~pfl_keywordenforcer.f90 sourcefile~wraparray.f90 WrapArray.F90 sourcefile~pflogger_stub.f90->sourcefile~wraparray.f90 sourcefile~binio.f90->sourcefile~base_base.f90 sourcefile~binio.f90->sourcefile~mapl_comms.f90 sourcefile~binio.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~binio.f90->sourcefile~fileioshared.f90 sourcefile~binio.f90->sourcefile~shmem.f90 sourcefile~mapl_sort.f90 MAPL_Sort.F90 sourcefile~binio.f90->sourcefile~mapl_sort.f90 sourcefile~fileioshared.f90->sourcefile~base_base.f90 sourcefile~fileioshared.f90->sourcefile~mapl_comms.f90 sourcefile~fileioshared.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~fileioshared.f90->sourcefile~shmem.f90 sourcefile~fileioshared.f90->sourcefile~mapl_sort.f90 sourcefile~mapl_errorhandling.f90->sourcefile~mapl_throw.f90 sourcefile~mapl_range.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~maplgrid.f90->sourcefile~constants.f90 sourcefile~maplgrid.f90->sourcefile~pflogger_stub.f90 sourcefile~maplgrid.f90->sourcefile~mapl_errorhandling.f90 sourcefile~maplgrid.f90->sourcefile~mapl_keywordenforcer.f90 sourcefile~maplgrid.f90->sourcefile~mapl_sort.f90 sourcefile~ncio.f90->sourcefile~base_base.f90 sourcefile~ncio.f90->sourcefile~mapl_comms.f90 sourcefile~ncio.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~ncio.f90->sourcefile~fileioshared.f90 sourcefile~ncio.f90->sourcefile~shmem.f90 sourcefile~ncio.f90->sourcefile~mapl_sort.f90 sourcefile~pfio.f90 pFIO.F90 sourcefile~ncio.f90->sourcefile~pfio.f90 sourcefile~physicalconstants.f90->sourcefile~mathconstants.f90 sourcefile~shmem.f90->sourcefile~constants.f90

Files dependent on this one

sourcefile~~mapl_sun_uc.f90~~AfferentGraph sourcefile~mapl_sun_uc.f90 MAPL_sun_uc.F90 sourcefile~base.f90 Base.F90 sourcefile~base.f90->sourcefile~mapl_sun_uc.f90 sourcefile~genericcplcomp.f90 GenericCplComp.F90 sourcefile~genericcplcomp.f90->sourcefile~mapl_sun_uc.f90 sourcefile~mapl_generic.f90 MAPL_Generic.F90 sourcefile~mapl_generic.f90->sourcefile~mapl_sun_uc.f90 sourcefile~mapl_generic.f90->sourcefile~genericcplcomp.f90 sourcefile~comp_testing_driver.f90 Comp_Testing_Driver.F90 sourcefile~comp_testing_driver.f90->sourcefile~mapl_generic.f90 sourcefile~extdatagridcompmod.f90 ExtDataGridCompMod.F90 sourcefile~extdatagridcompmod.f90->sourcefile~mapl_generic.f90 sourcefile~extdatagridcompng.f90 ExtDataGridCompNG.F90 sourcefile~extdatagridcompng.f90->sourcefile~mapl_generic.f90 sourcefile~mapl.f90 MAPL.F90 sourcefile~mapl.f90->sourcefile~base.f90 sourcefile~mapl.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_capgridcomp.f90 MAPL_CapGridComp.F90 sourcefile~mapl_capgridcomp.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_geosatmaskmod.f90 MAPL_GeosatMaskMod.F90 sourcefile~mapl_geosatmaskmod.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_historycollection.f90 MAPL_HistoryCollection.F90 sourcefile~mapl_historycollection.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_historygridcomp.f90 MAPL_HistoryGridComp.F90 sourcefile~mapl_historygridcomp.f90->sourcefile~genericcplcomp.f90 sourcefile~mapl_historygridcomp.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_nuopcwrappermod.f90 MAPL_NUOPCWrapperMod.F90 sourcefile~mapl_nuopcwrappermod.f90->sourcefile~base.f90 sourcefile~mapl_orbgridcompmod.f90 MAPL_OrbGridCompMod.F90 sourcefile~mapl_orbgridcompmod.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_stationsamplermod.f90 MAPL_StationSamplerMod.F90 sourcefile~mapl_stationsamplermod.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_trajectorymod.f90 MAPL_TrajectoryMod.F90 sourcefile~mapl_trajectorymod.f90->sourcefile~mapl_generic.f90 sourcefile~test_cfio_bundle.pf Test_CFIO_Bundle.pf sourcefile~test_cfio_bundle.pf->sourcefile~base.f90 sourcefile~tstqsat.f90 tstqsat.F90 sourcefile~tstqsat.f90->sourcefile~base.f90 sourcefile~ut_extdata.f90 ut_ExtData.F90 sourcefile~ut_extdata.f90->sourcefile~base.f90 sourcefile~utcfio_bundle.f90 utCFIO_Bundle.F90 sourcefile~utcfio_bundle.f90->sourcefile~base.f90

Source Code

!------------------------------------------------------------------------------
!               Global Modeling and Assimilation Office (GMAO)                !
!                    Goddard Earth Observing System (GEOS)                    !
!                                 MAPL Component                              !
!------------------------------------------------------------------------------
!
#include "MAPL_ErrLog.h"
!
!>
!### MODULE: `MAPL_SunMod`
!
! Author: GMAO SI-Team
!
! This class is intended to manage the sun`s position and provide
! the insolation at the top of the atmosphere. The main method
! is `GEOS_SunGetInsolation`, which depends on an Orbit object.
! The Orbit object defines this class and has public opaque type `GEOS_SunOrbit`.
! Methods are provided for creating it, destroying it, and making various queries.
!
module MAPL_SunMod

! !USES:

  use ESMF
  use MAPL_Constants
  use MAPL_BaseMod
  use MAPL_IOMod
  use MAPL_CommsMod
  use MAPL_ExceptionHandling
  use netcdf
  use, intrinsic :: iso_fortran_env, only: REAL64
  use pflogger, only: logging, Logger

  implicit none
  private

! !PUBLIC MEMBER FUNCTIONS:

  public MAPL_SunOrbitCreate
  public MAPL_SunOrbitCreateFromConfig
  public MAPL_SunOrbitCreated
  public MAPL_SunOrbitDestroy
  public MAPL_SunOrbitQuery
  public MAPL_SunGetInsolation
  public MAPL_SunGetSolarConstant
  public MAPL_SunGetDaylightDuration
  public MAPL_SunGetDaylightDurationMax
  public MAPL_SunGetLocalSolarHourAngle

! !PUBLIC TYPES:

  public MAPL_SunOrbit

!EOP

  integer, public, parameter :: MAPL_SunAutumnalEquinox = 1
  integer, public, parameter :: MAPL_SunWinterSolstice  = 2
  integer, public, parameter :: MAPL_SunVernalEquinox   = 3
  integer, public, parameter :: MAPL_SunSummerSolstice  = 4
  integer, public, parameter :: MAPL_SunDailyMean       = 5
  integer, public, parameter :: MAPL_SunAnnualMean      = 6

! Default solar orbital system parameters (private).
! Dont change these unless you know what you are doing.
! They are appropriate for the current modern epoch circa 2000.
! -------------------------------------------------------------

   ! Parameters of old orbital system (tabularized intercalation cycle)
   ! ------------------------------------------------------------------
   real, parameter    :: DEFAULT_ORBIT_ECCENTRICITY     = 0.0167
   real, parameter    :: DEFAULT_ORBIT_OBLIQUITY        = 23.45   ! degrees
   real, parameter    :: DEFAULT_ORBIT_PERIHELION       = 102.0   ! degrees
   integer, parameter :: DEFAULT_ORBIT_EQUINOX          = 80      ! days

   ! Parameters of new orbital system (analytic two-body), which allows some
   ! time-varying behavior, namely, linear variation in LAMBDAP, ECC, and OBQ.
   ! -------------------------------------------------------------------------

   ! Fixed anomalistic year length in mean solar days
   real, parameter    :: DEFAULT_ORB2B_YEARLEN          = 365.2596

   ! Reference date and time for orbital parameters
   ! (defaults to J2000 = 01Jan2000 12:00:00 TT = 11:58:56 UTC)
   integer, parameter :: DEFAULT_ORB2B_REF_YYYYMMDD     = 20000101
   integer, parameter :: DEFAULT_ORB2B_REF_HHMMSS       = 115856

   ! Orbital eccentricity at reference date
   real, parameter    :: DEFAULT_ORB2B_ECC_REF          = 0.016710
   ! Rate of change of orbital eccentricity per Julian century
   real, parameter    :: DEFAULT_ORB2B_ECC_RATE         = -4.2e-5

   ! Earth's obliquity (axial tilt) at reference date [degrees]
   real, parameter    :: DEFAULT_ORB2B_OBQ_REF          = 23.44
   ! Rate of change of obliquity [degrees per Julian century]
   real, parameter    :: DEFAULT_ORB2B_OBQ_RATE         = -1.3e-2

   ! Longitude of perihelion at reference date [degrees]
   !   (from March equinox to perihelion in direction of earth's motion)
   real, parameter    :: DEFAULT_ORB2B_LAMBDAP_REF      = 282.947
   ! Rate of change of LAMBDAP [degrees per Julian century]
   !   (Combines both equatorial and ecliptic precession)
   real, parameter    :: DEFAULT_ORB2B_LAMBDAP_RATE     = 1.7195

   ! March Equinox date and time
   ! (defaults to March 20, 2000 at 07:35:00 UTC)
   integer, parameter :: DEFAULT_ORB2B_EQUINOX_YYYYMMDD = 20000320
   integer, parameter :: DEFAULT_ORB2B_EQUINOX_HHMMSS   =  73500

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

  interface MAPL_SunGetInsolation
     module procedure SOLAR_1D
     module procedure SOLAR_2D
     module procedure SOLAR_ARR_INT
  end interface

  interface MAPL_SunGetSolarConstant
     module procedure MAPL_SunGetSolarConstantByYearDoY
     module procedure MAPL_SunGetSolarConstantByTime
     module procedure MAPL_SunGetSolarConstantFromNetcdfFile
     module procedure MAPL_SunGetSolarConstantFromNRLFile
  end interface

  type MAPL_SunOrbit
     private
     logical                          :: CREATED = .FALSE.
     type(ESMF_Clock)                 :: CLOCK
     real                             :: OB, ECC, PER, YEARLEN
     integer                          :: EQNX, YEARS_PER_CYCLE, DAYS_PER_CYCLE
     real, pointer, dimension(:)      :: ZC => null()
     real, pointer, dimension(:)      :: ZS => null()
     real, pointer, dimension(:)      :: PP => null()
     real, pointer, dimension(:)      :: TH => null()
     real, pointer, dimension(:)      :: ET => null()
     logical                          :: EOT
     logical                          :: FIX_SUN
     logical                          :: ANAL2B
     real                             :: ORB2B_YEARLEN
     type(ESMF_Time)                  :: ORB2B_TIME_REF
     real                             :: ORB2B_ECC_REF
     real                             :: ORB2B_ECC_RATE
     real                             :: ORB2B_OBQ_REF
     real                             :: ORB2B_OBQ_RATE
     real                             :: ORB2B_LAMBDAP_REF
     real                             :: ORB2B_LAMBDAP_RATE
     type(ESMF_Time)                  :: ORB2B_TIME_EQUINOX
     real                             :: ORB2B_OMG0
     type(ESMF_Time)                  :: ORB2B_TIME_PERI
  end type MAPL_SunOrbit

contains

!==========================================================================
      ! utlity functions

      ! rectify to [-pi,+pi)
      function RECT_PMPI(X)
        real :: X, RECT_PMPI
        RECT_PMPI = MODULO( X + MAPL_PI, 2*MAPL_PI ) - MAPL_PI
      end function

      ! true anomaly from eccentric anomaly
      function calcTAfromEA(EA,EAFAC) result(TA)
        real :: EA, EAFAC, TA
        TA = 2. * atan( tan(EA / 2.) / EAFAC)
      end function

      ! eccentric anomaly from true anomaly
      function calcEAfromTA(TA,EAFAC) result(EA)
        real :: TA, EAFAC, EA
        EA = 2. * atan( EAFAC * tan(TA / 2.))
      end function

      ! mean anomaly from eccentric anomaly (Kepler's equation)
      function calcMAfromEA(EA,ECC) result(MA)
        real :: EA, ECC, MA
        MA = EA - ECC * sin(EA)
      end function

      ! eccentric anomaly from mean anomaly
      ! (invert Kepler's equation by Newton-Raphson)
      subroutine invert_Keplers_Newton( &
          MA, ECC, &
          EA, dE, nits, &
          tol, max_its)

        real,              intent(in ) :: MA
        real,              intent(in ) :: ECC

        real,              intent(out) :: EA
        real,              intent(out) :: dE
        integer,           intent(out) :: nits

        real,    optional, intent(in ) :: tol
        integer, optional, intent(in ) :: max_its

        real :: f, df, tol_
        integer :: max_its_

        if (present(tol)) then
          tol_ = tol
        else
          tol_ = 1.e-10
        endif

        if (present(max_its)) then
          max_its_ = max_its
        else
          max_its_ = 10
        endif

        EA = MA
        do nits = 1, max_its_
          f = EA - ECC * sin(EA) - MA
          df = 1. - ECC * cos(EA)
          dE = -f / df
          EA = EA + dE
          if (abs(dE) < tol_) exit
        enddo

      end subroutine

      ! Earth-Sun distance from true anomaly
      function calcRadfromTA(TA,ECC,OMSQECC) result(Rad)
        real :: TA, ECC, OMSQECC, Rad
        Rad = OMSQECC / (1. + ECC * cos(TA))
      end function

!==========================================================================
!>
! Integrates the earth`s orbit and stores the necessary
! parameters to easily compute the earth`s position for each day
! of the full (usually 4-year) intercalation cycle.
! The orbital parameters are passed as arguments.
! The full calendar intercalation cycle is obtained from the
! ESMF clock passed as an argument. This becomes the orbit`s
! attached clock. Currently we assume a single intercalation.
!
! A good introduction to celestial mechanics for understanding
! this code can be found in Blanco & McCuskey, 1961: "Basic
! Physics of the Solar System", hereafter BM.
!
!```
! CLOCK: The orbit will depend on the calendar in this clock
!        This is used for the length of year, to set intercalation cycle.
!
! ECCENTRICITY: Eccentricity of the Earth`s orbit.
!
! PERIHELION: Longitude of perihelion, measured in degrees from
!             autumnal equinox in the direction of the Earth`s motion.
!
! OBLIQUITY: Tilt of the Earth`s rotation axis from a
!            normal to the plane of the orbit. In degrees.
!
! EQUINOX: Day of year of vernal equinox.
!          Equinox is assumed to occur at 0Z on this day on the
!          first year of the cycle.
!
! EOT: Apply Equation of Time correction?
!
! ORBIT_ANAL2B: New orbital system (analytic two-body) allows some
!               time-varying behavior, namely, linear time variation in LAMBDAP,
!               ECC, and OBQ. If .TRUE., the following ORB2B parameters are used
!               and only CLOCK and EOT above are used, i.e., the ECCENTRICITY,
!               OBLIQUITY, PERIHELION and EQUINOX above are NOT used and are
!               replaced by the relevant ORB2B parameters below.
!
! ORB2B_YEARLEN: Fixed anomalistic year length in mean solar days.
!
! ORB2B_REF_YYYYMMDD: Reference date for orbital parameters.
!
! ORB2B_REF_HHMMSS: Reference time for orbital parameters.
!
! ORB2B_ECC_REF: Orbital eccentricity at reference date.
!
! ORB2B_ECC_RATE: Rate of change of orbital eccentricity per Julian century.
!
! ORB2B_OBQ_REF: Earth's obliquity (axial tilt) at reference date [degrees].
!
! ORB2B_OBQ_RATE: Rate of change of obliquity [degrees per Julian century].
!
! ORB2B_LAMBDAP_REF: Longitude of perihelion at reference date [degrees]
!                    from March equinox to perihelion in direction of earth's motion).
!
! ORB2B_LAMBDAP_RATE: Rate of change of LAMBDAP [degrees per Julian century]
!                     (Combines both equatorial and ecliptic precession).
!
! ORB2B_EQUINOX_YYYYMMDD: March equinox date.
!
! ORB2B_EQUINOX_HHMMSS: March equinox time.
!```
!
type(MAPL_SunOrbit) function MAPL_SunOrbitCreate(CLOCK,                  &
                                                 ECCENTRICITY,           &
                                                 OBLIQUITY,              &
                                                 PERIHELION,             &
                                                 EQUINOX,                &
                                                 EOT,                    &
                                                 ORBIT_ANAL2B,           &
                                                 ORB2B_YEARLEN,          &
                                                 ORB2B_REF_YYYYMMDD,     &
                                                 ORB2B_REF_HHMMSS,       &
                                                 ORB2B_ECC_REF,          &
                                                 ORB2B_ECC_RATE,         &
                                                 ORB2B_OBQ_REF,          &
                                                 ORB2B_OBQ_RATE,         &
                                                 ORB2B_LAMBDAP_REF,      &
                                                 ORB2B_LAMBDAP_RATE,     &
                                                 ORB2B_EQUINOX_YYYYMMDD, &
                                                 ORB2B_EQUINOX_HHMMSS,   &
                                                 FIX_SUN,                &
                                                                      RC )

! !ARGUMENTS:

 type(ESMF_Clock)  , intent(IN ) :: CLOCK
 real              , intent(IN ) :: ECCENTRICITY
 real              , intent(IN ) :: OBLIQUITY
 real              , intent(IN ) :: PERIHELION
 integer           , intent(IN ) :: EQUINOX
 logical           , intent(IN ) :: EOT
 logical           , intent(IN ) :: ORBIT_ANAL2B
 real              , intent(IN ) :: ORB2B_YEARLEN
 integer           , intent(IN ) :: ORB2B_REF_YYYYMMDD
 integer           , intent(IN ) :: ORB2B_REF_HHMMSS
 real              , intent(IN ) :: ORB2B_ECC_REF
 real              , intent(IN ) :: ORB2B_ECC_RATE
 real              , intent(IN ) :: ORB2B_OBQ_REF
 real              , intent(IN ) :: ORB2B_OBQ_RATE
 real              , intent(IN ) :: ORB2B_LAMBDAP_REF
 real              , intent(IN ) :: ORB2B_LAMBDAP_RATE
 integer           , intent(IN ) :: ORB2B_EQUINOX_YYYYMMDD
 integer           , intent(IN ) :: ORB2B_EQUINOX_HHMMSS
 logical, optional , intent(IN ) :: FIX_SUN
 integer, optional , intent(OUT) :: RC

!EOPI

! =======================================
! PMN Dec 2019: Notes on Equation of Time
! =======================================
! (Part of a more complete analysis available in orbit.pdf)
!
! @ Introduction:
!
! The Earth rotates on its axis with a period T_S called the sidereal
! day (after the Latin for "star", since it is the rotation period of
! the Earth with respect to distant stars). T_S is slightly shorter
! than the so-called mean solar day, or clock day, of duration T_M =
! 86400 seconds. This is because the Earth is a prograde planet, i.e.,
! it rotates about its axis in the same sense (counterclockwise look-
! ing down on the North Pole) as it orbits the sun. Specifically, say
! the sun crosses the meridian of some location at a particular time.
! And imagine there is a distant star directly behind the sun at that
! moment. After one sidereal day the location will rotate 360 degrees
! about the earth's axis and the distant star will again cross its
! merdian. But during that period the earth will have moved a small
! counterclockwise distance around its orbit and so it will take a
! small additional rotation of the earth for the sun to also cross
! the meridian and thereby complete a solar day.
!
! Put another way, a solar day is slightly longer than a sidereal day
! because the sun appears to move slowly eastward across the celestial
! sphere with respect to distant stars as the year passes. The path
! of this motion is called the ecliptic. And clearly what governs the
! length of a solar day is the apparent velocity of the sun along the
! ecliptic, or, more particularly, the equatorial component of that
! velocity. But both the magnitude and equatorial component of the solar
! ecliptic velocity change during the year, the former because the earth's
! orbit is elliptical, not circular, and the latter because the earth's
! axis of rotation is tilted with respect to the orbital (ecliptic) plane.
! Thus the length of a solar day changes during the year. While these
! factors cause only a small perturbation to the length of the solar
! day (less than 30 seconds), the perturbations accumulate so that at
! different times of the year apparent solar time ("sundial time") and
! mean solar time ("clock time") can differ by as much as ~15 minutes.
! This difference is called the Equation of Time.
!
! To be more rigorous, consider a fictitious "Mean Sun" that moves at
! constant eastward speed around the celestial equator, completing a
! full revolution in a year, namely in the period Y * T_M, where Y is
! the number of mean solar days in a year (e.g., 365.25). Thus, in one
! mean solar day, T_M, the mean sun has moved an angle 2*PI/Y eastward.
! Hence, beyond one full earth revolution, period T_S, an additional
! earth rotation of (T_M-T_S) * 2*PI/T_S = 2*PI/Y is required to "catch
! up with the moving sun", as described earlier. Hence T_M - T_S = T_S / Y
! and so
!
! T_M = T_S (Y+1)/Y,..... (1)
!
! a constant (near unity) multiple of the fixed sidereal day. T_M is
! the length of the solar day for the "mean sun", or the "mean solar
! day". Because it is invariant during the year, it is convenient for
! timekeeping, and forms the basis for "mean solar time", which at
! Greenwich is essentially UTC. By *definition*, T_M = 24h = 86400s.
! That is, what we know as "hours", "minutes" and "seconds", are just
! convenient integer fractions of the mean solar day. In these units,
! the sidereal day T_S is approximately 23h 56m 4s.
!
! The solar zenith angle calculation (in MAPL_SunGetInsolation) needs
! the true local solar hour angle, h_T, which is the angle, measured
! westward along the equator, from the local meridian to the true sun.
! This is just the Greenwich solar hour angle, H_T, plus the longitude,
! so we will henceforth work exclusively with Greenwich hour angles.
! We should use the hour angle of the *true* sun, H_T, but a common
! approximation replaces this with the hour angle of the mean sun
!
! H_M = 2*PI*(u-0.5),..... (2)
!
! where u is UTC time (in days) and the offset is needed because the mean
! solar hour angle is zero at "noon". If more accuracy is required, the
! hour angle of the true sun is typically obtained as a small correction
! to H_M called the Equation of time, EOT:
!
! H_T = H_M + EOT,  where  EOT = H_T - H_M.
!
! As discussed above, EOT corrects for two factors:
! (a) the variable speed of the earth in its elliptical orbit about
!     the sun (e.g., moving fastest at perihelion), and
! (b) the tilt of the earth's axis of rotation wrt its orbital plane
!     (the "obliquity"), which causes the equatorial projection of
!     the sun's apparent ecliptic motion to vary with the season
!     (e.g., being parallel to the equator at the solstices.)
!
! @ Derivation of Equation of Time:
!
! We can write
!
! H_T  = H_1 - (H_1 - H_T) = H_1 - a_T,
!
! where H_1 is the Greenwich hour angle of the First Point of Aries (the
! location of the vernal equinox, denoted "1PoA"), and is also known as the
! Greenwich Sidereal hour angle, and where a_T is the right ascension of the
! true sun (since the right ascension of any object is just the difference
! between the hour angles of 1PoA and the object). Hence,
!
! EOT = H_1 - H_M - a_T...... (3)
!
! All three terms on the right of (3) are time variable: a_T changes slowly
! throughout the year, and is known from the earth-sun two-body elliptical
! orbit solution, while H_1 and H_M vary rapidly with earth's rotation. H_M
! has a period of one mean solar day, T_M, and H_1 has a period of one
! sidereal day, T_S.
!
! It may seem from from (2) that the mean sun and its hour angle are fully
! specified. That, in fact, is not yet the case: (2) is really just a def-
! inition of UTC, namely, that one UTC day is one mean solar day and that
! the time of culmination of the mean sun, what we call "noon", occurs at
! UTC 12h00m. What we are still at liberty to do is specify the phasing of
! the mean sun in its equatorial orbit, e.g., by specifying the time u_R
! at which the mean sun passes through 1PoA (both on the equator). At this
! time, H_1(u_R) = H_M(u_R), and so
!
! H_1(u) - H_M(u) = 2*PI*(u-u_R)*(Y+1)/Y - 2*PI*(u-u_R)
!                 = 2*PI * (u-u_R) / Y
!                       = MA(u) - MA(u_R),... (4)
!
! where MA(u) = 2*PI * (u-u_P) / Y is the so-called "mean anomaly", known
! from the earth-sun two-body orbital solution, and u_P is the time of
! perihelion. Thus, to fully determine EOT, through (3) and (4), we need
! only to specify MA(u_R).
!
! To understand the mean anomaly MA, consider the standard two-body earth-
! sun problem in which the earth E moves in an elliptical orbit about the
! sun S at one focus, all in the ecliptic plane. The point on the ellipse
! closest to S is called the perihelion P. Obviously, the center of the
! ellipse O, the focus S and the perihelion P are co-linear, the so-called
! major axis of the ellipse. Additionally, let C be the circumscribing circle
! of the ellipse, with center O and passing through P (and the corresponding
! aphelion A). By Kepler's Second Law, the sun-earth vector sweeps out equal
! areas in equal times, so the fractional area of the elliptical sector PSE
! is a linear function of time, being zero at perihelion and one a year later.
! Specifically, this fractional area is none other than the scaled mean anomaly
! MA(u)/(2*PI) = (u - u_P) / Y. Clearly MA(u) can also be interpreted as an
! angle, the angle POQ of a point Q orbiting on the circumcircle C at constant
! speed in the same direction as the earth, also with a yearly period, and
! passing through P at the same time u_P as the earth. Thus the point Q can
! be conceptualized as a sort of "mean earth" orbiting a "second mean sun"
! (different from M above) at O. Note that while the angle MA(u) = angle POQ
! of this mean earth at time u is a linear function of time, the corresponding
! angle of the real earth, namely TA(u) = angle PSE, called the true anomaly,
! is a non-linear function of time, since the real earth has a variable speed
! in its elliptical orbit, e.g., moving faster at perihelion, so that its areal
! fraction is linear in time. The relationship between MA(u) and TA(u) is known
! from the orbital solution. Finally, the ecliptic longitude of the earth,
! lambda = angle 1SE is the angle at the sun, measured in the same direction
! as the earth's motion, from the First Point of Aries to the earth. Then
!
!    TA(u) = angle PSE(u) = angle PS1 + angle 1SE(u) = lambda(u) - lambda_P,
!
! where lambda_P = lambda(u_P) = angle 1SP = -angle PS1 is known as the
! longitude of perihelion, and is currently about 283 deg, or equivalently
! -77 deg.
!
! With this background, we can understand the quantity MA(u_R) we are trying
! to specify. If we *choose*
!
!                      MA(u_R) = -lambda_P = angle PS1... (5)
!                      <==> angle POQ(u_R) = angle PS1,
!
! then at u_R, viewed from the mean earth Q, the second (ecliptic) mean sun
! O is in direction of 1PoA. And at that same time, by definition of u_R,
! the first (equatorial) mean sun M, as seen from the real earth E, is also
! in direction of 1PoA.
!
! I (PMN) have verified that the choice (5) for MA(u_R) leads to zero mean
! EOT to first order in the eccentricity, e. I cannot say, at this point,
! that it is generally true for any order in e. I add below a final explicit
! enforcement of zero mean EOT to the code. [I found that the order e^2 term
! for the mean EOT was NOT zero, and was larger than the mean EOT produced
! by this code (which is valid for all orders in e) before any explicity
! correction of the mean to zero. This suggests: (a) that I made a mistake
! in my calculations, or (b) that higer order e terms provide some cancel-
! ation.]
!
! Hence, finally,
!
! EOT(u) = MA(u) + PRHV - a_T(u)... (6)
!
! where PRHV is the name for lamba_P in the code.
! ===========================================================================

! Locals

      character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitCreate"

      real(kind=REAL64)  :: YEARLEN
      integer :: K, KP, YEARS_PER_CYCLE, DAYS_PER_CYCLE
      real(kind=REAL64)  :: TREL, T1, T2, T3, T4
      real(kind=REAL64)  :: SOB, COB, OMG0, OMG, PRH, PRHV
      real    :: OMECC, OPECC, OMSQECC, EAFAC
      real(kind=REAL64)  :: TA, EA, MA, TRRA, MNRA
      real    :: meanEOT
      type(MAPL_SunOrbit) :: ORBIT
      integer :: STATUS

      integer :: year, month, day, hour, minute, second
      real(ESMF_KIND_R8) :: days
      real :: ECC_EQNX, LAMBDAP_EQNX, EAFAC_EQNX
      real :: TA_EQNX, EA_EQNX, MA_EQNX
      type(ESMF_TimeInterval) :: DT

      ! record inputs needed by both orbit methods
      ORBIT%CLOCK  = CLOCK
      ORBIT%EOT    = EOT
      ORBIT%ANAL2B = ORBIT_ANAL2B

      ! Analytic two-body orbit?
      if (ORBIT_ANAL2B) then

        ! record inputs in ORBIT type
        associate(D2R => MAPL_DEGREES_TO_RADIANS)
          ORBIT%ORB2B_YEARLEN      = ORB2B_YEARLEN
          ORBIT%ORB2B_ECC_REF      = ORB2B_ECC_REF
          ORBIT%ORB2B_OBQ_REF      = ORB2B_OBQ_REF      * D2R           ! radians
          ORBIT%ORB2B_LAMBDAP_REF  = ORB2B_LAMBDAP_REF  * D2R           ! radians
          ORBIT%ORB2B_ECC_RATE     = ORB2B_ECC_RATE           / 36525.  ! per day
          ORBIT%ORB2B_OBQ_RATE     = ORB2B_OBQ_RATE     * D2R / 36525.  ! radians per day
          ORBIT%ORB2B_LAMBDAP_RATE = ORB2B_LAMBDAP_RATE * D2R / 36525.  ! radians per day
        end associate
        ! record MAPL Time object for REFerence time
        year   = ORB2B_REF_YYYYMMDD / 10000
        month  = mod(ORB2B_REF_YYYYMMDD, 10000) / 100
        day    = mod(ORB2B_REF_YYYYMMDD, 100)
        hour   = ORB2B_REF_HHMMSS / 10000
        minute = mod(ORB2B_REF_HHMMSS, 10000) / 100
        second = mod(ORB2B_REF_HHMMSS, 100)
        call ESMF_TimeSet(ORBIT%ORB2B_TIME_REF, &
          yy=year, mm=month, dd=day, h=hour, m=minute, s=second, rc=STATUS)
        _VERIFY(STATUS)
        ! record MAPL Time object for EQUINOX
        year   = ORB2B_EQUINOX_YYYYMMDD / 10000
        month  = mod(ORB2B_EQUINOX_YYYYMMDD, 10000) / 100
        day    = mod(ORB2B_EQUINOX_YYYYMMDD, 100)
        hour   = ORB2B_EQUINOX_HHMMSS / 10000
        minute = mod(ORB2B_EQUINOX_HHMMSS, 10000) / 100
        second = mod(ORB2B_EQUINOX_HHMMSS, 100)
        call ESMF_TimeSet(ORBIT%ORB2B_TIME_EQUINOX, &
          yy=year, mm=month, dd=day, h=hour, m=minute, s=second, rc=STATUS)
        _VERIFY(STATUS)

        ! time-invariant precalculations
        ORBIT%ORB2B_OMG0 = 2. * MAPL_PI / ORB2B_YEARLEN

        ! at provided ORB2B_TIME_EQUINOX LAMBDA=0 by definition
        call ESMF_TimeIntervalGet( &
          ORBIT%ORB2B_TIME_EQUINOX - ORBIT%ORB2B_TIME_REF, &
          d_r8=days, rc=STATUS)
        _VERIFY(STATUS)
        ECC_EQNX = ORBIT%ORB2B_ECC_REF + days * ORBIT%ORB2B_ECC_RATE
        LAMBDAP_EQNX = ORBIT%ORB2B_LAMBDAP_REF + days * ORBIT%ORB2B_LAMBDAP_RATE
        EAFAC_EQNX = sqrt((1.-ECC_EQNX)/(1.+ECC_EQNX))
        TA_EQNX = -LAMBDAP_EQNX  ! since LAMBDA=0
        EA_EQNX = calcEAfromTA(TA_EQNX,EAFAC_EQNX)
        MA_EQNX = calcMAfromEA(EA_EQNX,ECC_EQNX)

        ! Time of one perihelion (all subsequent ORB2B_YEARLEN apart)
        call ESMF_TimeIntervalSet(DT, d_r8 = dble(MA_EQNX / ORBIT%ORB2B_OMG0))
        ORBIT%ORB2B_TIME_PERI = ORBIT%ORB2B_TIME_EQUINOX - DT

      else

        ! ==================================
        ! Standard tabularized intercalation
        ! ==================================

        ! MJS:  This needs to come from the calendar when the time manager works right.
        YEARLEN = 365.25

        ! Factors involving the orbital parameters
        !-----------------------------------------
        OMECC = 1. - ECCENTRICITY
        OPECC = 1. + ECCENTRICITY
        OMSQECC = 1. - ECCENTRICITY**2 ! pmn: consider changing to line below when zero-diff not issue
!       OMSQECC = OMECC * OPECC
        EAFAC = sqrt(OMECC/OPECC)

        associate(D2R => MAPL_DEGREES_TO_RADIANS)
          OMG0 = 2.*MAPL_PI/YEARLEN
          OMG  = OMG0/sqrt(OMSQECC)**3
          PRH  = PERIHELION*D2R
          SOB  = sin(OBLIQUITY*D2R)
          COB  = cos(OBLIQUITY*D2R)
        end associate

        ! PRH is the ecliptic longitude of the perihelion, measured (at the Sun)
        ! from the autumnal equinox in the direction of the Earth`s orbital motion
        ! (counterclockwise as viewed from ecliptic north pole).
        ! For EOT calculations we will reference the perihelion wrt to the vernal
        ! equinox, PRHV. Of course, the difference is just PI.
        ! pmn: once the EOT code is established and zero-diff not an issue,
        !   consider removing original PRH and changing the original (non-EOT),
        !   code, which employs
        !       cos(Y \pm PI) = -COS(Y)
        ! to use PRHV, namely
        !       -cos(Y-PRH) = cos(Y-PRH-PI) = cos(Y-PRHV)
        PRHV = PRH + MAPL_PI_R8

        ! Compute length of leap cycle
        ! ----------------------------
        if(YEARLEN-int(YEARLEN) > 0.) then
         YEARS_PER_CYCLE = nint(1./(YEARLEN-int(YEARLEN)))
        else
         YEARS_PER_CYCLE = 1
        endif

        DAYS_PER_CYCLE=nint(YEARLEN*YEARS_PER_CYCLE)

        ! save inputs and intercalculation details
        ! ----------------------------------------
        ORBIT%OB              = OBLIQUITY
        ORBIT%ECC             = ECCENTRICITY
        ORBIT%PER             = PERIHELION
        ORBIT%EQNX            = EQUINOX
        ORBIT%YEARLEN         = YEARLEN
        ORBIT%YEARS_PER_CYCLE = YEARS_PER_CYCLE
        ORBIT%DAYS_PER_CYCLE  = DAYS_PER_CYCLE

        ! Allocate orbital cycle outputs
        ! ------------------------------
        ! TH: Ecliptic longitude of the true Sun, TREL [radians]
        ! ZS: Sine of declination
        ! ZC: Cosine of declination
        ! PP: Inverse of square of earth-sun distance [1/(AU**2)]
        ! ET: Equation of time [radians] =
        !     True solar hour angle - Mean solar hour angle

        if(associated(ORBIT%TH)) deallocate(ORBIT%TH)
        allocate(ORBIT%TH(DAYS_PER_CYCLE), stat=status)
        _VERIFY(STATUS)

        if(associated(ORBIT%ZC)) deallocate(ORBIT%ZC)
        allocate(ORBIT%ZC(DAYS_PER_CYCLE), stat=status)
        _VERIFY(STATUS)

        if(associated(ORBIT%ZS)) deallocate(ORBIT%ZS)
        allocate(ORBIT%ZS(DAYS_PER_CYCLE), stat=status)
        _VERIFY(STATUS)

        if(associated(ORBIT%PP)) deallocate(ORBIT%PP)
        allocate(ORBIT%PP(DAYS_PER_CYCLE), stat=status)
        _VERIFY(STATUS)

        if (ORBIT%EOT) then
          if(associated(ORBIT%ET)) deallocate(ORBIT%ET)
          allocate(ORBIT%ET(DAYS_PER_CYCLE), stat=status)
          _VERIFY(STATUS)
        end if

        ! Begin integration at the vernal equinox (K=1, KP=EQUINOX), at
        ! which, by defn, the ecliptic longitude of the true sun is zero.
        ! Right ascension of true sun at EQUINOX is also zero by defn.
        ! --------------------------------------------------------------
        KP           = EQUINOX
        TREL         = 0.
        ORBIT%ZS(KP) = sin(TREL)*SOB
        ORBIT%ZC(KP) = sqrt(1.-ORBIT%ZS(KP)**2)
        ORBIT%PP(KP) = ( (1.-ECCENTRICITY*cos(TREL-PRH)) / OMSQECC )**2
        ORBIT%TH(KP) = TREL

        if (ORBIT%EOT) then
          ! Right ascension of "mean sun" = MA(u) + PRHV [eqn (6) above].
          ! Calcn of True (TA), Eccentric (EA), and Mean Anomaly (MA).
          TA = TREL - PRHV                    ! by defn of TA and PRHV
          EA = 2.d0*atan(EAFAC*tan(TA/2.d0))  ! BM 4-55
          MA = EA - ECCENTRICITY*sin(EA)      ! Kepler's eqn (BM 4-49 ff.)
          MNRA = MA + PRHV                    ! see EOT notes above
          TRRA = 0.                           ! because TREL=0 at Equinox

          ! Equation of Time, ET [radians]
          ! True Solar hour angle = Mean Solar hour angle + ET
          ! (hour angle and right ascension are in reverse direction)
          ORBIT%ET(KP) = rect_pmpi(real(MNRA - TRRA))
        end if

        ! Integrate orbit for entire leap cycle using Runge-Kutta
        ! Mean sun moves at constant speed around Celestial Equator
        ! ---------------------------------------------------------
        do K=2,DAYS_PER_CYCLE
          T1 = dTRELdDAY(TREL       ,OMG,ECCENTRICITY,PRH)
          T2 = dTRELdDAY(TREL+T1*0.5,OMG,ECCENTRICITY,PRH)
          T3 = dTRELdDAY(TREL+T2*0.5,OMG,ECCENTRICITY,PRH)
          T4 = dTRELdDAY(TREL+T3    ,OMG,ECCENTRICITY,PRH)
          KP = mod(KP,DAYS_PER_CYCLE) + 1
          TREL = TREL + (T1 + 2.0*(T2 + T3) + T4) / 6.0
          ORBIT%ZS(KP) = sin(TREL)*SOB
          ORBIT%ZC(KP) = sqrt(1.-ORBIT%ZS(KP)**2)
          ORBIT%PP(KP) = ( (1.-ECCENTRICITY*cos(TREL-PRH)) / OMSQECC )**2
          ORBIT%TH(KP) = TREL
          if (ORBIT%EOT) then
            ! From BM 1-15 and 1-16 with beta=0 (Sun is on ecliptic),
            ! and dividing through by common cos(dec) since it does not
            ! affect the ratio of sin(RA) to cos(RA).
            TRRA = atan2(sin(TREL)*COB,cos(TREL))
            MNRA = MNRA + OMG0
            ORBIT%ET(KP) = rect_pmpi(real(MNRA - TRRA))
          end if
        enddo

        ! enforce zero mean EOT (just in case)
        if (ORBIT%EOT) then
          meanEOT = sum(ORBIT%ET)/DAYS_PER_CYCLE
          ORBIT%ET = ORBIT%ET - meanEOT
        end if

      end if

      if (present(FIX_SUN)) then
         ORBIT%FIX_SUN=FIX_SUN
      else
         ORBIT%FIX_SUN=.FALSE.
      end if

      ORBIT%CREATED = .TRUE.
      MAPL_SunOrbitCreate = ORBIT

      _RETURN(ESMF_SUCCESS)

      contains

         real(kind=REAL64) function dTRELdDAY(TREL,OMG,ECCENTRICITY,PRH)
            real(kind=REAL64), intent(in) :: TREL ! ecliptic longitude of true Sun
            real(kind=REAL64), intent(in) :: OMG
            real,              intent(in) :: ECCENTRICITY
            real(kind=REAL64), intent(in) :: PRH

            dTRELdDAY = OMG*(1.0-ECCENTRICITY*cos(TREL-PRH))**2
         end function dTRELdDAY

    end function MAPL_SunOrbitCreate

!==========================================================================
!>
! The function `MAPL_SunOrbitCreateFromConfig` works like `MAPL_SunOrbitCreate()`
!  but gets orbital parameters from Config CF.
!
   function MAPL_SunOrbitCreateFromConfig ( &
      CF, CLOCK, FIX_SUN, RC) result (ORBIT)

! !ARGUMENTS:

      type (ESMF_Config), intent(INOUT) :: CF
      type (ESMF_Clock),  intent(IN   ) :: CLOCK
      logical,            intent(IN   ) :: FIX_SUN
      integer, optional,  intent(OUT  ) :: RC

      type (MAPL_SunOrbit)              :: ORBIT

      character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitCreateFromConfig"
      integer :: STATUS

      real :: ECC, OB, PER
      integer :: EQNX

      logical :: EOT, ORBIT_ANAL2B
      integer :: ORB2B_REF_YYYYMMDD, ORB2B_REF_HHMMSS, &
           ORB2B_EQUINOX_YYYYMMDD, ORB2B_EQUINOX_HHMMSS
      real :: ORB2B_YEARLEN, &
           ORB2B_ECC_REF, ORB2B_ECC_RATE, &
           ORB2B_OBQ_REF, ORB2B_OBQ_RATE, &
           ORB2B_LAMBDAP_REF, ORB2B_LAMBDAP_RATE

      ! pmn: There is one orbit is per STATE, so, for example, the MAPL states of the
      ! solar and land gridded components can potentially have independent solar orbits.
      ! Usually these "independent orbits" will be IDENTICAL because the configuration
      ! resources such as "ECCENTRICITY:" or "EOT:" will not be qualified by the name
      ! of the gridded component. But for example, if the resource file specifies
      !   "EOT: .FALSE."
      ! but
      !   "SOLAR_EOT: .TRUE."
      ! then only SOLAR will have an EOT correction. The same goes for the new orbital
      ! system choice ORBIT_ANAL2B.
      !   A state's orbit is actually created in this routine by requesting the ORBIT
      ! object. If its not already created then it will be made below. GridComps that
      ! don't needed an orbit and dont request one will not have one.

      ! Parameters of standard orbital system (tabularized intercalation cycle)
      ! -----------------------------------------------------------------------
      call ESMF_ConfigGetAttribute (CF, &
         ECC, label="ECCENTRICITY:", &
         default=DEFAULT_ORBIT_ECCENTRICITY, _RC)

      call ESMF_ConfigGetAttribute (CF, &
         OB, label="OBLIQUITY:", &
         default=DEFAULT_ORBIT_OBLIQUITY, _RC)

      call ESMF_ConfigGetAttribute (CF, &
         PER, label="PERIHELION:", &
         default=DEFAULT_ORBIT_PERIHELION, _RC)

      call ESMF_ConfigGetAttribute (CF, &
         EQNX, label="EQUINOX:", &
         default=DEFAULT_ORBIT_EQUINOX, _RC)

      ! Apply Equation of Time correction?
      ! ----------------------------------
      call ESMF_ConfigGetAttribute (CF, &
         EOT, label="EOT:", &
         default=.FALSE., _RC)

      ! New orbital system (analytic two-body) allows some time-varying
      ! behavior, namely, linear variation in LAMBDAP, ECC, and OBQ.
      ! ---------------------------------------------------------------

      call ESMF_ConfigGetAttribute (CF, &
         ORBIT_ANAL2B, label="ORBIT_ANAL2B:", &
         default=.FALSE., _RC)

      ! Fixed anomalistic year length in mean solar days
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_YEARLEN, label="ORB2B_YEARLEN:", &
         default=DEFAULT_ORB2B_YEARLEN, _RC)

      ! Reference date and time for orbital parameters
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_REF_YYYYMMDD, label="ORB2B_REF_YYYYMMDD:", &
         default=DEFAULT_ORB2B_REF_YYYYMMDD, _RC)
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_REF_HHMMSS, label="ORB2B_REF_HHMMSS:", &
         default=DEFAULT_ORB2B_REF_HHMMSS, _RC)

      ! Orbital eccentricity at reference date
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_ECC_REF, label="ORB2B_ECC_REF:", &
         default=DEFAULT_ORB2B_ECC_REF, _RC)

      ! Rate of change of orbital eccentricity per Julian century
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_ECC_RATE, label="ORB2B_ECC_RATE:", &
         default=DEFAULT_ORB2B_ECC_RATE, _RC)

      ! Earth's obliquity (axial tilt) at reference date [degrees]
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_OBQ_REF, label="ORB2B_OBQ_REF:", &
         default=DEFAULT_ORB2B_OBQ_REF, _RC)

      ! Rate of change of obliquity [degrees per Julian century]
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_OBQ_RATE, label="ORB2B_OBQ_RATE:", &
         default=DEFAULT_ORB2B_OBQ_RATE, _RC)

      ! Longitude of perihelion at reference date [degrees]
      !   (from March equinox to perihelion in direction of earth's motion)
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_LAMBDAP_REF, label="ORB2B_LAMBDAP_REF:", &
         default=DEFAULT_ORB2B_LAMBDAP_REF, _RC)

      ! Rate of change of LAMBDAP [degrees per Julian century]
      !   (Combines both equatorial and ecliptic precession)
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_LAMBDAP_RATE, label="ORB2B_LAMBDAP_RATE:", &
         default=DEFAULT_ORB2B_LAMBDAP_RATE, _RC)

      ! March Equinox date and time
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_EQUINOX_YYYYMMDD, label="ORB2B_EQUINOX_YYYYMMDD:", &
         default=DEFAULT_ORB2B_EQUINOX_YYYYMMDD, _RC)
      call ESMF_ConfigGetAttribute (CF, &
         ORB2B_EQUINOX_HHMMSS, label="ORB2B_EQUINOX_HHMMSS:", &
         default=DEFAULT_ORB2B_EQUINOX_HHMMSS, _RC)

      ! create the orbit object
      ORBIT = MAPL_SunOrbitCreate ( &
         CLOCK, ECC, OB, PER, EQNX, &
         EOT, ORBIT_ANAL2B, ORB2B_YEARLEN, &
         ORB2B_REF_YYYYMMDD, ORB2B_REF_HHMMSS, &
         ORB2B_ECC_REF, ORB2B_ECC_RATE, &
         ORB2B_OBQ_REF, ORB2B_OBQ_RATE, &
         ORB2B_LAMBDAP_REF, ORB2B_LAMBDAP_RATE, &
         ORB2B_EQUINOX_YYYYMMDD, ORB2B_EQUINOX_HHMMSS, &
         FIX_SUN=FIX_SUN,_RC)

      _RETURN(ESMF_SUCCESS)

   end function MAPL_SunOrbitCreateFromConfig

!==========================================================================
!>
! The routine `MAPL_SunOrbitDestroy` Destroys a *GEOS_SunOrbit* object,
! deallocating the space used to save the ephemeris.
!
    subroutine MAPL_SunOrbitDestroy(ORBIT, RC)

! !ARGUMENTS:

       type(MAPL_SunOrbit),    intent(INOUT) :: ORBIT
       integer, optional,      intent(  OUT) :: RC

       character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitDestroy"

       if(associated(ORBIT%TH)) deallocate(ORBIT%TH)
       if(associated(ORBIT%ZC)) deallocate(ORBIT%ZC)
       if(associated(ORBIT%ZS)) deallocate(ORBIT%ZS)
       if(associated(ORBIT%PP)) deallocate(ORBIT%PP)
       if(associated(ORBIT%ET)) deallocate(ORBIT%ET)
       ORBIT%CREATED = .FALSE.

       _RETURN(ESMF_SUCCESS)

     end subroutine MAPL_SunOrbitDestroy



!==========================================================================
!>
! The function `MAPL_SunOrbitCreated` returns `.true.`
! if the given orbit object has been initilized.
!
       logical function  MAPL_SunOrbitCreated(ORBIT, RC)

! !ARGUMENTS:

       type(MAPL_SunOrbit),    intent(IN ) :: ORBIT
       integer, optional,      intent(OUT) :: RC

       character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitCreated"

       MAPL_SunOrbitCreated = ORBIT%CREATED
       _RETURN(ESMF_SUCCESS)
       return

     end function MAPL_SunOrbitCreated

!==========================================================================
!>
! Query for quantities in an orbit object.
! Optionally returns the parameters of the orbit and its
! associated `ESMF_Clock`. It fails
! if the orbit has not been created.
!
! @bug
! Not updated for ORBIT_ANAL2B option, which does not precalc
! many of the above outputs.
!@endbug
!
subroutine  MAPL_SunOrbitQuery(ORBIT,           &
                               ECCENTRICITY,    &
                               OBLIQUITY,       &
                               PERIHELION,      &
                               EQUINOX,         &
                               YEAR_LENGTH,     &
                               YEARS_PER_CYCLE, &
                               DAYS_PER_CYCLE,  &
                               CLOCK,           &
                               ZS,              &
                               ZC,              &
                               TH,              &
                               PP,              &
                               ET,              &
                                             RC )

! !ARGUMENTS:

       type(MAPL_SunOrbit),           intent(IN ) :: ORBIT
       real,                optional, intent(OUT) :: OBLIQUITY
       real,                optional, intent(OUT) :: ECCENTRICITY
       real,                optional, intent(OUT) :: PERIHELION
       real,                optional, intent(OUT) :: YEAR_LENGTH
       integer,             optional, intent(OUT) :: EQUINOX
       integer,             optional, intent(OUT) :: YEARS_PER_CYCLE
       integer,             optional, intent(OUT) :: DAYS_PER_CYCLE
       type(ESMF_Clock   ), optional, intent(OUT) :: CLOCK
       real,                optional, pointer, dimension(:) :: ZS
       real,                optional, pointer, dimension(:) :: ZC
       real,                optional, pointer, dimension(:) :: TH
       real,                optional, pointer, dimension(:) :: PP
       real,                optional, pointer, dimension(:) :: ET
       integer,             optional, intent(OUT) :: RC

       character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitQuery"
       integer :: STATUS

       _ASSERT(MAPL_SunOrbitCreated(ORBIT,RC=STATUS),'MAPL_SunOrbit not yet created!')

       if(present(CLOCK          )) CLOCK           = ORBIT%CLOCK
       if(present(OBLIQUITY      )) OBLIQUITY       = ORBIT%OB
       if(present(ECCENTRICITY   )) ECCENTRICITY    = ORBIT%ECC
       if(present(PERIHELION     )) PERIHELION      = ORBIT%PER
       if(present(EQUINOX        )) EQUINOX         = ORBIT%EQNX
       if(present(YEAR_LENGTH    )) YEAR_LENGTH     = ORBIT%YEARLEN
       if(present(DAYS_PER_CYCLE )) DAYS_PER_CYCLE  = ORBIT%DAYS_PER_CYCLE
       if(present(YEARS_PER_CYCLE)) YEARS_PER_CYCLE = ORBIT%YEARS_PER_CYCLE
       if(present(ZS             )) ZS => ORBIT%ZS
       if(present(ZC             )) ZC => ORBIT%ZC
       if(present(TH             )) TH => ORBIT%TH
       if(present(PP             )) PP => ORBIT%PP
       if(present(ET             )) ET => ORBIT%ET

       _RETURN(ESMF_SUCCESS)

     end subroutine MAPL_SunOrbitQuery

!==========================================================================
!>
! `GEOS_SunGetInsolation` returns the cosine of the solar zenith angle and the
! insolation at the top of the atmosphere for the given reference time, latitudes,
! longitudes, and orbit.  It is overloaded to accept either 1d or 2d
! FORTRAN arrays or ESMF arrays of lats and lons and to produce the
! corresponding outputs.
!
! The reference time is obtained as follows. If CurrTime is specified, it is used;
! otherwise, if the optional clock is specified, it is set to the time
! on that clock. If neither currTime nor clock are given, the time on the attached
! clock is used.
!
! If the optional time interval is specified, the return values are averages
! over that interval following the reference time. In this case, the cosine of
! the solar zenith angle is an insolation-weighted average. The straight average
! of the zenith angle and the average over the daylight part of the
! interval (ZTHB,ZTHD) and the values at the beginning and end of the interval
! (ZTH1,ZTHN) are also optionally available. The last two are supported only for
! first two (non-ESMF) overloads.
! PMN Jun 2020: Added optional ZTHP, like ZTHB a pure average of the cosine of the
! solar zenith angle, except that for ZTHP the value averaged is allowed to be
! negative (below the horizon). It will be used for Photolysis calculations.
! Not implemented for ESMF overload.
!
! If the interval is not specified, the values are instantaneous values valid at
! the reference time.
!
! The optional *TIME* argument is used to return some specialized
! insolations. For example, the orbit at any of four Equinox or Solstice
! positions. If *TIME* is present, only the time of day is used from the clock,
! and a time interval, if specified, must be less than 24 hours. It can also be
! used to return daily-mean insolation for the date on the clock, or the annual-mean
! insolation for the year on the clock.
!
! The *TIME* argument can be any of the following:
!```
!      MAPL_SunAutumnalEquinox
!      MAPL_SunWinterSolstice
!      MAPL_SunVernalEquinox
!      MAPL_SunSummerSolstice
!      MAPL_SunDailyMean
!      MAPL_SunAnnualMean
!```
!
! @note
! If ORBIT%EOT is .TRUE., an Equation of Time correction will be
! applied. This shifts the actual daylight period w.r.t. to mean solar
! noon, to account for small but cumulative eccentricity and obliquity
! effects on the actual length of the solar day.
!@endnote
!
#undef  DIMENSIONS
#define DIMENSIONS (:)
#define THE_SIZE   (size(LONS,1))
      recursive subroutine SOLAR_1D(LONS, LATS, ORBIT,ZTH,SLR,INTV,CLOCK, &
                                    TIME,currTime,DIST,ZTHB,ZTHD,ZTH1,ZTHN,ZTHP,&
                                    STEPSIZE,RC)
#include "sun.H"
      end subroutine SOLAR_1D

!==========================================================================

#undef  DIMENSIONS
#undef  THE_SIZE
#define DIMENSIONS (:,:)
#define THE_SIZE   (size(LONS,1),size(LONS,2))
      recursive subroutine SOLAR_2D(LONS, LATS, ORBIT,ZTH,SLR,INTV,CLOCK, &
                                    TIME,currTime,DIST,ZTHB,ZTHD,ZTH1,ZTHN,ZTHP,&
                                    STEPSIZE,RC)
#include "sun.H"
      end subroutine SOLAR_2D
#undef  DIMENSIONS
#undef  THE_SIZE

!BOP



!EOP


!==========================================================================

      subroutine SOLAR_ARR_INT(LONS, LATS, ORBIT, ZTH, SLR, INTV, CLOCK,  &
                               TIME, currTime, DIST, ZTHB, ZTHD,  RC)

      type (MAPL_SunOrbit),               intent(IN ) :: ORBIT
      type (ESMF_Array),                  intent(IN ) :: LATS
      type (ESMF_Array),                  intent(IN ) :: LONS
      type (ESMF_Array),                  intent(OUT) :: ZTH
      type (ESMF_Array),                  intent(OUT) :: SLR
      type (ESMF_TimeInterval), optional, intent(INout) :: INTV
      type (ESMF_Clock),        optional, intent(IN ) :: CLOCK
      type (ESMF_Time),         optional, intent(IN ) :: currTime
      integer,                  optional, intent(IN ) :: TIME
      real,                     optional, intent(OUT) :: DIST
      type (ESMF_Array),        optional, intent(OUT) :: ZTHB
      type (ESMF_Array),        optional, intent(OUT) :: ZTHD
      integer,                  optional, intent(OUT) :: RC

!   Locals


!!$      character(len=ESMF_MAXSTR)      :: IAm = "SunGetInsolationArr"
      integer                         :: STATUS

      real, pointer, dimension (:  )  :: LONS1, LATS1, ZTH1, SLR1, ZTHB1, ZTHD1
      real, pointer, dimension (:,:)  :: LONS2, LATS2, ZTH2, SLR2, ZTHB2, ZTHD2

      integer                         :: RANK

!   Begin

      _FAIL(" pmn: this routine is not up to date, is it even used anywhere?")

      call ESMF_ArrayGet(LONS, RANK=RANK, RC=STATUS)
      _VERIFY(STATUS)

      select case(RANK)

      case(1)
         call ESMF_ArrayGet(LATS, localDE=0, farrayptr=LATS1, RC=STATUS)
         _VERIFY(STATUS)
         call ESMF_ArrayGet(LONS,localDE=0, farrayptr=LONS1, RC=STATUS)
         _VERIFY(STATUS)
         call ESMF_ArrayGet(ZTH ,localDE=0, farrayptr=ZTH1, RC=STATUS)
         _VERIFY(STATUS)
         call ESMF_ArrayGet(SLR ,localDE=0, farrayptr=SLR1, RC=STATUS)
         _VERIFY(STATUS)

         if(present(ZTHB) .and. present(ZTHD)) then
            call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB1 ,RC=STATUS)
            _VERIFY(STATUS)
            call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD1 ,RC=STATUS)
            _VERIFY(STATUS)
            call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,ZTHB=ZTHB1,ZTHD=ZTHD1,RC=STATUS)
         elseif(present(ZTHB)) then
            call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB1 ,RC=STATUS)
            _VERIFY(STATUS)
            call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,ZTHB=ZTHB1,RC=STATUS)
         elseif(present(ZTHD)) then
            call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD1 ,RC=STATUS)
            _VERIFY(STATUS)
            call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,ZTHD=ZTHD1,RC=STATUS)
         else
            call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,RC=STATUS)
         endif
         _VERIFY(STATUS)

      case(2)
         call ESMF_ArrayGet(LATS,localDE=0, farrayptr=LATS2,RC=STATUS)
         _VERIFY(STATUS)
         call ESMF_ArrayGet(LONS,localDE=0, farrayptr=LONS2,RC=STATUS)
         _VERIFY(STATUS)
         call ESMF_ArrayGet(ZTH ,localDE=0, farrayptr=ZTH2 ,RC=STATUS)
         _VERIFY(STATUS)
         call ESMF_ArrayGet(SLR ,localDE=0, farrayptr=SLR2 ,RC=STATUS)
         _VERIFY(STATUS)

         if(present(ZTHB) .and. present(ZTHD)) then
            call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB2 ,RC=STATUS)
            _VERIFY(STATUS)
            call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD2 ,RC=STATUS)
            _VERIFY(STATUS)
            call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,ZTHB=ZTHB2,ZTHD=ZTHD2,RC=STATUS)
         elseif(present(ZTHB)) then
            call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB2 ,RC=STATUS)
            _VERIFY(STATUS)
            call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,ZTHB=ZTHB2,RC=STATUS)
         elseif(present(ZTHD)) then
            call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD2 ,RC=STATUS)
            _VERIFY(STATUS)
            call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,ZTHD=ZTHD2,RC=STATUS)
         else
            call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,&
                 TIME,currTime,DIST=DIST,RC=STATUS)
         endif
         _VERIFY(STATUS)

      case default
         _RETURN(ESMF_FAILURE)

      end select

      _RETURN(ESMF_SUCCESS)

    end subroutine SOLAR_ARR_INT

    subroutine GETIDAY(IDAY,TIME,ORBIT,RC)
      integer,              intent(OUT) :: IDAY
      integer,              intent(IN ) :: TIME
      type(MAPL_SunORBIT),  intent(IN ) :: ORBIT
      integer, optional,    intent(OUT) :: RC

      real :: ANOMALY

      select case(TIME)
      case (MAPL_SunAutumnalEquinox)
         ANOMALY = MAPL_PI
      case (MAPL_SunWinterSolstice )
         ANOMALY = (MAPL_PI*3.0)/2.0
      case (MAPL_SunVernalEquinox  )
         ANOMALY = 0.0
      case (MAPL_SunSummerSolstice )
         ANOMALY = MAPL_PI/2.0
      case  default
         _RETURN(ESMF_FAILURE)
      end select

      do IDAY=1,ORBIT%DAYS_PER_CYCLE-1
         if(ORBIT%TH(IDAY)<=ANOMALY .and. ORBIT%TH(IDAY+1)>ANOMALY) then
            _RETURN(ESMF_SUCCESS)
         end if
      end do

      _RETURN(ESMF_FAILURE)
    end subroutine GETIDAY
!==========================================================================

    subroutine MAPL_SunGetSolarConstantByTime(Time,SC,HK,rc)

      type (ESMF_Time),  intent(INOUT) :: Time
      real,              intent(OUT)   :: SC
      real,   optional,  intent(OUT)   :: HK(:)
      integer, optional, intent(OUT)   :: rc

      integer    :: YY, DOY
      integer    :: STATUS

      call ESMF_TimeGet (TIME, YY=YY, DayOfYear=DOY, RC=STATUS)
      _VERIFY(STATUS)

      call MAPL_SunGetSolarConstantByYearDoY(YY,DOY,SC,HK, RC=STATUS)
      _VERIFY(STATUS)

      _RETURN(ESMF_SUCCESS)
    end subroutine MAPL_SunGetSolarConstantByTime

!-------------------------------------------------------------------------
!NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1, GEOS/DAS!
!-------------------------------------------------------------------------
!>
! Given the year and day-of-year, this function returns the CMIP5 solar constant
! and the 8 band fractions for the Chou solar code. These are based on annual
! values from 1610 to 2008 and a repeating cycle after 2008. For dates prior
! tp 1 January 1610, it returns the value at the start of 1610.
!
! DayOfYear is expected to have a value of 1.00 at 0:00 UTC Jan 1 and leap years
! are acounted for by repeating the 366th day on January first of the following year.
!
! The SC values have been multiplied by .9965 to calibrate to the SOURCE/TIM scale.
! The SC values include the "background" variation.
!
! @bug
! Band values for RRTMG not implemented
!@endbug
!
   subroutine MAPL_SunGetSolarConstantByYearDoY(year,dayofyear,SC,HK, rc)

   integer,           intent(IN)  :: Year
   integer,           intent(IN)  :: DayOfYear
   real,              intent(OUT) :: SC
   real,    optional, intent(OUT) :: HK(:)
   integer, optional, intent(OUT) :: rc

   real    :: F
   integer :: i1,i2,Current

   integer, parameter :: firstYear = 1610
   integer, parameter :: finalYear = 2008
   integer, parameter :: SolCycle  = 12  ! Cycle 23 was 12.2 years

   real, parameter ::  SO(firstYear:finalYear) = (/  &
         1360.768, 1360.751, 1361.147, 1361.259, 1361.333,&
         1361.051, 1360.636, 1360.510, 1360.500, 1360.590,&
         1360.586, 1360.582, 1360.578, 1360.574, 1360.534,&
         1360.751, 1360.640, 1360.563, 1360.573, 1360.514,&
         1360.365, 1360.373, 1360.340, 1360.438, 1360.358,&
         1360.363, 1360.330, 1360.327, 1360.792, 1360.844,&
         1360.415, 1360.304, 1360.622, 1360.404, 1360.349,&
         1360.268, 1360.267, 1360.267, 1360.265, 1360.262,&
         1360.261, 1360.261, 1360.288, 1360.266, 1360.258,&
         1360.248, 1360.247, 1360.242, 1360.237, 1360.235,&
         1360.248, 1360.240, 1360.235, 1360.235, 1360.235,&
         1360.235, 1360.235, 1360.235, 1360.235, 1360.235,&
         1360.235, 1360.241, 1360.237, 1360.235, 1360.236,&
         1360.234, 1360.246, 1360.236, 1360.235, 1360.234,&
         1360.239, 1360.234, 1360.234, 1360.234, 1360.244,&
         1360.234, 1360.238, 1360.235, 1360.238, 1360.237,&
         1360.237, 1360.238, 1360.240, 1360.240, 1360.240,&
         1360.241, 1360.240, 1360.240, 1360.240, 1360.241,&
         1360.244, 1360.246, 1360.249, 1360.267, 1360.282,&
         1360.297, 1360.286, 1360.302, 1360.288, 1360.284,&
         1360.286, 1360.296, 1360.302, 1360.308, 1360.325,&
         1360.352, 1360.391, 1360.451, 1360.393, 1360.564,&
         1360.507, 1360.510, 1360.443, 1360.406, 1360.490,&
         1360.479, 1360.648, 1360.654, 1360.844, 1360.569,&
         1360.881, 1360.403, 1360.530, 1360.420, 1360.430,&
         1360.560, 1360.769, 1360.604, 1360.558, 1360.800,&
         1360.503, 1360.831, 1360.546, 1360.491, 1360.436,&
         1360.438, 1360.444, 1360.452, 1360.875, 1360.895,&
         1360.859, 1360.689, 1360.661, 1360.643, 1360.556,&
         1360.543, 1360.571, 1360.694, 1360.807, 1360.870,&
         1360.846, 1361.010, 1360.872, 1360.802, 1360.778,&
         1360.628, 1360.597, 1360.802, 1361.054, 1361.238,&
         1361.241, 1361.135, 1361.063, 1360.845, 1360.807,&
         1360.670, 1360.725, 1360.886, 1361.115, 1361.173,&
         1360.995, 1361.106, 1360.842, 1360.749, 1360.618,&
         1360.677, 1360.988, 1361.153, 1361.100, 1361.079,&
         1360.980, 1360.830, 1360.817, 1360.796, 1360.702,&
         1360.586, 1360.556, 1360.503, 1360.478, 1360.485,&
         1360.525, 1360.796, 1360.675, 1360.543, 1360.547,&
         1360.557, 1360.460, 1360.393, 1360.368, 1360.340,&
         1360.332, 1360.344, 1360.385, 1360.434, 1360.457,&
         1360.512, 1360.608, 1360.581, 1360.533, 1360.518,&
         1360.463, 1360.441, 1360.465, 1360.470, 1360.511,&
         1360.584, 1360.687, 1360.798, 1360.885, 1360.903,&
         1360.943, 1360.779, 1360.678, 1360.582, 1360.618,&
         1360.878, 1361.258, 1361.331, 1361.119, 1361.037,&
         1360.921, 1360.770, 1360.715, 1360.645, 1360.675,&
         1360.803, 1360.900, 1361.001, 1361.197, 1361.188,&
         1360.990, 1360.994, 1360.915, 1360.805, 1360.682,&
         1360.604, 1360.606, 1360.715, 1360.908, 1361.125,&
         1361.191, 1361.082, 1360.932, 1360.849, 1360.787,&
         1360.696, 1360.633, 1360.581, 1360.740, 1360.976,&
         1361.213, 1361.154, 1361.107, 1360.907, 1360.779,&
         1360.631, 1360.581, 1360.581, 1360.552, 1360.575,&
         1360.721, 1360.865, 1360.888, 1360.935, 1360.837,&
         1360.728, 1360.563, 1360.520, 1360.484, 1360.456,&
         1360.491, 1360.741, 1360.839, 1360.980, 1361.075,&
         1360.959, 1360.779, 1360.634, 1360.611, 1360.559,&
         1360.529, 1360.451, 1360.459, 1360.669, 1360.938,&
         1360.750, 1360.946, 1360.830, 1360.895, 1360.785,&
         1360.652, 1360.569, 1360.523, 1360.539, 1360.605,&
         1360.909, 1361.118, 1361.267, 1361.229, 1361.000,&
         1360.838, 1360.733, 1360.616, 1360.648, 1360.679,&
         1360.882, 1360.983, 1361.166, 1361.044, 1361.003,&
         1360.985, 1360.864, 1360.757, 1360.637, 1360.748,&
         1360.964, 1361.352, 1361.286, 1361.222, 1361.206,&
         1361.143, 1361.065, 1360.962, 1360.805, 1360.834,&
         1361.100, 1361.198, 1361.437, 1361.567, 1361.474,&
         1361.232, 1360.996, 1360.987, 1360.849, 1360.877,&
         1360.997, 1361.529, 1361.885, 1361.850, 1361.601,&
         1361.495, 1361.139, 1360.968, 1360.916, 1360.918,&
         1360.954, 1361.137, 1361.333, 1361.383, 1361.466,&
         1361.461, 1361.177, 1361.271, 1361.019, 1360.947,&
         1360.755, 1360.866, 1361.053, 1361.493, 1361.852,&
         1361.865, 1361.912, 1361.504, 1361.418, 1361.030,&
         1360.862, 1360.858, 1361.010, 1361.301, 1361.865,&
         1361.770, 1361.663, 1361.520, 1361.247, 1361.017,&
         1360.920, 1360.832, 1360.960, 1361.321, 1361.603,&
         1361.882, 1361.819, 1361.897, 1361.448, 1361.267,&
         1361.074, 1361.030, 1360.944, 1360.912 /)

   real, parameter :: CHOUBAND1(firstYear:finalYear) = (/ &
         0.00516,   0.00515,   0.00519,   0.00519,   0.00520,&
         0.00518,   0.00515,   0.00514,   0.00514,   0.00514,&
         0.00514,   0.00514,   0.00514,   0.00514,   0.00514,&
         0.00516,   0.00515,   0.00514,   0.00514,   0.00514,&
         0.00513,   0.00513,   0.00513,   0.00514,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00517,   0.00517,&
         0.00514,   0.00513,   0.00515,   0.00514,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00514,   0.00513,   0.00515,&
         0.00514,   0.00514,   0.00514,   0.00513,   0.00514,&
         0.00514,   0.00515,   0.00515,   0.00517,   0.00514,&
         0.00517,   0.00513,   0.00514,   0.00513,   0.00513,&
         0.00514,   0.00516,   0.00515,   0.00514,   0.00516,&
         0.00514,   0.00516,   0.00514,   0.00514,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00517,   0.00517,&
         0.00517,   0.00515,   0.00515,   0.00515,   0.00514,&
         0.00514,   0.00514,   0.00515,   0.00516,   0.00516,&
         0.00516,   0.00517,   0.00516,   0.00516,   0.00515,&
         0.00514,   0.00514,   0.00516,   0.00518,   0.00519,&
         0.00519,   0.00518,   0.00518,   0.00516,   0.00515,&
         0.00514,   0.00515,   0.00516,   0.00518,   0.00518,&
         0.00517,   0.00518,   0.00516,   0.00515,   0.00514,&
         0.00515,   0.00517,   0.00518,   0.00518,   0.00518,&
         0.00517,   0.00516,   0.00516,   0.00516,   0.00515,&
         0.00514,   0.00514,   0.00514,   0.00514,   0.00514,&
         0.00514,   0.00516,   0.00515,   0.00514,   0.00514,&
         0.00514,   0.00514,   0.00513,   0.00513,   0.00513,&
         0.00513,   0.00513,   0.00513,   0.00514,   0.00514,&
         0.00514,   0.00515,   0.00515,   0.00514,   0.00514,&
         0.00514,   0.00513,   0.00514,   0.00513,   0.00514,&
         0.00514,   0.00515,   0.00516,   0.00517,   0.00517,&
         0.00517,   0.00516,   0.00515,   0.00514,   0.00514,&
         0.00516,   0.00519,   0.00520,   0.00518,   0.00517,&
         0.00517,   0.00515,   0.00515,   0.00514,   0.00515,&
         0.00516,   0.00516,   0.00517,   0.00519,   0.00518,&
         0.00517,   0.00517,   0.00517,   0.00516,   0.00515,&
         0.00514,   0.00514,   0.00515,   0.00516,   0.00518,&
         0.00519,   0.00518,   0.00517,   0.00516,   0.00516,&
         0.00515,   0.00514,   0.00514,   0.00515,   0.00517,&
         0.00519,   0.00518,   0.00518,   0.00517,   0.00516,&
         0.00515,   0.00514,   0.00514,   0.00514,   0.00514,&
         0.00515,   0.00516,   0.00516,   0.00517,   0.00517,&
         0.00516,   0.00515,   0.00514,   0.00514,   0.00514,&
         0.00514,   0.00516,   0.00517,   0.00518,   0.00518,&
         0.00517,   0.00516,   0.00515,   0.00515,   0.00514,&
         0.00514,   0.00514,   0.00514,   0.00515,   0.00516,&
         0.00517,   0.00517,   0.00517,   0.00517,   0.00516,&
         0.00515,   0.00514,   0.00514,   0.00514,   0.00515,&
         0.00517,   0.00518,   0.00520,   0.00519,   0.00518,&
         0.00516,   0.00516,   0.00515,   0.00515,   0.00515,&
         0.00517,   0.00518,   0.00518,   0.00519,   0.00518,&
         0.00517,   0.00516,   0.00515,   0.00515,   0.00515,&
         0.00517,   0.00519,   0.00521,   0.00520,   0.00520,&
         0.00518,   0.00517,   0.00517,   0.00516,   0.00516,&
         0.00517,   0.00520,   0.00522,   0.00522,   0.00522,&
         0.00519,   0.00518,   0.00517,   0.00516,   0.00516,&
         0.00517,   0.00522,   0.00525,   0.00525,   0.00523,&
         0.00521,   0.00517,   0.00517,   0.00516,   0.00516,&
         0.00516,   0.00517,   0.00520,   0.00520,   0.00520,&
         0.00521,   0.00518,   0.00519,   0.00517,   0.00516,&
         0.00515,   0.00516,   0.00517,   0.00520,   0.00523,&
         0.00523,   0.00523,   0.00522,   0.00519,   0.00517,&
         0.00516,   0.00516,   0.00516,   0.00519,   0.00523,&
         0.00522,   0.00522,   0.00520,   0.00518,   0.00517,&
         0.00516,   0.00515,   0.00516,   0.00518,   0.00520,&
         0.00522,   0.00522,   0.00522,   0.00519,   0.00518,&
         0.00517,   0.00516,   0.00516,   0.00516     /)

   real, parameter :: CHOUBAND2(firstYear:finalYear) = (/ &
         0.00567,   0.00567,   0.00569,   0.00570,   0.00570,&
         0.00569,   0.00567,   0.00566,   0.00566,   0.00567,&
         0.00567,   0.00567,   0.00567,   0.00567,   0.00566,&
         0.00567,   0.00567,   0.00567,   0.00567,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00568,   0.00568,&
         0.00566,   0.00566,   0.00567,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00567,&
         0.00567,   0.00567,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00567,   0.00567,   0.00568,   0.00567,&
         0.00568,   0.00566,   0.00567,   0.00566,   0.00566,&
         0.00567,   0.00568,   0.00567,   0.00567,   0.00568,&
         0.00566,   0.00568,   0.00567,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00568,   0.00568,&
         0.00568,   0.00567,   0.00567,   0.00567,   0.00566,&
         0.00566,   0.00566,   0.00567,   0.00568,   0.00568,&
         0.00568,   0.00568,   0.00568,   0.00567,   0.00567,&
         0.00567,   0.00566,   0.00567,   0.00569,   0.00569,&
         0.00569,   0.00569,   0.00569,   0.00568,   0.00567,&
         0.00567,   0.00567,   0.00568,   0.00569,   0.00569,&
         0.00568,   0.00569,   0.00568,   0.00567,   0.00567,&
         0.00567,   0.00568,   0.00569,   0.00569,   0.00569,&
         0.00568,   0.00568,   0.00568,   0.00567,   0.00567,&
         0.00567,   0.00567,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00568,   0.00567,   0.00567,   0.00567,&
         0.00567,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00566,   0.00567,   0.00567,   0.00567,   0.00566,&
         0.00566,   0.00566,   0.00566,   0.00566,   0.00566,&
         0.00567,   0.00567,   0.00568,   0.00568,   0.00568,&
         0.00568,   0.00567,   0.00567,   0.00566,   0.00567,&
         0.00568,   0.00569,   0.00570,   0.00569,   0.00568,&
         0.00568,   0.00567,   0.00567,   0.00567,   0.00567,&
         0.00567,   0.00568,   0.00568,   0.00569,   0.00569,&
         0.00568,   0.00568,   0.00568,   0.00567,   0.00567,&
         0.00567,   0.00567,   0.00567,   0.00568,   0.00569,&
         0.00569,   0.00569,   0.00568,   0.00568,   0.00567,&
         0.00567,   0.00567,   0.00566,   0.00567,   0.00568,&
         0.00569,   0.00569,   0.00569,   0.00568,   0.00567,&
         0.00567,   0.00567,   0.00567,   0.00566,   0.00566,&
         0.00567,   0.00568,   0.00568,   0.00568,   0.00568,&
         0.00568,   0.00567,   0.00567,   0.00566,   0.00566,&
         0.00567,   0.00568,   0.00568,   0.00569,   0.00569,&
         0.00568,   0.00568,   0.00567,   0.00567,   0.00567,&
         0.00567,   0.00566,   0.00566,   0.00567,   0.00568,&
         0.00568,   0.00568,   0.00568,   0.00568,   0.00568,&
         0.00567,   0.00567,   0.00567,   0.00567,   0.00567,&
         0.00568,   0.00569,   0.00570,   0.00569,   0.00569,&
         0.00568,   0.00567,   0.00567,   0.00567,   0.00567,&
         0.00568,   0.00569,   0.00569,   0.00569,   0.00569,&
         0.00568,   0.00568,   0.00567,   0.00567,   0.00567,&
         0.00568,   0.00570,   0.00570,   0.00570,   0.00570,&
         0.00569,   0.00568,   0.00568,   0.00568,   0.00567,&
         0.00568,   0.00570,   0.00571,   0.00571,   0.00571,&
         0.00569,   0.00569,   0.00568,   0.00568,   0.00567,&
         0.00568,   0.00571,   0.00573,   0.00573,   0.00572,&
         0.00570,   0.00569,   0.00568,   0.00568,   0.00568,&
         0.00568,   0.00569,   0.00570,   0.00570,   0.00570,&
         0.00570,   0.00569,   0.00569,   0.00568,   0.00568,&
         0.00567,   0.00568,   0.00568,   0.00570,   0.00572,&
         0.00572,   0.00572,   0.00571,   0.00570,   0.00568,&
         0.00568,   0.00567,   0.00568,   0.00570,   0.00572,&
         0.00571,   0.00571,   0.00570,   0.00569,   0.00568,&
         0.00568,   0.00567,   0.00568,   0.00569,   0.00570,&
         0.00571,   0.00571,   0.00572,   0.00570,   0.00569,&
         0.00568,   0.00568,   0.00568,   0.00567    /)

   real, parameter :: CHOUBAND3(firstYear:finalYear) =  (/ &
         0.01202,   0.01202,   0.01203,   0.01203,   0.01203,&
         0.01202,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01202,   0.01202,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01202,   0.01202,&
         0.01201,   0.01201,   0.01202,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01202,   0.01202,   0.01202,   0.01201,&
         0.01202,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01202,   0.01201,   0.01201,   0.01202,&
         0.01201,   0.01202,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01202,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01201,&
         0.01201,   0.01201,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01201,   0.01201,   0.01202,   0.01202,   0.01203,&
         0.01203,   0.01203,   0.01202,   0.01202,   0.01202,&
         0.01201,   0.01202,   0.01202,   0.01203,   0.01203,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01201,&
         0.01202,   0.01202,   0.01203,   0.01203,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01202,   0.01202,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01201,   0.01201,&
         0.01202,   0.01203,   0.01203,   0.01203,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01201,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01203,   0.01203,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01201,   0.01201,   0.01202,   0.01202,   0.01203,&
         0.01203,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01201,   0.01201,   0.01202,   0.01202,&
         0.01203,   0.01203,   0.01203,   0.01202,   0.01202,&
         0.01201,   0.01201,   0.01201,   0.01201,   0.01201,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01202,   0.01201,   0.01201,   0.01201,&
         0.01201,   0.01202,   0.01202,   0.01203,   0.01203,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01201,&
         0.01201,   0.01201,   0.01201,   0.01202,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01201,   0.01201,   0.01201,   0.01202,&
         0.01202,   0.01203,   0.01203,   0.01203,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01202,   0.01203,   0.01203,   0.01203,&
         0.01202,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01203,   0.01203,   0.01203,   0.01203,&
         0.01203,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01203,   0.01204,   0.01204,   0.01204,&
         0.01203,   0.01202,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01204,   0.01205,   0.01205,   0.01204,&
         0.01203,   0.01203,   0.01202,   0.01202,   0.01202,&
         0.01202,   0.01203,   0.01203,   0.01203,   0.01203,&
         0.01203,   0.01203,   0.01203,   0.01202,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01203,   0.01204,&
         0.01204,   0.01204,   0.01204,   0.01203,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01203,   0.01204,&
         0.01204,   0.01204,   0.01203,   0.01203,   0.01202,&
         0.01202,   0.01202,   0.01202,   0.01203,   0.01203,&
         0.01204,   0.01204,   0.01204,   0.01203,   0.01203,&
         0.01202,   0.01202,   0.01202,   0.01202     /)

   real, parameter :: CHOUBAND4(firstYear:finalYear) = (/ &
         0.05687,   0.05687,   0.05690,   0.05690,   0.05691,&
         0.05689,   0.05686,   0.05685,   0.05685,   0.05686,&
         0.05686,   0.05686,   0.05686,   0.05686,   0.05685,&
         0.05687,   0.05686,   0.05686,   0.05686,   0.05685,&
         0.05684,   0.05685,   0.05684,   0.05685,   0.05684,&
         0.05685,   0.05684,   0.05684,   0.05687,   0.05688,&
         0.05685,   0.05684,   0.05686,   0.05685,   0.05685,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05684,   0.05684,   0.05684,   0.05684,&
         0.05684,   0.05685,   0.05685,   0.05685,   0.05686,&
         0.05685,   0.05685,   0.05685,   0.05685,   0.05685,&
         0.05685,   0.05686,   0.05686,   0.05688,   0.05686,&
         0.05688,   0.05685,   0.05686,   0.05685,   0.05685,&
         0.05686,   0.05687,   0.05686,   0.05686,   0.05687,&
         0.05685,   0.05688,   0.05686,   0.05685,   0.05685,&
         0.05685,   0.05685,   0.05685,   0.05688,   0.05688,&
         0.05688,   0.05687,   0.05686,   0.05686,   0.05686,&
         0.05685,   0.05686,   0.05686,   0.05687,   0.05688,&
         0.05687,   0.05689,   0.05688,   0.05687,   0.05687,&
         0.05686,   0.05686,   0.05687,   0.05689,   0.05690,&
         0.05690,   0.05689,   0.05689,   0.05687,   0.05687,&
         0.05686,   0.05686,   0.05688,   0.05689,   0.05690,&
         0.05688,   0.05689,   0.05687,   0.05687,   0.05686,&
         0.05686,   0.05688,   0.05690,   0.05689,   0.05689,&
         0.05688,   0.05687,   0.05687,   0.05687,   0.05687,&
         0.05686,   0.05686,   0.05685,   0.05685,   0.05685,&
         0.05685,   0.05687,   0.05686,   0.05686,   0.05686,&
         0.05686,   0.05685,   0.05685,   0.05685,   0.05684,&
         0.05684,   0.05684,   0.05685,   0.05685,   0.05685,&
         0.05685,   0.05686,   0.05686,   0.05686,   0.05685,&
         0.05685,   0.05685,   0.05685,   0.05685,   0.05685,&
         0.05686,   0.05686,   0.05687,   0.05688,   0.05688,&
         0.05688,   0.05687,   0.05686,   0.05686,   0.05686,&
         0.05688,   0.05690,   0.05691,   0.05689,   0.05689,&
         0.05688,   0.05687,   0.05686,   0.05686,   0.05686,&
         0.05687,   0.05688,   0.05688,   0.05690,   0.05690,&
         0.05688,   0.05688,   0.05688,   0.05687,   0.05686,&
         0.05686,   0.05686,   0.05686,   0.05688,   0.05689,&
         0.05690,   0.05689,   0.05688,   0.05687,   0.05687,&
         0.05686,   0.05686,   0.05686,   0.05687,   0.05688,&
         0.05690,   0.05690,   0.05689,   0.05688,   0.05687,&
         0.05686,   0.05686,   0.05686,   0.05685,   0.05686,&
         0.05687,   0.05688,   0.05688,   0.05688,   0.05688,&
         0.05687,   0.05686,   0.05686,   0.05686,   0.05685,&
         0.05686,   0.05687,   0.05688,   0.05689,   0.05690,&
         0.05689,   0.05687,   0.05687,   0.05687,   0.05686,&
         0.05686,   0.05685,   0.05686,   0.05687,   0.05688,&
         0.05688,   0.05689,   0.05688,   0.05688,   0.05688,&
         0.05687,   0.05686,   0.05686,   0.05686,   0.05686,&
         0.05688,   0.05689,   0.05691,   0.05690,   0.05689,&
         0.05688,   0.05687,   0.05686,   0.05686,   0.05687,&
         0.05688,   0.05689,   0.05690,   0.05690,   0.05689,&
         0.05688,   0.05688,   0.05687,   0.05686,   0.05687,&
         0.05688,   0.05691,   0.05691,   0.05691,   0.05691,&
         0.05690,   0.05689,   0.05688,   0.05687,   0.05687,&
         0.05689,   0.05691,   0.05693,   0.05693,   0.05692,&
         0.05690,   0.05689,   0.05688,   0.05688,   0.05688,&
         0.05689,   0.05693,   0.05695,   0.05695,   0.05694,&
         0.05692,   0.05689,   0.05688,   0.05688,   0.05688,&
         0.05688,   0.05689,   0.05691,   0.05691,   0.05692,&
         0.05692,   0.05690,   0.05690,   0.05689,   0.05688,&
         0.05687,   0.05688,   0.05689,   0.05692,   0.05694,&
         0.05694,   0.05695,   0.05693,   0.05691,   0.05689,&
         0.05688,   0.05687,   0.05688,   0.05691,   0.05694,&
         0.05694,   0.05693,   0.05692,   0.05690,   0.05688,&
         0.05688,   0.05687,   0.05688,   0.05690,   0.05692,&
         0.05694,   0.05694,   0.05694,   0.05691,   0.05690,&
         0.05689,   0.05689,   0.05688,   0.05688    /)

   real, parameter :: CHOUBAND5(firstYear:finalYear) = (/ &
         0.37831,   0.37831,   0.37831,   0.37832,   0.37832,&
         0.37831,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37831,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37829,   0.37829,   0.37829,   0.37830,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37830,   0.37830,&
         0.37829,   0.37829,   0.37830,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37829,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37830,   0.37829,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37831,   0.37830,&
         0.37831,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37831,   0.37830,   0.37830,   0.37831,&
         0.37830,   0.37831,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37831,   0.37831,&
         0.37831,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37831,   0.37832,&
         0.37832,   0.37832,   0.37832,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37832,   0.37832,&
         0.37831,   0.37832,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37832,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37831,   0.37831,&
         0.37830,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37831,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37829,   0.37829,&
         0.37829,   0.37829,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37830,   0.37831,&
         0.37831,   0.37832,   0.37832,   0.37832,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37832,   0.37832,&
         0.37831,   0.37831,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37831,   0.37831,   0.37832,&
         0.37832,   0.37832,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37830,   0.37831,   0.37831,&
         0.37832,   0.37832,   0.37832,   0.37831,   0.37831,&
         0.37830,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37831,   0.37831,   0.37831,   0.37831,   0.37831,&
         0.37830,   0.37830,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37831,   0.37831,   0.37831,   0.37832,&
         0.37831,   0.37831,   0.37830,   0.37830,   0.37830,&
         0.37830,   0.37830,   0.37830,   0.37831,   0.37832,&
         0.37830,   0.37832,   0.37831,   0.37831,   0.37831,&
         0.37831,   0.37830,   0.37830,   0.37830,   0.37831,&
         0.37832,   0.37832,   0.37832,   0.37833,   0.37832,&
         0.37831,   0.37831,   0.37830,   0.37831,   0.37831,&
         0.37831,   0.37831,   0.37832,   0.37831,   0.37831,&
         0.37832,   0.37832,   0.37831,   0.37831,   0.37831,&
         0.37832,   0.37833,   0.37832,   0.37832,   0.37832,&
         0.37832,   0.37832,   0.37832,   0.37831,   0.37832,&
         0.37833,   0.37832,   0.37832,   0.37833,   0.37833,&
         0.37832,   0.37831,   0.37832,   0.37832,   0.37832,&
         0.37832,   0.37833,   0.37834,   0.37833,   0.37832,&
         0.37833,   0.37833,   0.37832,   0.37832,   0.37832,&
         0.37832,   0.37833,   0.37833,   0.37833,   0.37833,&
         0.37833,   0.37832,   0.37833,   0.37832,   0.37832,&
         0.37831,   0.37832,   0.37833,   0.37834,   0.37835,&
         0.37835,   0.37835,   0.37833,   0.37834,   0.37832,&
         0.37832,   0.37832,   0.37832,   0.37833,   0.37835,&
         0.37834,   0.37833,   0.37834,   0.37833,   0.37832,&
         0.37832,   0.37832,   0.37832,   0.37833,   0.37834,&
         0.37835,   0.37835,   0.37835,   0.37833,   0.37833,&
         0.37832,   0.37832,   0.37832,   0.37832    /)

   real, parameter :: CHOUBAND6(firstYear:finalYear) = (/ &
         0.33338,   0.33338,   0.33334,   0.33332,   0.33332,&
         0.33335,   0.33339,   0.33340,   0.33341,   0.33340,&
         0.33340,   0.33340,   0.33340,   0.33340,   0.33340,&
         0.33338,   0.33339,   0.33340,   0.33340,   0.33340,&
         0.33342,   0.33341,   0.33342,   0.33341,   0.33342,&
         0.33341,   0.33342,   0.33342,   0.33337,   0.33336,&
         0.33341,   0.33342,   0.33339,   0.33341,   0.33341,&
         0.33342,   0.33342,   0.33342,   0.33342,   0.33342,&
         0.33342,   0.33342,   0.33342,   0.33342,   0.33342,&
         0.33342,   0.33342,   0.33342,   0.33343,   0.33343,&
         0.33342,   0.33342,   0.33343,   0.33343,   0.33343,&
         0.33343,   0.33343,   0.33343,   0.33343,   0.33343,&
         0.33343,   0.33342,   0.33342,   0.33343,   0.33343,&
         0.33343,   0.33342,   0.33343,   0.33343,   0.33343,&
         0.33342,   0.33343,   0.33343,   0.33343,   0.33342,&
         0.33343,   0.33342,   0.33343,   0.33342,   0.33343,&
         0.33343,   0.33342,   0.33342,   0.33342,   0.33342,&
         0.33342,   0.33342,   0.33342,   0.33342,   0.33342,&
         0.33342,   0.33342,   0.33342,   0.33342,   0.33342,&
         0.33342,   0.33342,   0.33342,   0.33342,   0.33342,&
         0.33342,   0.33342,   0.33342,   0.33342,   0.33342,&
         0.33342,   0.33341,   0.33341,   0.33341,   0.33339,&
         0.33340,   0.33340,   0.33341,   0.33341,   0.33340,&
         0.33340,   0.33339,   0.33339,   0.33337,   0.33340,&
         0.33336,   0.33341,   0.33340,   0.33341,   0.33341,&
         0.33340,   0.33337,   0.33339,   0.33340,   0.33337,&
         0.33340,   0.33337,   0.33340,   0.33340,   0.33341,&
         0.33341,   0.33341,   0.33341,   0.33336,   0.33336,&
         0.33337,   0.33338,   0.33339,   0.33339,   0.33340,&
         0.33340,   0.33340,   0.33339,   0.33337,   0.33337,&
         0.33337,   0.33335,   0.33337,   0.33338,   0.33338,&
         0.33339,   0.33340,   0.33338,   0.33335,   0.33333,&
         0.33333,   0.33334,   0.33335,   0.33337,   0.33338,&
         0.33339,   0.33339,   0.33337,   0.33334,   0.33334,&
         0.33336,   0.33335,   0.33337,   0.33338,   0.33340,&
         0.33339,   0.33336,   0.33334,   0.33334,   0.33334,&
         0.33336,   0.33337,   0.33337,   0.33338,   0.33338,&
         0.33340,   0.33340,   0.33340,   0.33341,   0.33341,&
         0.33340,   0.33337,   0.33338,   0.33340,   0.33340,&
         0.33340,   0.33341,   0.33341,   0.33341,   0.33342,&
         0.33342,   0.33342,   0.33341,   0.33341,   0.33341,&
         0.33340,   0.33339,   0.33339,   0.33340,   0.33340,&
         0.33341,   0.33341,   0.33341,   0.33341,   0.33340,&
         0.33340,   0.33339,   0.33337,   0.33336,   0.33336,&
         0.33336,   0.33338,   0.33339,   0.33340,   0.33340,&
         0.33337,   0.33333,   0.33332,   0.33334,   0.33335,&
         0.33336,   0.33338,   0.33339,   0.33339,   0.33339,&
         0.33338,   0.33337,   0.33336,   0.33333,   0.33334,&
         0.33336,   0.33336,   0.33336,   0.33338,   0.33339,&
         0.33340,   0.33340,   0.33339,   0.33337,   0.33334,&
         0.33334,   0.33335,   0.33336,   0.33337,   0.33338,&
         0.33339,   0.33339,   0.33340,   0.33338,   0.33336,&
         0.33333,   0.33334,   0.33334,   0.33336,   0.33338,&
         0.33339,   0.33340,   0.33340,   0.33340,   0.33340,&
         0.33338,   0.33337,   0.33337,   0.33336,   0.33336,&
         0.33337,   0.33339,   0.33340,   0.33340,   0.33340,&
         0.33340,   0.33337,   0.33336,   0.33334,   0.33334,&
         0.33335,   0.33337,   0.33338,   0.33339,   0.33339,&
         0.33340,   0.33340,   0.33340,   0.33338,   0.33336,&
         0.33336,   0.33336,   0.33336,   0.33336,   0.33337,&
         0.33338,   0.33339,   0.33340,   0.33340,   0.33339,&
         0.33336,   0.33334,   0.33332,   0.33333,   0.33335,&
         0.33337,   0.33338,   0.33339,   0.33339,   0.33338,&
         0.33336,   0.33334,   0.33333,   0.33334,   0.33334,&
         0.33336,   0.33337,   0.33338,   0.33339,   0.33338,&
         0.33336,   0.33332,   0.33331,   0.33331,   0.33332,&
         0.33334,   0.33335,   0.33336,   0.33337,   0.33337,&
         0.33335,   0.33332,   0.33329,   0.33329,   0.33329,&
         0.33333,   0.33334,   0.33336,   0.33337,   0.33337,&
         0.33336,   0.33329,   0.33325,   0.33325,   0.33327,&
         0.33330,   0.33334,   0.33336,   0.33337,   0.33337,&
         0.33337,   0.33335,   0.33332,   0.33331,   0.33331,&
         0.33330,   0.33333,   0.33333,   0.33336,   0.33336,&
         0.33338,   0.33337,   0.33336,   0.33331,   0.33327,&
         0.33327,   0.33326,   0.33329,   0.33332,   0.33335,&
         0.33337,   0.33337,   0.33336,   0.33332,   0.33327,&
         0.33328,   0.33328,   0.33331,   0.33334,   0.33336,&
         0.33337,   0.33338,   0.33336,   0.33333,   0.33330,&
         0.33328,   0.33328,   0.33327,   0.33332,   0.33334,&
         0.33335,   0.33336,   0.33337,   0.33337    /)

   real, parameter :: CHOUBAND7(firstYear:finalYear) = (/ &
         0.16512,   0.16512,   0.16508,   0.16507,   0.16506,&
         0.16509,   0.16513,   0.16515,   0.16515,   0.16514,&
         0.16514,   0.16514,   0.16514,   0.16514,   0.16515,&
         0.16512,   0.16513,   0.16514,   0.16514,   0.16515,&
         0.16516,   0.16516,   0.16516,   0.16515,   0.16516,&
         0.16516,   0.16517,   0.16517,   0.16511,   0.16511,&
         0.16516,   0.16517,   0.16513,   0.16516,   0.16516,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16517,   0.16517,   0.16517,   0.16517,   0.16517,&
         0.16516,   0.16516,   0.16515,   0.16516,   0.16514,&
         0.16515,   0.16515,   0.16515,   0.16516,   0.16515,&
         0.16515,   0.16513,   0.16513,   0.16511,   0.16514,&
         0.16511,   0.16516,   0.16514,   0.16516,   0.16516,&
         0.16514,   0.16512,   0.16514,   0.16514,   0.16512,&
         0.16515,   0.16511,   0.16514,   0.16515,   0.16516,&
         0.16516,   0.16515,   0.16515,   0.16511,   0.16510,&
         0.16511,   0.16513,   0.16513,   0.16513,   0.16514,&
         0.16515,   0.16514,   0.16513,   0.16512,   0.16511,&
         0.16511,   0.16509,   0.16511,   0.16512,   0.16512,&
         0.16514,   0.16514,   0.16512,   0.16509,   0.16507,&
         0.16507,   0.16508,   0.16509,   0.16511,   0.16512,&
         0.16513,   0.16513,   0.16511,   0.16508,   0.16508,&
         0.16510,   0.16508,   0.16511,   0.16512,   0.16514,&
         0.16513,   0.16510,   0.16508,   0.16508,   0.16509,&
         0.16510,   0.16511,   0.16512,   0.16512,   0.16513,&
         0.16514,   0.16514,   0.16515,   0.16515,   0.16515,&
         0.16515,   0.16512,   0.16513,   0.16514,   0.16514,&
         0.16514,   0.16515,   0.16516,   0.16516,   0.16516,&
         0.16517,   0.16516,   0.16516,   0.16515,   0.16515,&
         0.16515,   0.16514,   0.16514,   0.16514,   0.16515,&
         0.16515,   0.16515,   0.16515,   0.16515,   0.16515,&
         0.16514,   0.16513,   0.16512,   0.16511,   0.16510,&
         0.16510,   0.16512,   0.16513,   0.16514,   0.16514,&
         0.16511,   0.16507,   0.16506,   0.16508,   0.16509,&
         0.16510,   0.16512,   0.16513,   0.16514,   0.16513,&
         0.16512,   0.16511,   0.16510,   0.16507,   0.16508,&
         0.16510,   0.16510,   0.16511,   0.16512,   0.16513,&
         0.16514,   0.16514,   0.16513,   0.16511,   0.16508,&
         0.16507,   0.16509,   0.16510,   0.16511,   0.16512,&
         0.16513,   0.16514,   0.16514,   0.16512,   0.16510,&
         0.16507,   0.16508,   0.16508,   0.16511,   0.16512,&
         0.16514,   0.16514,   0.16514,   0.16514,   0.16514,&
         0.16513,   0.16511,   0.16511,   0.16510,   0.16510,&
         0.16511,   0.16513,   0.16514,   0.16514,   0.16515,&
         0.16514,   0.16512,   0.16510,   0.16508,   0.16508,&
         0.16509,   0.16511,   0.16513,   0.16513,   0.16514,&
         0.16514,   0.16515,   0.16515,   0.16512,   0.16510,&
         0.16511,   0.16509,   0.16510,   0.16510,   0.16511,&
         0.16513,   0.16514,   0.16514,   0.16514,   0.16513,&
         0.16510,   0.16508,   0.16506,   0.16507,   0.16509,&
         0.16511,   0.16512,   0.16513,   0.16513,   0.16512,&
         0.16510,   0.16509,   0.16507,   0.16508,   0.16508,&
         0.16509,   0.16511,   0.16512,   0.16513,   0.16512,&
         0.16510,   0.16505,   0.16505,   0.16505,   0.16506,&
         0.16507,   0.16509,   0.16510,   0.16511,   0.16511,&
         0.16509,   0.16506,   0.16503,   0.16502,   0.16503,&
         0.16506,   0.16509,   0.16510,   0.16511,   0.16511,&
         0.16509,   0.16502,   0.16498,   0.16499,   0.16501,&
         0.16503,   0.16508,   0.16510,   0.16510,   0.16511,&
         0.16510,   0.16508,   0.16505,   0.16505,   0.16504,&
         0.16504,   0.16507,   0.16506,   0.16509,   0.16510,&
         0.16512,   0.16511,   0.16509,   0.16504,   0.16500,&
         0.16500,   0.16499,   0.16503,   0.16505,   0.16509,&
         0.16511,   0.16511,   0.16510,   0.16506,   0.16500,&
         0.16501,   0.16501,   0.16504,   0.16507,   0.16509,&
         0.16510,   0.16511,   0.16510,   0.16506,   0.16503,&
         0.16500,   0.16501,   0.16500,   0.16505,   0.16507,&
         0.16509,   0.16510,   0.16510,   0.16511    /)

   real, parameter :: CHOUBAND8(firstYear:finalYear) = (/ &
         0.04348,   0.04348,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04349,   0.04349,   0.04349,&
         0.04349,   0.04349,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04347,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04347,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04347,   0.04347,&
         0.04347,   0.04347,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04347,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04347,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04347,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04347,   0.04347,&
         0.04347,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04347,&
         0.04348,   0.04347,   0.04348,   0.04347,   0.04348,&
         0.04348,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04347,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04348,   0.04348,   0.04348,   0.04348,&
         0.04347,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04347,   0.04348,   0.04348,   0.04348,&
         0.04347,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04348,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04346,   0.04346,&
         0.04347,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04346,   0.04346,   0.04346,   0.04346,&
         0.04346,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04347,   0.04346,&
         0.04346,   0.04347,   0.04347,   0.04347,   0.04347,&
         0.04348,   0.04347,   0.04347,   0.04346,   0.04346,&
         0.04346,   0.04346,   0.04346,   0.04346,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04347,   0.04346,&
         0.04346,   0.04346,   0.04346,   0.04347,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04347,   0.04346,&
         0.04346,   0.04346,   0.04346,   0.04346,   0.04347,&
         0.04347,   0.04347,   0.04347,   0.04347    /)

   ! Establish location in table for current, previous and next year.
   ! ----------------------------------------------------------------

   if(Year>finalYear .or. (Year==finalYear .and. dayOfYear > 182) ) then
      Current = (finalYear - SolCycle) + mod(Year - finalYear, SolCycle)
   else
      Current = Year
   endif

   ! Divide the year into halves.
   ! ---------------------------

   if(dayOfYear <= 182 .or. Current<firstYear) then
      i1 = max(Current-1,firstYear)
      i2 = max(Current  ,firstYear)
      F  = (dayOfYear+183)/365.00
   else
      i1 = Current
      i2 = Current+1
      F  = (dayOfYear-183)/365.00
   end if

   ! Linear interpolation to the given day-of-year.
   ! ----------------------------------------------

   SC    = so (i1)*(1.-F) + so (i2)*F

   if(present(HK)) then
      if(size(HK)==8) then ! Chou Bands
         HK(1) = ChouBand1(i1)*(1.-F) + ChouBand1(i2)*F
         HK(2) = ChouBand2(i1)*(1.-F) + ChouBand2(i2)*F
         HK(3) = ChouBand3(i1)*(1.-F) + ChouBand3(i2)*F
         HK(4) = ChouBand4(i1)*(1.-F) + ChouBand4(i2)*F
         HK(5) = ChouBand5(i1)*(1.-F) + ChouBand5(i2)*F
         HK(6) = ChouBand6(i1)*(1.-F) + ChouBand6(i2)*F
         HK(7) = ChouBand7(i1)*(1.-F) + ChouBand7(i2)*F
         HK(8) = ChouBand8(i1)*(1.-F) + ChouBand8(i2)*F
         _ASSERT(abs(1.0-sum(HK))<1.e-4,'Chou Solar band weightings do not sum to unity!')
      else
         _FAIL('HK: Solar band weightings only available for Chou')
      endif
   end if

   _RETURN(ESMF_SUCCESS)
   end subroutine MAPL_SunGetSolarConstantByYearDoY

!-------------------------------------------------------------------------
!NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1, GEOS/DAS!
!-------------------------------------------------------------------------
!>
! Acquire the solar constant and the eight band fractions for the Chou solar code
! from an external NetCDF file.  The initial configuration assumes monthly values
! beginning in January, and does not interpolate beyond the range of available data.
!
!@bug
! Band values for RRTMG not implemented
!@endbug
!
   subroutine MAPL_SunGetSolarConstantFromNetcdfFile(CLOCK,fileName,SC,HK,MESOPHOT,JCALC4,rc)

      implicit none

      type(ESMF_Clock), intent(in)           :: CLOCK
      character(len=*), intent(in)           :: fileName
      real, intent(out)                      :: SC
      real, optional, intent(out)            :: HK(:)
      real, optional, intent(out)            :: MESOPHOT(:)
      real, optional, intent(out)            :: JCALC4(:)
      integer, optional, intent(out)         :: rc

      type(ESMF_Time)            :: time
      integer                    :: N
      integer                    :: begYear, endYear
      integer                    :: INDX1, INDX2
      integer                    :: MM, YY, DD, CCYY

      integer                    :: STATUS

      real                       :: FAC

      integer :: ncid

      integer :: ndate,          dimid_ndate
      integer :: nbin_sorad,     dimid_nbin_sorad
      integer :: nbin_meso_phot, dimid_nbin_meso_phot
      integer :: nbin_jcalc4,    dimid_nbin_jcalc4

      integer, dimension(:  ), allocatable :: date_year
      integer                              :: varid_date_year

      integer, dimension(:  ), allocatable :: date_month
      integer                              :: varid_date_month

      real,    dimension(:  ), allocatable :: tsi
      integer                              :: varid_tsi

      real,    dimension(:,:), allocatable :: coef_sorad
      integer                              :: varid_coef_sorad

      real,    dimension(:,:), allocatable :: coef_meso_phot
      integer                              :: varid_coef_meso_phot

      real,    dimension(:,:), allocatable :: coef_jcalc4
      integer                              :: varid_coef_jcalc4

      character(len=ESMF_MAXSTR) :: IAm = "MAPL_SunGetSolarConstantFromNetcdfFile"

      ! Open the file
      ! -------------
      STATUS = nf90_open(trim(fileName), NF90_NOWRITE, ncid)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error opening file ', trim(fileName), status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      ! Read in dimensions
      ! ------------------

      status = nf90_inq_dimid(ncid, 'ndate', dimid_ndate)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting ndate dimid', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      status = nf90_inquire_dimension(ncid, dimid_ndate, len = ndate)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting ndate length', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      if (present(HK)) then
         status = nf90_inq_dimid(ncid, 'nbin_sorad', dimid_nbin_sorad)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting nbin_sorad dimid', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

         status = nf90_inquire_dimension(ncid, dimid_nbin_sorad, len = nbin_sorad)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting nbin_sorad length', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if
      end if

      if (present(MESOPHOT)) then
         status = nf90_inq_dimid(ncid, 'nbin_meso_phot', dimid_nbin_meso_phot)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting nbin_meso_phot dimid', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

         status = nf90_inquire_dimension(ncid, dimid_nbin_meso_phot, len = nbin_meso_phot)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting nbin_meso_phot length', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if
      end if

      if (present(JCALC4)) then
         status = nf90_inq_dimid(ncid, 'nbin_jcalc4', dimid_nbin_jcalc4)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting nbin_jcalc4 dimid', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

         status = nf90_inquire_dimension(ncid, dimid_nbin_jcalc4, len = nbin_jcalc4)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting nbin_jcalc4 length', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if
      end if

      ! Allocate our arrays
      ! -------------------

      allocate(date_year(ndate), source=0, stat=status)
      _VERIFY(STATUS)

      allocate(date_month(ndate), source=0, stat=status)
      _VERIFY(STATUS)

      allocate(tsi(ndate), source=0.0, stat=status)
      _VERIFY(STATUS)

      if (present(HK) ) then
         allocate(coef_sorad(nbin_sorad,ndate), source=0.0, stat=status)
         _VERIFY(STATUS)
      end if

      if (present(MESOPHOT) ) then
         allocate(coef_meso_phot(nbin_meso_phot,ndate), source=0.0, stat=status)
         _VERIFY(STATUS)
      end if

      if (present(JCALC4) ) then
         allocate(coef_jcalc4(nbin_jcalc4,ndate), source=0.0, stat=status)
         _VERIFY(STATUS)
      end if

      ! Read in date_year
      ! -----------------

      status = nf90_inq_varid(ncid, 'date_year', varid_date_year)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting date_year varid', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      status = nf90_get_var(ncid, varid_date_year, date_year)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting date_year variable', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      ! Read in date_month
      ! ------------------

      status = nf90_inq_varid(ncid, 'date_month', varid_date_month)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting date_month varid', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      status = nf90_get_var(ncid, varid_date_month, date_month)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting date_month variable', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      ! Read in tsi
      ! -----------

      status = nf90_inq_varid(ncid, 'tsi', varid_tsi)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting tsi varid', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      status = nf90_get_var(ncid, varid_tsi, tsi)
      if (STATUS /= NF90_NOERR) then
         write (*,*) trim(Iam)//': Error getting tsi variable', status
         write (*,*) nf90_strerror(status)
         _FAIL('needs informative message')
      end if

      ! Read in coef_sorad
      ! ------------------

      if (present(HK)) then

         status = nf90_inq_varid(ncid, 'coef_sorad', varid_coef_sorad)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting coef_sorad varid', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

         status = nf90_get_var(ncid, varid_coef_sorad, coef_sorad)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting coef_sorad variable', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

      end if

      ! Read in coef_meso_phot
      ! ----------------------

      if (present(MESOPHOT)) then

         status = nf90_inq_varid(ncid, 'coef_meso_phot', varid_coef_meso_phot)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting coef_meso_phot varid', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

         status = nf90_get_var(ncid, varid_coef_meso_phot, coef_meso_phot)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting coef_meso_phot variable', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

      end if

      ! Read in coef_jcalc4
      ! -------------------

      if (present(JCALC4)) then

         status = nf90_inq_varid(ncid, 'coef_jcalc4', varid_coef_jcalc4)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting coef_jcalc4 varid', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

         status = nf90_get_var(ncid, varid_coef_jcalc4, coef_jcalc4)
         if (STATUS /= NF90_NOERR) then
            write (*,*) trim(Iam)//': Error getting coef_jcalc4 variable', status
            write (*,*) nf90_strerror(status)
            _FAIL('needs informative message')
         end if

      end if

      ! Time interpolation parameters
      ! -----------------------------
      call MAPL_ClimInterpFac(CLOCK, INDX1, INDX2, FAC, RC=STATUS)
      _VERIFY(STATUS)

      ! Time: Year and month
      ! --------------------
      call ESMF_ClockGet(CLOCK, CURRTIME=time, RC=STATUS)
      _VERIFY(STATUS)
      call ESMF_TimeGet(time, YY=YY, MM=MM, DD=DD, RC=STATUS)
      _VERIFY(STATUS)

      ! This index search routine below is adaptedfrom PChemGridComp
      !
      ! Read bracketing months, make sure INDX1 and INDX2 are in range. Annual
      ! cycle is preserved for years that precede and succeed the climatology.
      ! ----------------------------------------------------------------------

       N = 12*ndate
       begYear = date_year(1)
       endYear = date_year(ndate)

       CCYY = YY
       IF(CCYY < begYear) CCYY = begYear
       IF(CCYY > endYear) CCYY = endYear

       INDX1 = INDX1+(CCYY-begYear)*12
       INDX2 = INDX2+(CCYY-begYear)*12

       IF(MM ==  1 .AND. INDX1 > INDX2) INDX1 = INDX1-12
       IF(MM == 12 .AND. INDX2 < INDX1) INDX2 = INDX2+12

       IF(INDX1 == 0) INDX1 = 12
       IF(YY < begYear .AND. INDX2 > 12) INDX2 = 1

       IF(INDX2 == N+1) INDX2 = N-11
       IF(YY > endYear .AND. INDX1 == N-12) INDX1 = N


      ! Linear Interpolation to the given day-of-month
      ! ----------------------------------------------

      SC = tsi(INDX1)*FAC + tsi(INDX2)*(1.0-FAC)

      if (present(HK)) then
         HK = coef_sorad(:,INDX1)*FAC + coef_sorad(:,INDX2)*(1.0-FAC)
      end if

      if (present(MESOPHOT)) then
         HK = coef_meso_phot(:,INDX1)*FAC + coef_meso_phot(:,INDX2)*(1.0-FAC)
      end if

      if (present(JCALC4)) then
         HK = coef_jcalc4(:,INDX1)*FAC + coef_jcalc4(:,INDX2)*(1.0-FAC)
      end if

      ! Close the file
      ! --------------
      STATUS = nf90_close(ncid)
      _VERIFY(STATUS)

      ! Bounds check
      ! ------------
      if (PRESENT(HK)) then
         _ASSERT(ABS(1.0-SUM(HK)) < 1.E-4,'needs informative message')
      end if

      ! Deallocate our arrays
      ! ---------------------

      deallocate(date_year, stat=status)
      _VERIFY(STATUS)

      deallocate(date_month, stat=status)
      _VERIFY(STATUS)

      deallocate(tsi, stat=status)
      _VERIFY(STATUS)

      if (present(HK)) then
         deallocate(coef_sorad, stat=status)
         _VERIFY(STATUS)
      end if

      if (present(MESOPHOT)) then
         deallocate(coef_meso_phot, stat=status)
         _VERIFY(STATUS)
      end if

      if (present(JCALC4)) then
         deallocate(coef_jcalc4, stat=status)
         _VERIFY(STATUS)
      end if

      _RETURN(ESMF_SUCCESS)

   end subroutine MAPL_SunGetSolarConstantFromNetcdfFile

!-------------------------------------------------------------------------
!NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1, GEOS/DAS!
!-------------------------------------------------------------------------
!>
! Acquire the TSI, Mg, and SB from file.
!
! @bug
! Band values for RRTMG not implemented
!@endbug
!
   subroutine MAPL_SunGetSolarConstantFromNRLFile(CLOCK,filename_in,SC,MG,SB,PersistSolar,rc)

      implicit none

      type(ESMF_Clock), intent(in)       :: CLOCK
      character(len=*), intent(in)       :: filename_in
      real, intent(out)                  :: SC
      real, intent(out)                  :: MG
      real, intent(out)                  :: SB
      logical, optional, intent(in)      :: PersistSolar
      integer, optional, intent(out)     :: rc

      type(ESMF_VM)              :: VM
      type(ESMF_Time)            :: noonCurrentDay, prevNoon, nextNoon
      type(ESMF_TimeInterval)    :: intToNextNoon, oneDayInterval

      type(ESMF_Time)            :: currentTime, origTime
      type(ESMF_Time)            :: timeBasedOnCycle24
      type(ESMF_Time)            :: startCycle24, startCycle25
      type(ESMF_TimeInterval)    :: timeSinceStartOfCycle25

      integer :: currentYear, currentMon, currentDay, currentDOY
      integer :: prevDOY, nextDOY, prevNoonYear, nextNoonYear
      integer :: originalYear, originalMon, originalDay, origDOY

      real(ESMF_KIND_R8)         :: days_r8

      integer                    :: INDX1, INDX2, i

      integer                    :: STATUS, UNIT

      real                       :: FAC

      character(len=ESMF_MAXPATHLEN) :: FILENAME
      logical :: amIRoot
      integer :: deId, NPES
      logical :: outOfTable

      logical :: PersistSolar_

      integer, parameter :: YEAR_NOT_FOUND = -99

      character(len=256) :: line
      character(len=1), parameter :: commentChar = "#"
      integer :: stat

      logical, save :: TableCreated = .FALSE.
      integer, save, allocatable, dimension(:) :: yearTable, doyTable
      real,    save, allocatable, dimension(:) :: tsi, mgindex, sbindex
      integer, save :: numlines

      class(Logger), pointer :: lgr

      lgr => logging%get_logger('MAPL.SUN')

      if (present(PersistSolar)) then
         PersistSolar_ = PersistSolar
      else
         PersistSolar_ = .TRUE.
      end if

      call ESMF_VMGetCurrent(vm, rc=status)

      call ESMF_VmGet(VM, localPet=deId, petCount=npes, rc=status)
      _VERIFY(STATUS)
      amIRoot = (deId == 0)

      CREATE_TABLE: if (.not. TableCreated) then

         ! First we open the file on root to get the
         ! number of lines so we can allocate our tables
         ! ---------------------------------------------

         if (amIRoot) then

            ! Open the file
            ! -------------
            filename = trim(filename_in)
            open(newunit=unit, file=filename, form="formatted", status="old", iostat=status)
            _ASSERT(status==0,'Could not find NRL data file '// trim(filename ))

            ! Determine length of file
            ! ------------------------
            call lgr%debug("Scanning the Solar Table to determine number of data points")
            numlines = num_lines_in_file(UNIT)
            call lgr%debug("Solar Table Data Points: %i0", numlines)

         end if

         ! Broadcast the number of lines
         ! -----------------------------
         call MAPL_CommsBcast(vm, DATA=numlines, N=1, ROOT=0, _RC)

         ! Allocate our arrays on all processes
         ! ------------------------------------

         allocate(yearTable(numlines), source=0,   _STAT)
         allocate(doyTable(numlines),  source=0,   _STAT)
         allocate(tsi(numlines),       source=0.0, _STAT)
         allocate(mgindex(numlines),   source=0.0, _STAT)
         allocate(sbindex(numlines),   source=0.0, _STAT)

         ! Back to root to read in the values
         ! ----------------------------------

         if (amIRoot) then

            call lgr%debug("Reading the Solar Table")
            i = 1
            do
               read(unit,'(A)',iostat=stat) line
               if (is_iostat_end(stat)) exit
               ! Skip lines starting with #-comments
               if (index (line, commentChar) /= 0) cycle
               read(line,*) yearTable(i), doyTable(i), tsi(i), mgindex(i), sbindex(i)
               i = i + 1
            end do

            ! Belt and suspenders check that all data was read
            _ASSERT(size(yearTable) == numlines,"Inconsistency in NRL number of lines")

            close(unit, _IOSTAT)

         end if

         ! Broadcast the tables
         ! --------------------

         call MAPL_CommsBcast(vm, DATA=yearTable, N=numlines, ROOT=0, _RC)
         call MAPL_CommsBcast(vm, DATA=doyTable,  N=numlines, ROOT=0, _RC)
         call MAPL_CommsBcast(vm, DATA=tsi,       N=numlines, ROOT=0, _RC)
         call MAPL_CommsBcast(vm, DATA=mgindex,   N=numlines, ROOT=0, _RC)
         call MAPL_CommsBcast(vm, DATA=sbindex,   N=numlines, ROOT=0, _RC)

         TableCreated = .TRUE.

      end if CREATE_TABLE

      ! Now we need to find the two bracketing days
      ! -------------------------------------------

      ! Get current time
      ! ----------------
      call ESMF_ClockGet(CLOCK, CURRTIME=currentTime, _RC)

      call ESMF_TimeGet( currentTime, YY = currentYear,   &
                                      MM = currentMon,    &
                                      DD = currentDay,    &
                               dayOfYear = currentDOY, _RC)

      ! Test if current time is outside our file
      ! ----------------------------------------

      outOfTable = .FALSE.

      ! First is current year higher than last in file...
      if ( currentYear > yearTable(numlines) ) then
         outOfTable = .TRUE.
      ! ...or if a partial year, are we near the end
      else if ( currentYear == yearTable(numlines) .and. currentDOY >= doyTable(numlines)) then
         outOfTable = .TRUE.
      end if

      ! If we are out of the table...
      ! -----------------------------

      OUT_OF_TABLE: if ( outOfTable ) then

         PERSIST_SOLAR: if ( PersistSolar_ ) then

            ! If we are outOfTable and we have the PersistSolar
            ! option we just use the last value in the table...
            ! -------------------------------------------------

            SC =     tsi(numlines)
            MG = mgindex(numlines)
            SB = sbindex(numlines)

            call lgr%debug('Off the end of table, persisting last values')
            call lgr%debug('  tsi at end of table: %F8.3', tsi(numlines))
            call lgr%debug('  mgindex at end of table: %F8.6', mgindex(numlines))
            call lgr%debug('  sbindex at end of table: %F9.4', sbindex(numlines))

            _RETURN(ESMF_SUCCESS)

         else

            ! If we are out of the table and not persisting, we must
            ! recenter our day to be based on the last complete Solar Cycle
            ! -------------------------------------------------------------

            ! Create an ESMF_Time at start of Cycle 24
            ! ----------------------------------------
            call ESMF_TimeSet( startCycle24, YY = 2008,  &
                                             MM = 12,    &
                                             DD = 1,     &
                                              H = 12,    &
                                              M = 00,    &
                                              S = 00, _RC)

            ! Create an ESMF_Time at start of Cycle 25
            ! ----------------------------------------
            call ESMF_TimeSet( startCycle25, YY = 2019,  &
                                             MM = 12,    &
                                             DD = 1,     &
                                              H = 12,    &
                                              M = 00,    &
                                              S = 00, _RC)

            ! Create TimeInterval based on interval
            ! from start of latest Cycle 25
            ! -------------------------------------

            timeSinceStartOfCycle25 = currentTime - startCycle25

            ! Make a new time based on that
            ! interval past start of Cycle 24
            ! -------------------------------

            timeBasedOnCycle24 = startCycle24 + timeSinceStartOfCycle25

            ! Store our original time just in case
            ! ------------------------------------
            origTime     = currentTime
            originalYear = currentYear
            originalMon  = currentMon
            originalDay  = currentDay
            origDOY      = currentDOY

            ! Make our "current" time the one calculated above
            ! ------------------------------------------------
            currentTime = timeBasedOnCycle24

            ! Get new currentYear, currentMon, currentDay
            ! -------------------------------------------

            call ESMF_TimeGet( currentTime, YY = currentYear,   &
                                            MM = currentMon,    &
                                            DD = currentDay,    &
                                     dayOfYear = currentDOY, _RC)

            call lgr%debug('Off the end of table, moving into last complete cycle')
            call lgr%debug('  Original Year-Mon-Day to Find: %i0.4~-%i0.2~-%i0.2', originalYear,originalMon,originalDay)
            call lgr%debug('  Original Day of Year: %i0', origDOY)
            call lgr%debug('  New Year-Mon-Day to Find: %i0.4~-%i0.2~-%i0.2', currentYear,currentMon,currentDay)
            call lgr%debug('  New Day of Year: %i0', currentDOY)

         end if PERSIST_SOLAR

      end if OUT_OF_TABLE

      ! Create an ESMF_Time on noon of current day
      ! ------------------------------------------
      call ESMF_TimeSet( noonCurrentDay, YY = currentYear, &
                                         MM = currentMon,  &
                                         DD = currentDay,  &
                                          H = 12,          &
                                          M = 00,          &
                                          S = 00,       _RC)

      ! Figure out bracketing days for interpolation
      ! NOTE: nextNoon is mainly for debugging purposes
      ! -----------------------------------------------
      call ESMF_TimeIntervalSet(oneDayInterval, D=1, _RC)
      if (currentTime <= noonCurrentDay) then
         prevNoon = noonCurrentDay - oneDayInterval
         nextNoon = noonCurrentDay
      else
         prevNoon = noonCurrentDay
         nextNoon = noonCurrentDay + oneDayInterval
      end if

      ! Get the DOYs
      ! ------------
      call ESMF_TimeGet( prevNoon, YY = prevNoonYear, dayOfYear = prevDOY, _RC)
      call ESMF_TimeGet( nextNoon, YY = nextNoonYear, dayOfYear = nextDOY, _RC)

      ! Our interpolation factor is based of when we are compared to the next noon
      ! --------------------------------------------------------------------------
      intToNextNoon = nextNoon-currentTime

      ! The FAC for interpolating is just the real version
      ! of the size of the timeinterval to the next noon
      ! --------------------------------------------------
      call ESMF_TimeIntervalGet(intToNextNoon, d_r8=days_r8, _RC)
      FAC = real(days_r8)

      ! Use our find_file_index function to get the index for previous noon
      ! -------------------------------------------------------------------
      INDX1 = find_file_index(numlines, yearTable, prevNoonYear, prevDOY)
      INDX2 = INDX1 + 1

      ! Linear Interpolation to the given day-of-month
      ! ----------------------------------------------

      SC =     tsi(INDX1)*FAC +     tsi(INDX2)*(1.0-FAC)
      MG = mgindex(INDX1)*FAC + mgindex(INDX2)*(1.0-FAC)
      SB = sbindex(INDX1)*FAC + sbindex(INDX2)*(1.0-FAC)

      call lgr%debug('  First DOY to Find:     %i3', prevDOY)
      call lgr%debug('    file_index for date: %i6', INDX1)
      call lgr%debug('    yearTable(date):     %i4', yearTable(INDX1))
      call lgr%debug('    doyTable(date):      %i3', doyTable(INDX1))
      call lgr%debug('    tsi(date):           %f8.3', tsi(INDX1))
      call lgr%debug('    mgindex(date):       %f8.6', mgindex(INDX1))
      call lgr%debug('    sbindex(date):       %f9.4', sbindex(INDX1))

      call lgr%debug('  Second DOY to Find:    %i3', nextDOY)
      call lgr%debug('    file_index for date: %i6', INDX2)
      call lgr%debug('    yearTable(date):     %i4', yearTable(INDX2))
      call lgr%debug('    doyTable(date):      %i3', doyTable(INDX2))
      call lgr%debug('    tsi(date):           %f8.3', tsi(INDX2))
      call lgr%debug('    mgindex(date):       %f8.6', mgindex(INDX2))
      call lgr%debug('    sbindex(date):       %f9.4', sbindex(INDX2))

      call lgr%debug('  Interpolation Factor:  %f8.6', FAC)

      _RETURN(ESMF_SUCCESS)

      contains

         integer function num_lines_in_file(unit) result(count)

            implicit none

            integer, intent(in) :: unit

            character(len=256) :: line
            integer :: stat
            character(len=1), parameter :: commentChar = "#"

            count = 0

            do
               read(unit,'(A)',iostat=stat) line
               if (is_iostat_end(stat)) exit
               ! Skip lines starting with #-comments
               if (index (line, commentChar) /= 0) cycle
               count = count + 1
            end do

            rewind(unit)

         end function num_lines_in_file

         integer function find_file_index(numlines, year, year_to_find, doy_to_find)

            implicit none

            integer, intent(in) :: numlines
            integer, intent(in), dimension(numlines) :: year
            integer, intent(in) :: year_to_find, doy_to_find

            integer :: i, yr_index

            ! Find the index for the year
            ! ---------------------------

            yr_index = YEAR_NOT_FOUND

            do i = 1, numlines
               if (year(i) == year_to_find) then
                  yr_index = i
                  exit
               end if
            end do

            if (yr_index == YEAR_NOT_FOUND) then
               find_file_index = YEAR_NOT_FOUND
            else
               find_file_index = yr_index + doy_to_find - 1
            end if

         end function find_file_index


   end subroutine MAPL_SunGetSolarConstantFromNRLFile

!==========================================================================
!>
! Return the daylight duration in seconds (i.e, the time between sunrise and sunset) for
! a day around the specified time. The routine is accurate enough for most purposes, but
! does not solve for precise sunrise/sunset times influenced by changes in the orbital
! parameters between those times. The time input does NOT need to be noon --- it is used
! simply to evaluate the solar declination needed for the daylight duration calculation.
! In order of preference, time is taken from currTime, if present, or else the currTime
! of CLOCK, if present, or else the currTime of the ORBIT's associated clock.
!
   subroutine MAPL_SunGetDaylightDuration(ORBIT,LATS,DAYL,currTime,CLOCK,RC)

! !ARGUMENTS:

      type(MAPL_SunOrbit),           intent(IN ) :: ORBIT
      real, dimension(:) ,           intent(IN ) :: LATS
      real, dimension(:) ,           intent(OUT) :: DAYL
      type(ESMF_Time)    , optional, intent(IN ) :: currTime
      type(ESMF_Clock)   , optional, intent(IN ) :: CLOCK
      integer,             optional, intent(OUT) :: RC

!     Locals

      character(len=ESMF_MAXSTR), parameter :: IAm = "MAPL_SunGetDaylightDuration"
      integer                               :: STATUS

      type(ESMF_Time) :: CURRENTTIME
      integer         :: YEAR, SEC_OF_DAY, DAY_OF_YEAR, IDAY, IDAYP1
      real            :: FAC, ZS, ZC

      real(ESMF_KIND_R8) :: days
      real :: ECC, OBQ, LAMBDAP
      real :: OMECC, OPECC, OMSQECC, EAFAC
      real :: MA, EA, dE, TA, LAMBDA
      integer :: nits

!     Begin

      _ASSERT(MAPL_SunOrbitCreated(ORBIT),'MAPL_SunOrbit not yet created!')

      ! which current time?
      if (present(currTime)) then
         CURRENTTIME = CURRTIME
      else
         if (present(CLOCK)) then
            call ESMF_ClockGet(      CLOCK, currTime=CURRENTTIME, RC=STATUS)
         else
            call ESMF_ClockGet(ORBIT%CLOCK, currTime=CURRENTTIME, RC=STATUS)
         end if
         _VERIFY(STATUS)
      end if

      if (ORBIT%ANAL2B) then

         ! include time variation in orbit from reference time
         call ESMF_TimeIntervalGet( &
           CURRENTTIME - ORBIT%ORB2B_TIME_REF, &
           d_r8=days, rc=STATUS)
         _VERIFY(STATUS)
         ECC     = ORBIT%ORB2B_ECC_REF     + days * ORBIT%ORB2B_ECC_RATE
         OBQ     = ORBIT%ORB2B_OBQ_REF     + days * ORBIT%ORB2B_OBQ_RATE
         LAMBDAP = ORBIT%ORB2B_LAMBDAP_REF + days * ORBIT%ORB2B_LAMBDAP_RATE
         ! derived quantities
         OMECC = 1. - ECC
         OPECC = 1. + ECC
         OMSQECC = OMECC * OPECC
         EAFAC = sqrt(OMECC/OPECC)
         ! time interval since perhelion in days
         call ESMF_TimeIntervalGet( &
           CURRENTTIME - ORBIT%ORB2B_TIME_PERI, &
           d_r8=days, rc=STATUS)
         _VERIFY(STATUS)
         ! mean anomaly
         MA = ORBIT%ORB2B_OMG0 * days
         ! eccentric anomaly
         call invert_Keplers_Newton(MA,ECC,EA,dE,nits)
         ! true anomaly
         TA = calcTAfromEA(EA,EAFAC)
         ! celestial longitude
         LAMBDA = TA + LAMBDAP
         ! sin and cos of solar declination
         ZS = sin(LAMBDA) * sin(OBQ)
         ZC = sqrt(1. - ZS**2)

      else

         call ESMF_TimeGet(CURRENTTIME, YY=YEAR, S=SEC_OF_DAY, &
            dayOfYear=DAY_OF_YEAR, RC=STATUS)
         _VERIFY(STATUS)

         YEAR = mod(YEAR-1,ORBIT%YEARS_PER_CYCLE)
         IDAY = YEAR*int(ORBIT%YEARLEN)+DAY_OF_YEAR
         IDAYP1 = mod(IDAY,ORBIT%DAYS_PER_CYCLE) + 1

         FAC = real(SEC_OF_DAY)/86400.
         ZS = ORBIT%ZS(IDAYP1)*FAC + ORBIT%ZS(IDAY)*(1.-FAC)
         ZC = ORBIT%ZC(IDAYP1)*FAC + ORBIT%ZC(IDAY)*(1.-FAC)

      endif

      ! dayligt duration [secs]
      DAYL = (86400./MAPL_PI)*acos(min(1.,max(-1.,-tan(LATS)*ZS/ZC)))

      _RETURN(ESMF_SUCCESS)

   end subroutine MAPL_SunGetDaylightDuration

!==========================================================================
!>
! Return the daylight duration in seconds (i.e, the time between sunrise and sunset) for
! its MAXIMUM at the summer solstice. The routine is accurate enough for most purposes,
! but does not solve for precise sunrise/sunset times influenced by changes in the orbital
! parameters between those times. The time input does NOT need to be noon --- it is used
! simply to evaluate the obliquity needed for the maximum daylight duration calculation.
! In order of preference, time is taken from currTime, if present, or else the currTime
! of CLOCK, if present, or else the currTime of the ORBIT's associated clock.
! Note: Unless ORBIT_ANAL2B, the obliquity is fixed and the time is irrelevant.
!
   subroutine MAPL_SunGetDaylightDurationMax(ORBIT,LATS,DAYL,currTime,CLOCK,RC)

! !ARGUMENTS:

      type(MAPL_SunOrbit),           intent(IN ) :: ORBIT
      real, dimension(:) ,           intent(IN ) :: LATS
      real, dimension(:) ,           intent(OUT) :: DAYL
      type(ESMF_Time)    , optional, intent(IN ) :: currTime
      type(ESMF_Clock)   , optional, intent(IN ) :: CLOCK
      integer,             optional, intent(OUT) :: RC

!EOPI

      character(len=ESMF_MAXSTR), parameter :: IAm = "MAPL_SunGetDaylightDurationMax"
      integer                               :: STATUS

      type(ESMF_Time) :: CURRENTTIME
      real(ESMF_KIND_R8) :: days
      real :: OBQ

      _ASSERT(MAPL_SunOrbitCreated(ORBIT),'MAPL_SunOrbit not yet created!')

      ! Which time?
      if (present(currTime)) then
         CURRENTTIME = CURRTIME
      else
         if (present(CLOCK)) then
            call ESMF_ClockGet(      CLOCK, currTime=CURRENTTIME, RC=STATUS)
         else
            call ESMF_ClockGet(ORBIT%CLOCK, currTime=CURRENTTIME, RC=STATUS)
         end if
         _VERIFY(STATUS)
      end if

      if (ORBIT%ANAL2B) then
         ! time variation in obliquity from ref time
         call ESMF_TimeIntervalGet( &
           CURRENTTIME - ORBIT%ORB2B_TIME_REF, &
           d_r8=days, rc=STATUS)
         _VERIFY(STATUS)
         OBQ = ORBIT%ORB2B_OBQ_REF + days * ORBIT%ORB2B_OBQ_RATE
      else
         ! obliquity fixed in this case
         OBQ = ORBIT%OB * (MAPL_PI/180.)
      endif
      _ASSERT(OBQ >= 0, 'Obliquity less than 0 detected!')
      _ASSERT(OBQ < MAPL_PI, 'Obliquity greater than pi detected!')

      ! Maximum daylight duration at summer solstice [secs]
      !   (an even function of latitude)
      DAYL = (86400./MAPL_PI)*acos(min(1.,max(-1., &
                -tan(ABS(LATS))*tan(OBQ))))

      _RETURN(ESMF_SUCCESS)

   end subroutine MAPL_SunGetDaylightDurationMax

!==========================================================================
!>
! Returns the local solar hour angle (in radians) at the single time and
! multiple longitudes specified. In order of preference, time is taken
! from TIME, if present, or else the CURRTIME of CLOCK, if present, or
! else the CURRTIME of the ORBIT's associated clock.
!
! NB: For accurate results, namely to receive the TRUE local solar hour
! angle, ensure the ORBIT has the EOT flag set true. Conversely, to get
! only the MEAN local solar hour angle, use the optional argument
! FORCE_MLSHA=.TRUE.. This will turn off the Equation of Time correction
! (for this LSHA calculation only) even if the ORBIT includes it. For
! example, in the local noon detection in the EXAMPLE below, this will
! give mean local noons that are exactly 24 hours apart at a particular
! location. But they will no longer exactly be the solar culmination
! times (the TRUE local noon) in that case. TRUE local noons are not
! exactly 24h apart because of orbital variations in length of day
! throughout the year, as described by the Equation of Time.
!
! @note
! `Example of use:`:
!
!     type (ESMF_Time) :: NOW
!     type (ESMF_TimeInterval) :: DELT
!     real,    dimension(size(LONS)) :: LSHA0, LSHA1
!     logical, dimension(size(LONS)) :: isNoon
!     call ESMF_ClockGet (CLOCK, CURRTIME=NOW, TIMESTEP=DELT, __RC__)
!     call MAPL_SunGetLocalSolarHourAngle (ORBIT, LONS, LSHA0, TIME=NOW,      __RC__)
!     call MAPL_SunGetLocalSolarHourAngle (ORBIT, LONS, LSHA1, TIME=NOW+DELT, __RC__)
!     isnoon = (LSHA0 <= 0. .and. LSHA1 > 0.)
!@endnote
!
   subroutine MAPL_SunGetLocalSolarHourAngle (ORBIT,LONS,LSHA, &
      TIME,CLOCK,FORCE_MLSHA,RC)

! !ARGUMENTS:

      type (MAPL_SunOrbit),            intent(IN ) :: ORBIT
      real,    dimension(:),           intent(IN ) :: LONS   ! [radians]
      real,    dimension(:),           intent(OUT) :: LSHA
      type (ESMF_Time),      optional, intent(IN ) :: TIME
      type (ESMF_Clock),     optional, intent(IN ) :: CLOCK
      logical,               optional, intent(IN ) :: FORCE_MLSHA
      integer,               optional, intent(OUT) :: RC

!EOPI

!     Locals

      character(len=ESMF_MAXSTR), parameter :: IAm = "MAPL_SunGetLocalSolarHourAngle"
      integer                               :: STATUS

      type (ESMF_Time) :: T
      real (ESMF_KIND_R8) :: days
      integer :: YEAR, SEC_OF_DAY, DAY_OF_YEAR, IDAY, IDAYP1
      real :: DFRAC, GSHA
      real :: ECC, OBQ, LAMBDAP
      real :: OMECC, OPECC, OMSQECC, EAFAC
      real :: MA, EA, dE, TA, LAMBDA
      real :: RT, RM, ET
      integer :: i, nits
      logical :: do_EOT

      _ASSERT (MAPL_SunOrbitCreated(ORBIT),'MAPL_SunOrbit not yet created!')

      ! Which time?
      if (present(TIME)) then
         T = TIME
      else
         if (present(CLOCK)) then
            call ESMF_ClockGet (      CLOCK, CURRTIME=T, _RC)
         else
            call ESMF_ClockGet (ORBIT%CLOCK, CURRTIME=T, _RC)
         end if
      end if

      ! NB: include YY and dayOfYear here so that S is seconds WITHIN a day.
      ! YEAR and DAY_OF_YEAR are used within the non-ANAL2B branch anyway.
      call ESMF_TimeGet (T, YY=YEAR, dayOfYear=DAY_OF_YEAR, S=SEC_OF_DAY, RC=STATUS)
      _VERIFY(STATUS)

      ! fraction of day (0 at midnight)
      DFRAC = real(SEC_OF_DAY) / 86400.

      ! Greenwich MEAN solar hour angle (zero at noon)
      GSHA = 2. * MAPL_PI * (DFRAC - 0.5)

      ! Apply equation of time correction?
      do_EOT = ORBIT%EOT
      if (present(FORCE_MLSHA)) then
         if (FORCE_MLSHA) do_EOT = .FALSE.
      endif
      if (do_EOT) then

         if (ORBIT%ANAL2B) then

            ! include time variation in orbit from reference time
            call ESMF_TimeIntervalGet (T - ORBIT%ORB2B_TIME_REF, d_r8=days, _RC)
            ECC     = ORBIT%ORB2B_ECC_REF     + days * ORBIT%ORB2B_ECC_RATE
            OBQ     = ORBIT%ORB2B_OBQ_REF     + days * ORBIT%ORB2B_OBQ_RATE
            LAMBDAP = ORBIT%ORB2B_LAMBDAP_REF + days * ORBIT%ORB2B_LAMBDAP_RATE
            ! derived quantities
            OMECC = 1. - ECC
            OPECC = 1. + ECC
            OMSQECC = OMECC * OPECC
            EAFAC = sqrt(OMECC/OPECC)
            ! time interval since perhelion in days
            call ESMF_TimeIntervalGet (T - ORBIT%ORB2B_TIME_PERI, d_r8=days, _RC)
            ! mean anomaly
            MA = ORBIT%ORB2B_OMG0 * days
            ! eccentric anomaly
            call invert_Keplers_Newton (MA,ECC,EA,dE,nits)
            ! true anomaly
            TA = calcTAfromEA (EA,EAFAC)
            ! celestial longitude
            LAMBDA = TA + LAMBDAP
            ! solar right ascension (true and mean)
            RT = atan2(sin(LAMBDA)*cos(OBQ),cos(LAMBDA))
            RM = MA + LAMBDAP
            ! equation of time
            ET = RECT_PMPI (RM - RT)

         else

            ! get equation of time by table interpolation
            YEAR = mod(YEAR-1,ORBIT%YEARS_PER_CYCLE)
            IDAY = YEAR*int(ORBIT%YEARLEN)+DAY_OF_YEAR
            IDAYP1 = mod(IDAY,ORBIT%DAYS_PER_CYCLE) + 1
            ET = ORBIT%ET(IDAYP1)*DFRAC + ORBIT%ET(IDAY)*(1.-DFRAC)

         endif

         ! Gives Greenwich TRUE solar hour angle
         GSHA = GSHA + ET

      end if  ! EOT correction

      ! LOCAL solar hour angle
      do i = 1, size(LONS)
         LSHA(i) = RECT_PMPI (GSHA + LONS(i))
      end do

      _RETURN(ESMF_SUCCESS)

   end subroutine MAPL_SunGetLocalSolarHourAngle

!==========================================================================

end module MAPL_SunMod