MAPL_SunOrbitCreate Function

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